Discrete survivor function

Questions
1

I’m following along a tutorial written by Marketer, Bruce Hardie. He doesn’t use Bayesian methods, so I’m trying to convert his MLE estimates into PyMC3 code.

P(churn) = Bernoulli(\theta)
P(T=t) = Geometric(\theta) = \theta(1-\theta)^{t-1}
P(T>t) = S(t) = (1-\theta)^{t-1}

I’m not sure how best to implement this survivor function (as the likelihood) in PyMC3. The two most obvious choices are:

  1. S(t) = (1-\theta)^{t-1} which in PyMC3 might just be
    lik = pm.Bernoulli(1-theta)**(t-1)

  2. S(t) = \frac{Geometric(\theta)}{Bernoulli(\theta)} , which in PyMC3 might be
    lik = pm.Geometric(theta)/Bernoulli(theta)

My two concerns are: (A) I’ve never combined likelihoods in PyMC3 before, so I’m not sure if option 2 would even work. However, (B) I’ve never multiplied or raised a likelihood to a power before and so I’m not sure if this would work either.

Thoughts?

Edit:
I’ve taken inspiration from this post and have used the below model:

with pm.Model() as model:    
    beta_alpha = pm.Uniform('beta_alpha', 0.0001, 5)
    beta_beta = pm.Uniform('beta_beta', 0.0001, 5)
    
    churn = pm.Beta('churn',
                   alpha=beta_alpha,
                   beta=beta_beta)
    renewal = pm.Deterministic('renewal', 1-churn)
    renewals = pm.Geometric('renewals',p=churn, observed=geo)/churn
    
    trace = pm.sample(chains=4)

Whether or not I divide by churn, the posterior estimates are the same. Philosophically, I’m curious as to why this is? And second, is PyMC3 perhaps ignoring this division as it’s not part of what PyMC3 perceives as the likelihood?

And when I try to write a custom likelihood function using pm.potential the sampler goes haywire…

with pm.Model() as model:    
    beta_alpha = pm.Uniform('beta_alpha', 0.0001, 5)
    beta_beta = pm.Uniform('beta_beta', 0.0001, 5)
    
    churn = pm.Beta('churn',
                   alpha=beta_alpha,
                   beta=beta_beta)
    renewal = pm.Deterministic('renewal', 1-churn)

    def likelihood(theta, t):
      return theta**(t-1)
    
    lik = pm.Potential('like', likelihood(theta=renewal, t=geo))
    trace = pm.sample(chains=4)