<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 The aim is to find the parameters of the convolution of two geometric distributions with different probabilities, that is the two probabilities $r$ and $s$ which are the parameters of the two underlying geometric distributions. This problem is related to anactual problem involving the process of some items, so we'll name the parameters to fit that problem; $r$ is the probability that a item is ready for selection, and $s$ is the probability that a ready item is selected to send out. 
The hypothesized distributions are
$$\text{Ready} \sim \text{Geo}(r)$$ 
and
$$\text{Selected} \sim \text{Geo}(s)$$
And the time it takes from an item being started to when an item is sent out is 
$$\text{Time} = \text{Ready} + \text{Selected}.$$
I want to see if PyMC3 can be used to distinguish the two parts of the convolution, and return $r$ and $s$. To start, we going to create a data set of 50000 draws from the $\text{Time}$ distribution to run the PyMC model against. 
Create Fake Data
Step1:  Infer the Readiness and Selection parameters
Now, we declare the Bayesian model we'll use to try to infer $r$ and $s$. 
Try1
We start by instantiating a Model object, and adding $r$ and $s$ as variableds with Uniform prior distributions.
Step2:  I didn't understand this error message at the time, and I moved on to another attempt. But the problem was I had not declared an observed argument for any of the random variables, so PyMC3 had nothing to work with. I also discovered later that the method supplied to DensityDist must be the log probability of the probability distribution, and while I could pass 
python
ready + selected
as I did, the answers provided would not be correlated to the input x, so the inference would have completely wrong. 
Try 2
I had an inckling of the uncorrelated problem above, so I tried using the geometric probability mass function in scipy.stats.
Step3:  But $r$ and $s$ are not numbers; they are random variables, and the pmf mthod cannot accept random variables as an input for p. 
Try 3
Next, I tried subtracting out a geometric distribution, and fitting the residual.
Step4:  The $s$ parameters getting fit with little uncertainty, but it's uncertainty region doesn't extend to cover either parameters true value.
Step5:  For $r$, it never converges a particular value, and it has uncertainty spread from 0 to 1. There is clearly not enough of a "forcing" function using the clip to force $r$ to converge. 
Try 4
My google searching led me to Potentials, so I thought I'd try one. The idea is to reduce the liklihood of random variable values which don't "make sense".
Step6:  The potential here defines not making sense to be value of $ready+selected$ which is not equal to a data point.
Step7:  This is clearly not correct either. The main thing is that liklihood is calculated on the log p scale, so log probabilities are added, not multiplied. The potential I chose therefore doesn't affect much, as adding 0 or 1 is not a large penalty to incur for being "unlikly." I think this led to the results we see here. 
Try 5
Next, I tried to calculate the geometric convolution explicitly for the two parameters, and use that as the liklihood function in a Potential. I think I flipped back and forth a few times with a DensityDist as well. Nothing worked.
Step8:  Try 6
Step9:  Try 7
Step10:  Try 8
<ASSISTANT_TASK:>
Python Code:
# Fake data parameters
r_actual = 1/7
s_actual = 2/17
data_set_size = 50000
# a number of draws from the summed geometric distributions
time_draws = (np.random.geometric(r_actual, size=data_set_size) 
              + np.random.geometric(s_actual, size=data_set_size)
              )
time_draws[:15]
# 0 is not a valid draw from either geometric idstribution
assert 0 not in time_draws
# plot the histogram of the count of each integer
bins = max(time_draws)-min(time_draws)
plt.hist(time_draws, bins=bins);
Explanation: The aim is to find the parameters of the convolution of two geometric distributions with different probabilities, that is the two probabilities $r$ and $s$ which are the parameters of the two underlying geometric distributions. This problem is related to anactual problem involving the process of some items, so we'll name the parameters to fit that problem; $r$ is the probability that a item is ready for selection, and $s$ is the probability that a ready item is selected to send out. 
The hypothesized distributions are
$$\text{Ready} \sim \text{Geo}(r)$$ 
and
$$\text{Selected} \sim \text{Geo}(s)$$
And the time it takes from an item being started to when an item is sent out is 
$$\text{Time} = \text{Ready} + \text{Selected}.$$
I want to see if PyMC3 can be used to distinguish the two parts of the convolution, and return $r$ and $s$. To start, we going to create a data set of 50000 draws from the $\text{Time}$ distribution to run the PyMC model against. 
Create Fake Data
End of explanation
with pm.Model() as try1:
    r = pm.Uniform('r', 0,1)
    s = pm.Uniform('s', 0,1)
# declare the two related geometric distributions
with try1:
    ready = pm.Geometric('ready', p=r)
    selected = pm.Geometric('selected', p=s)
# Define a custom density distribution using the two random variables
with try1:
    time = pm.DensityDist('time', lambda x: ready+selected, )
with try1:
    try:
        trace = pm.sample(10000, step=pm.Metropolis())
    except ValueError as err:
        print(type(err).__name__,':', err)
Explanation: Infer the Readiness and Selection parameters
Now, we declare the Bayesian model we'll use to try to infer $r$ and $s$. 
Try1
We start by instantiating a Model object, and adding $r$ and $s$ as variableds with Uniform prior distributions.
End of explanation
with pm.Model() as try2:
    r = pm.Uniform('r', 0,1)
    s = pm.Uniform('s', 0,1)
    try:
        time = pm.DensityDist('time', lambda x: ss.geom.pmf(x, p=r) + ss.geom.pmf(x, p=s))
    except TypeError as err:
        print(type(err).__name__,':', err)
Explanation: I didn't understand this error message at the time, and I moved on to another attempt. But the problem was I had not declared an observed argument for any of the random variables, so PyMC3 had nothing to work with. I also discovered later that the method supplied to DensityDist must be the log probability of the probability distribution, and while I could pass 
python
ready + selected
as I did, the answers provided would not be correlated to the input x, so the inference would have completely wrong. 
Try 2
I had an inckling of the uncorrelated problem above, so I tried using the geometric probability mass function in scipy.stats.
End of explanation
with pm.Model() as try3:
    r = pm.Uniform('r', 0,1, testval=1/10)
    s = pm.Uniform('s', 0,1, testval=1/10)
with try3:
    ready = pm.Geometric('ready', p=r, )
    
    # create residual data
    adjusted_data = time_draws - ready
    
    # force the residuals to be non-negative by clipping the numbers in a particular interval
    pm.math.clip(adjusted_data, 0, np.inf)
    
    # fit the 2nd geometric to the residuals
    selected = pm.Geometric('selected', p=s, observed=adjusted_data)
with try3:
    try:
        trace = pm.sample(10000, step=pm.Metropolis())
    except ValueError as err:
        print(type(err).__name__,':', err)
pm.traceplot(trace=trace[:2000],);
Explanation: But $r$ and $s$ are not numbers; they are random variables, and the pmf mthod cannot accept random variables as an input for p. 
Try 3
Next, I tried subtracting out a geometric distribution, and fitting the residual.
End of explanation
r_actual, s_actual
Explanation: The $s$ parameters getting fit with little uncertainty, but it's uncertainty region doesn't extend to cover either parameters true value.
End of explanation
with pm.Model() as try4:
    r = pm.Uniform('r', 0,1, testval=1/10)
    s = pm.Uniform('s', 0,1, testval=1/10)
    ready = pm.Geometric('ready', p=r, shape=data_set_size)
    selected = pm.Geometric('selected', p=s, shape=data_set_size)
Explanation: For $r$, it never converges a particular value, and it has uncertainty spread from 0 to 1. There is clearly not enough of a "forcing" function using the clip to force $r$ to converge. 
Try 4
My google searching led me to Potentials, so I thought I'd try one. The idea is to reduce the liklihood of random variable values which don't "make sense".
End of explanation
with try4:    
    pot = pm.Potential('pot', pm.math.sum(pm.math.switch( pm.math.eq(time_draws, ready+selected), 1, 0))/data_set_size)
with try4:
    try:
        trace = pm.sample(10000, step=pm.Metropolis())
    except ValueError as err:
        print(type(err).__name__,':', err)
pm.traceplot(trace=trace,varnames=['r', 's',]);
Explanation: The potential here defines not making sense to be value of $ready+selected$ which is not equal to a data point.
End of explanation
def geo_convolution_likelihood(p,q):
    probabilities = [p*q*pm.math.sum( [((1-p)**inner) * ((1-q)**(outer-inner))
                                     for inner in range(outer+1)]) 
                     for outer in time_draws]
    intersection_prob = 1
    for prob in probabilities:
        intersection_prob *= prob
    return intersection_prob
with pm.Model() as try5:
    r = pm.Uniform('r', 0,1, testval=1/10)
    s = pm.Uniform('s', 0,1, testval=1/10)
with try5:
    pm.Potential('TotalDays', geo_convolution_likelihood(r,s))
with try5:
    try:
        trace = pm.sample(10000, step=pm.Metropolis())
    except ValueError as err:
        print(type(err).__name__,':', err)
pm.traceplot(trace=trace,varnames=['r', 's',]);
Explanation: This is clearly not correct either. The main thing is that liklihood is calculated on the log p scale, so log probabilities are added, not multiplied. The potential I chose therefore doesn't affect much, as adding 0 or 1 is not a large penalty to incur for being "unlikly." I think this led to the results we see here. 
Try 5
Next, I tried to calculate the geometric convolution explicitly for the two parameters, and use that as the liklihood function in a Potential. I think I flipped back and forth a few times with a DensityDist as well. Nothing worked.
End of explanation
with pm.Model() as try6:
    r = pm.Uniform('r', 0,1, testval=1/7)
    s = pm.Uniform('s', 0,1, testval=1/17)
with try6:
    ready = pm.Geometric('ready', p=r, shape=data_set_size)
    selected = pm.Geometric('selected', p=s, shape=data_set_size)
    
    selected_draws = time_draws - ready
    pm.math.clip(selected_draws, 0, np.inf)
    selected_obs = pm.Geometric('selected_obs', p=s, observed=selected_draws)
    
    ready_draws = time_draws - selected
    pm.math.clip(ready_draws, 0, np.inf)
    ready_obs = pm.Geometric('ready_obs', p=r, observed=ready_draws)
with try6:
    try:
        trace = pm.sample(10000, step=pm.Metropolis())
    except ValueError as err:
        print(type(err).__name__,':', err)
pm.traceplot(trace=trace[5000:], varnames=['r', 's']);
math.sqrt((1/17)* (1/7))
Explanation: Try 6
End of explanation
with pm.Model() as try7:
    r = pm.Uniform('r', 0,1, testval=1/7)
    s = pm.Uniform('s', 0,1, testval=1/17)
with try7:
    ready = pm.Geometric('ready', p=r, shape=data_set_size)
    selected = pm.Geometric('selected', p=s, shape=data_set_size)
    
    th try7:
    try:
        trace = pm.sample(10000, step=pm.Metropolis())
    except ValueError as err:
        print(type(err).__name__,':', err)otal = ready + selected - time_draws
    
    pm.Potential('TotalDiff', 1/pm.math.dot(total, total))
with try7:
    try:
        trace = pm.sample(10000, step=pm.Metropolis())
    except ValueError as err:
        print(type(err).__name__,':', err)
pm.traceplot(trace=trace[5000:], varnames=['r', 's']);
Explanation: Try 7
End of explanation
with pm.Model() as try8:
    r = pm.Uniform('r', 0,1, testval=1/7)
    s = pm.Uniform('s', 0,1, testval=1/17)
    
    all_convolutions = [r*s*pm.math.sum([r**i * s**(conv_num-2-i) for i in range(conv_num-2)])
                        for conv_num in range(1,25)
                        ]
    def log_density(x):
        convol = pm.Deterministic('convol', all_convolutions[x])
        return pm.math.log(r) + pm.math.log(s) + pm.math.log(convol)
    
    sum_obs = pm.DensityDist('sum_obs', log_density, observed=time_draws)
with try8:
    try:
        trace = pm.sample(10000, step=pm.Metropolis())
    except ValueError as err:
        print(type(err).__name__,':', err)
pm.traceplot(trace=trace[:2000],);
# Example starting code in this cell
n1 = 1000
n = n1*2
r = 1/17
q = 1/3
data1 = np.random.geometric(r,size=n1)
data2 = np.random.geometric(q,size=n1)
data = data1 + data2
datatt = shared(data)
with pm.Model() as model:
    
    t = pm.DiscreteUniform( "t", 1, data-1, shape=len(data)) 
    
    r = pm.Uniform('r',0,1, testval=1/17)
    s = pm.Uniform('s',0,1, testval=1/3)
    
    Ready = pm.Geometric('Ready', r, observed=t)
    Selected = pm.Geometric('Selected', s, observed=data-t)
    
    ordering = pm.Potential('ordering',pm.math.switch(r >= s, -10000, 0))
with model:
    step1 = pm.NUTS([r,s, ordering])
    step2 = pm.Metropolis([t])
    trace = pm.sample(20000, step=[step1, step2] )
pm.traceplot(trace[1000:],);
trace['t']
ind = len(data)*[0] + len(data)*[1]
with pm.Model() as model:
    t = pm.DiscreteUniform( "t", 1, data-1, shape=len(data)) 
    
    p = pm.Uniform('p',0,1,shape=2, testval=np.array([r,q]))
    
  
    testmodel = pm.Deterministic('testmodel',tt.concatenate([t,datatt-t]))
    
    process = pm.Geometric('process',p=p[ind],observed=testmodel)
    
with model:
    
    step1 = pm.Metropolis([t])
    
    step2 = pm.NUTS([p])
    step3 = pm.Metropolis([testmodel, process])
    trace = pm.sample(50000, [step1, step2, step3])
    pm.traceplot(trace[1000:][::10], )
trace['testmodel']
trace['t']
with pm.Model() as model:
    t = pm.DiscreteUniform( "t", 1 , data, shape=len(data)) 
    
    p = pm.Uniform('p',0,1,shape=2, testval=np.array([r,q]))
    
    prob_diff = pm.Deterministic('diff', pm.math.switch(p[0] + 0.20 > p[1], -10000, 0))
    sepe = pm.Potential('sepe', var=prob_diff)
    
    ind = len(data)*[0] + len(data)*[1]
  
    testmodel = pm.Deterministic('testmodel',rateFunc(t))
    
    process = pm.Geometric('process',p=p[ind],observed=testmodel)
  
with model:
    trace = pm.sample(50000, pm.Metropolis())
pm.traceplot(trace[1000:][::10], ['p']);
with pm.Model() as model:
    t = pm.DiscreteUniform( "t", 1 , data, shape=len(data)) 
    
    p = pm.Uniform('p',np.array([0,0.3]),np.array([0.1,1]),shape=2, testval=np.array([r,q]))
    prob_diff = pm.Deterministic('diff', pm.math.switch(p[0] + 0.20 > p[1], -10000, 0))
    sepe = pm.Potential('sepe', var=prob_diff)
    
    ind = len(data)*[0] + len(data)*[1]
  
    testmodel = pm.Deterministic('testmodel',rateFunc(t))
    
    process = pm.Geometric('process',p=p[ind],observed=testmodel)
  
with model:
    trace = pm.sample(50000, pm.Metropolis())
pm.traceplot(trace[1000:][::10], ['p']);
Explanation: Try 8
End of explanation
<END_TASK>