Shifted Beta Geometric (sBG) distribution in pymc3

Questions
2

Your implementation doesnt seems correct. For example, checking the log-likelihood it does not match the one in the post:

retention.logp(model.test_point)
array(-16884.090589033225)

The problem is that you need to rewrite the numpy part into theano. Something like this should work (sample with NUTS, which is :100:):

with pm.Model() as model:
    alpha = pm.Uniform('alpha', 0.00001, 1000.0, testval=1)
    beta = pm.Uniform('beta', 0.00001, 1000.0, testval=1)

    p = [0., alpha / (alpha + beta)]
    s = [0., 1 - p[1]]
    for t in range(2, n):
        pt = ((beta + t - 2) / (alpha + beta + t - 1)) * p[t-1]
        p.append(pt)
        s.append(s[t-1] - p[t])
    p = tt.stack(p)
    s = tt.stack(s)
    
    def logp(value):
        active = value[0,:]
        lost = value[1,:]

        # Those who've churned along the way...
        died = tt.mul(tt.log(p[1:]), lost[1:])

        # and those still active in last period
        still_active = tt.log(s[-1]) * active[-1]
        return  tt.sum(died) + still_active

    retention = pm.DensityDist('retention', logp, observed=data)
    trace = pm.sample(1000)
pm.traceplot(trace)