# Bayesian Discrete Survival Analysis

I've been reading over Bruce Hardie's tutorial on survival analysis of customers in a discrete-time subscription model. Two design choices he delineates early on include the Geometric PMF: $$P(T=t|\theta) = \theta^{1}*(1-\theta)^{t-1}$$. And the Survival function: $$P(T>t|\theta) = (1-\theta)^{t}$$.

The fictitious data he provides:

t, customers, losses
_,__________,_______
0, 1000,     0
1, 631,      369
2, 468,      163
3, 382,      86
4, 326,      56


And so the likelihood is function he provides (See page 14 of link for more info) can be expressed as $$Survival(t|theta)^{survivors} * \prod_{1}^{t} Geometric(t|theta)^{Losses[t]}$$.

His MLE estimate of $$\theta = 0.27$$, however, he notes that there is no variance around this estimate and so improves the model by incorporating the Beta PDF to handle "heterogeneity around theta". He's still using MLE methods so he doesn't introduce the use of $$Beta(\theta|\alpha, \beta)$$ as a prior. Rather, he describes the new likelihood function as a "Geometri-Beta" distribution. I would argue that he is indeed adding a Beta prior to his likelihood and arriving at a MAP estimate, but that's my naive two-cents.

I implemented a model in PyMC3, needing to specify the custom likelihood function as described.

s = 326
L = [369, 163, 86, 56]

with pm.Model() as model:
alpha = pm.HalfNormal('alpha',sigma=1)
beta = pm.HalfNormal('beta',sigma=1)

churn = pm.Beta('churn',
alpha=alpha,
beta=beta)
renewal = pm.Deterministic('renewal', 1-churn)

def likelihood(theta):
def logp_(s,L):
"""
@s: number of survivors whose death is censored
@L: array of losses over time (discrete sequence)
@target: summation of log loss (output)
"""
L1 = list(L)
target = 0
# step1: compute geometric(t|theta) for l in Losses
for i in range(len(L1)):
t = i+1
l = L1[i]

log_p_churn = np.log(theta)
log_p_renew = (t-1) * np.log(1-theta)
log_geometric = log_p_churn + log_p_renew
log_loss = l * log_geometric
target += log_loss

# step 2: compute survival(s|theta)
t = len(L1)

log_survive = t * np.log(1-theta)
log_loss = s * log_survive
target += log_loss

return target
return logp_

lik = pm.DensityDist('likelihood', logp=likelihood(churn), observed=dict(s=s, L=L))
trace = pm.sample(chains=4)


The author demonstrates that this is what the distribution over churn probabilities, $$\theta$$ should look like (See page 24.)

However, my posterior estimate of $$\theta$$ looks like this, with mean parameter values of $$\alpha=0.82, \beta=1.08, churn=0.27, renewal=1-churn=0.72$$

I decided to plot what churn should look like given $$Beta(\alpha=0.82,\beta=1.08)$$ and I was surprised to find that it did not resemble the posterior plot of turn much at all!

So, putting two and two together; my posterior estimate of churn was very close to the author's MLE estimate (0.27). And so, it seemed that my custom likelihood function needs to incorporate the prior. In other words, PyMC3 doesn't automatically augment likelihoods by priors defined in the model.

So I decided to rewrite the likelihood function so that the I could account for the log pdf value of $$Beta(\theta|\alpha, \beta)$$. Note, the betaln function returns gammaln(x) + gammaln(y) - gammaln(x + y).

from pymc3.distributions.dist_math import betaln

with pm.Model() as model:
alpha = pm.HalfNormal('alpha',sigma=0.5)
beta = pm.HalfNormal('beta',sigma=0.5)

churn = pm.Beta('churn',
alpha=alpha,
beta=beta)
renewal = pm.Deterministic('renewal', 1-churn)

def likelihood(theta,a,b):

def logbetapdf(theta, a,b):
return (a-1)*tt.log(theta) + (b-1)*tt.log(1-theta) - betaln(a,b)

def logp_(s,L,):
"""
@s: number of survivors whose death is censored
@L: array of losses over time (discrete sequence)
@target: summation of log loss (output)
"""
L1 = list(L)
target = 0
# step1: compute geometric(t|theta) for l in Losses
for i in range(len(L1)):
t = i+1
l = L1[i]

log_p_churn = np.log(theta)
log_p_renew = (t-1) * np.log(1-theta)
log_geometric = log_p_churn + log_p_renew + logbetapdf(theta,a,b)
log_loss = l * log_geometric
target += log_loss

# step 2: compute survival(s|theta)
t = len(L1)

log_survive = t * np.log(1-theta) + logbetapdf(theta,a,b)
log_loss = s * log_survive
target += log_loss

return target
return logp_

lik = pm.DensityDist('likelihood', logp=likelihood(churn, alpha, beta), observed=dict(s=s, L=L))
trace = pm.sample(chains=4)


So, after all this, I have a new posterior plot, but this time churn is distributed according to $$Beta(\theta|\alpha,\beta)$$. .

However, these posterior predictions are quite different from what the shape that the author inferred... which is troubling.

A few questions:

1. I've used the half normal as I wanted to keep $$\alpha$$ and $$\beta$$ close-ish to 0. this didn't really work that well. Should I be using different priors?
2. Is it my priors that are causing my estimates to differ so much from the author's Beta-Geomeric mixture MLE estimates?
3. Is there a bug in my code from the latter model? I went back and forth deciding where to account for the prior via logbetapdf(theta,a,b). Ultimately, I chose to account for it along with the geometric and survivor functions, respectively, and scale these outputs as a function of the specific numbers of people who churned/survived respectively. Was this the right choice?
4. Anything else you'd like to add?

Thank you!!