PyMC3中的离散生存函数(自定义似然)

我正试图跟随一个营销教程。本教程使用Frequentist/MLE方法;我喜欢PyMC3并决定使用它。本教程的作者使用了一个生存函数,即

S(t|churn_rate) = (1-churn_rate)**(t-1) S(t|churn_rate) = churn_rate*(1-churn_rate)**(t-1) .

PyMC3内置了几何分布,所以我的问题不在这里。而是找到一种把生存函数写成可能性的方法。

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)

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.

likelihood log_likelihood

我这边有几个疑点:

1. 目前还不清楚 pm.Potential 或 pm.DensityDist . 社会似乎认为 这是一个更好的选择。
2. 我正在将一个名为geo的数组传递到 . 也许它期望的是一个标量,不太确定数组的组成。。。

资料来源:

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

Blackbox likelihood example