Hello,
I am working on converting the Fader & Hardie shifted Beta-Geometric model into PyMC3 and have run into a hitch.
I’m using the data in the appendix to and then attempting to write the custom likelihood function. I’ve seen a post where this was done having the probability recursively defined, but I would like to calculate it directly (following Eq (5) and (6) in the paper). I wrote a test code to calculate the log-likelihood and verify that I get the expected result from the paper (-2115.545…) which I do:
import numpy as np
import pymc3 as pm
import arviz as az
from scipy.special import beta as beta_f
survived = np.array([869, 743, 653, 593, 551, 517, 491])
churned = np.array([131, 126, 90, 60, 42, 34, 26])
n_obs = len(survived)
def logp_test(survived, churned, alpha, beta):
ln_surv_prob = np.log(beta_f(alpha, beta + n_obs)/beta_f(alpha, beta))
ln_prob_vec = np.ones(n_obs)
for i in range(1, n_obs + 1):
ln_prob_vec[i - 1] = beta_f(alpha + 1, beta + i - 1)/beta_f(alpha, beta)
ln_prob_vec = np.log(ln_prob_vec)
return np.dot(churned, ln_prob_vec) + survived[-1] * ln_surv_prob
My issue comes when I try and run this within the context of a PyMC3 model. I initially tried to insert this into a model and was basically told that the scipy.special.beta function cannot be run in this way. I dug around for a while before ending up on the idea of using as_op. Following this example on “as_op”, I wrote the following:
import numpy as np
import pymc3 as pm
import arviz as az
import theano.tensor as tt
from math import gamma as gamma_f
from theano.compile.ops import as_op
survived = np.array([869, 743, 653, 593, 551, 517, 491])
churned = np.array([131, 126, 90, 60, 42, 34, 26])
n_obs = len(survived)
@as_op(itypes = [tt.dscalar, tt.dscalar], otypes = [tt.dscalar])
def beta_func(a, b):
return (gamma_f(a) * gamma_f(b) / gamma_f(a + b))
with pm.Model() as sbg_model:
alpha = pm.Uniform('alpha', 0.0001, 10) # MLE est was 0.668
beta = pm.Uniform('beta', 0.0001, 10) # MLE est was 3.806
def logp(survived, churned, alpha = alpha, beta = beta):
# Calculate the final surivival probability (eq. 6)
ln_surv_prob = np.log(beta_func(alpha, beta + n_obs)/beta_func(alpha, beta))
# Find the probability of churn for all values prior
ln_prob_vec = np.ones(n_obs)
for i in range(1, n_obs + 1):
ln_prob_vec[i - 1] = beta_func(alpha + 1, beta + i - 1)/beta_func(alpha, beta)
ln_prob_vec = np.log(ln_prob_vec)
return pm.math.dot(churned, ln_prob_vec) + survived[-1] * ln_surv_prob
likelihood = pm.DensityDist('likelihood', logp, observed = {'survived': survived, 'churned': churned})
sbg_trace = pm.sample(200) # just a sanity check to see if it even runs...
However, when I run it, I get:
TypeErrorTraceback (most recent call last)
TypeError: float() argument must be a string or a number, not 'TensorVariable'
The above exception was the direct cause of the following exception:
ValueErrorTraceback (most recent call last)
<ipython-input-105-e24d51fa3067> in <module>
18 return pm.math.dot(churned, ln_prob_vec) + survived[-1] * ln_surv_prob
19
---> 20 likelihood = pm.DensityDist('likelihood', logp, observed = {'survived': survived, 'churned': churned})
21 sbg_trace = pm.sample(200) # just a sanity check to see if it even runs...
/usr/local/lib/python3.8/site-packages/pymc3/distributions/distribution.py in __new__(cls, name, *args, **kwargs)
120 else:
121 dist = cls.dist(*args, **kwargs)
--> 122 return model.Var(name, dist, data, total_size, dims=dims)
123
124 def __getnewargs__(self):
/usr/local/lib/python3.8/site-packages/pymc3/model.py in Var(self, name, dist, data, total_size, dims)
1160 elif isinstance(data, dict):
1161 with self:
-> 1162 var = MultiObservedRV(
1163 name=name,
1164 data=data,
/usr/local/lib/python3.8/site-packages/pymc3/model.py in __init__(self, name, data, distribution, total_size, model)
1871 datum.missing_values for datum in self.data.values() if datum.missing_values is not None
1872 ]
-> 1873 self.logp_elemwiset = distribution.logp(**self.data)
1874 # The logp might need scaling in minibatches.
1875 # This is done in `Factor`.
<ipython-input-105-e24d51fa3067> in logp(survived, churned, alpha, beta)
13 ln_prob_vec = np.ones(n_obs)
14 for i in range(1, n_obs + 1):
---> 15 ln_prob_vec[i - 1] = beta_func(alpha + 1, beta + i - 1)/beta_func(alpha, beta)
16 ln_prob_vec = np.log(ln_prob_vec)
17
ValueError: setting an array element with a sequence.
I’m quite frankly stuck. Can anyone point me in a good direction? I am running pymc3 = 3.11.2
Thank you in advance!