I am trying to build a model to estimate the parameters of a Compound Poisson process with a Binomial distribution as the underlying distribution. The process Y(t) is defined as: Y(t) = \sum_{m=1}^{N(t)}{Z_m}, where Z_m ~ i.i.d. Binom(n0, p0) and N(t) ~ Poisson(\lambda*t). From this definition, one can see that the number of Z_m’s is variable and has to be dynamically set based on N(t). I have come up with this basic framework, but am having trouble with dynamically setting the shape of pm.Binomial
.
import numpy as np
import pymc as pm
import matplotlib.pyplot as plt
def sampling(data):
"""
Markov Chain Monte Carlo estimation of the Compound Poisson process
"""
with pm.Model() as cp_model:
# Priors
rate = pm.Gamma('rate', alpha=1, beta=1)
n0 = pm.Gamma('n0', alpha=5, beta=1)
p0 = pm.Beta('p0', alpha=2, beta=2)
n_t = pm.Poisson.dist(mu=rate)
trunc_n_t = pm.Truncated('trunc_n_t', n_t, lower=1)
n_t_draw = pm.draw(trunc_n_t)
print(n_t_draw)
z_m = pm.Binomial('z_m', n=n0, p=p0, shape=(n_t_draw,))
# Likelihood
y_t = pm.Deterministic('y_t', pm.math.sum(z_m))
likelihood = pm.Normal('likelihood', mu=y_t, observed=data)
with cp_model:
trace = pm.sample(
draws=30000,
chains=4,
tune=1000,
step=pm.Metropolis(),
cores=2
)
res = pm.summary(trace)
print(res)
pm.plot_trace(trace, rug=False, divergences=None)
plt.tight_layout()
plt.show()
return res
I tried to use pm.draw
to set the shape of pm.Binomial
, but from the summary table and plot trace, the shape seems to be set on model definition instead of dynamically changed during training. Any help would be appreciated, thanks!