Discrete survival function (custom likelihood) in PyMC3

Question

I'm trying to follow along with a Marketing tutorial. The tutorial uses Frequentist/MLE methods; I'm fond of PyMC3 and decided to use that. The author of the tutorial uses a Survival function, which is

S(t|churn_rate) = (1-churn_rate)**(t-1) This is in contrast to the Geometric distribution, which simply adds one more term to the above: S(t|churn_rate) = churn_rate*(1-churn_rate)**(t-1).

PyMC3 has the geometric distribution built in, so my problem isn't there. But rather finding a way to write the survival function as a likelihood.

import arviz as az
import pymc3 as pm
import numpy as np
from pipe import traverse
wins =   [1000, 631, 468, 382, 326]
geo = [[idx+1 for i in range(n)] for idx,n in enumerate(wins)]
geo = np.array(list(geo | traverse)) #flattens the array

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 log_likelihood(theta, t):
      return (t-1)*np.log(theta)
    
    lik = pm.Potential('like', log_likelihood(theta=renewal, t=geo))
    trace = pm.sample(chains=4)

Unfortunately, the sampler has gone haywire...

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Sequential sampling (4 chains in 1 job)
NUTS: [churn, beta_beta, beta_alpha]

 100.00% [2000/2000 00:02<00:00 Sampling chain 0, 821 divergences]

 100.00% [2000/2000 00:03<00:00 Sampling chain 1, 562 divergences]

 100.00% [2000/2000 00:02<00:00 Sampling chain 2, 628 divergences]

 100.00% [2000/2000 00:03<00:00 Sampling chain 3, 364 divergences]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 12 seconds.
There were 822 divergences after tuning. Increase `target_accept` or reparameterize.
The acceptance probability does not match the target. It is 0.498263496658037, but should be close to 0.8. Try to increase the number of tuning steps.
There were 1385 divergences after tuning. Increase `target_accept` or reparameterize.
There were 2013 divergences after tuning. Increase `target_accept` or reparameterize.
The acceptance probability does not match the target. It is 0.6553072990104106, but should be close to 0.8. Try to increase the number of tuning steps.
There were 2378 divergences after tuning. Increase `target_accept` or reparameterize.
The estimated number of effective samples is smaller than 200 for some parameters.

I had previously just written a likelihood function, not a log_likelihood function but the sampler wasn't happy with this either.

A couple of suspicions on my end:

1. It's unclear if pm.Potential or pm.DensityDist. The SO community seems to think that pm.Potential is a better choice.
2. I'm passing an array called geo into log_likelihood. Perhaps it's expecting a scalar and isn't quite to sure what to make of an array...

Sources:

http://www.brucehardie.com/papers/021/sbg_2006-05-30.pdf

Blackbox likelihood example