@jessegrabowski, many thanks for checking! I really like your solution.
However, I tried to run a similar code, but I’ve got many divergences 20221016-pymc-question.ipynb · GitHub
Just to confirm, is it not about priors or something like that?
with pm.Model() as model:
# Data
T = tt.as_tensor_variable(counts.shape[0], dtype="int64")
D = tt.as_tensor_variable(counts.shape[1] - 1, dtype="int64")
D_grid = at.arange(0, D + 1)
counts_matrix = tt.as_tensor_variable(counts, dtype="int64")
# Priors
μ = pm.HalfNormal('μ', 5.0, initval=2.0)
σ = pm.HalfNormal('σ', 5.0, initval=2.0)
# Likelihood
### Parameters of the lognormal distribution
σlog = tt.sqrt(tt.log((σ / μ)**2 + 1))
μlog = tt.log(μ) - σlog**2 / 2
log_norm = pm.Lognormal('log_norm', mu = μlog, sigma = σlog)
def cdf(x):
return tt.exp(pm.logcdf(log_norm, x))
π, _ = ae.scan(lambda t: cdf(t + 0.5) - cdf(t - 0.5), sequences = [D_grid])
pm.Deterministic('p', π)
log_like = pm.Multinomial('obs',
n=counts_matrix.sum(axis=1),
p = π / π.sum(),
observed=counts_matrix)
trace = pm.sample(tune=3000, chains=4, init='advi')