Hi all!
First thanks for this very helpful forum.
I’m having a hard time tying to model the occurrence of some events.
When an event occurs at time t (following lets say a gamma distribution of average u1), this event will can induce 0->n events occurs (following a poisson distribution, probability p) at the time t + t_diff (t_diff following another gamma, of mean u_diff).
My data is a collection of the time of those events, I don’t know if observation comes from the first or the secondary process.
I’m mostly interested in estimating u1, p, and u_diff from my data.
I’m quite new to Bayesian statistic, my data looks like a bimodal mixture, but to me it is not the right way to look at it because the second distribution depends on the first one, and the formulation above is closer to the actual process of my subject of observation.
Here is a trial to implement this model:
import numpy as np
import pymc as pm
# generate some data
n = 2000
# first process of events
t1 = np.random.gamma(144, 0.83 size=(n,)) # u: 120, sigm: 10
p = np.random.poisson(0.5, n)
# second process of events
t_tmp = np.repeat(t1, p)
t2 = t_tmp + np.random.gamma(225, 0.27, t_tmp.shape) # u: 60, sigma: 4
t = np.hstack((t1, t2))
with pm.Model() as model_g:
μ1 = pm.Uniform('μ1', lower=50, upper=300)
σ = pm.HalfNormal('σ', sd=10)
p = pm.Poisson('p', 0.5) # Don't know how to include this in y2
μ_diff = pm.Uniform('μ_diff', mu=30, upper=70)
y1 = pm.Gamma('y', mu=μ1, sigma=σ)
y2 = pm.Gamma('y', mu=μ1+μ_diff, sigma=σ) # not exactly the formulation of the problem, but good enough
# ... ?
I read on other post that there is no way to “combine” distributions in pymc. I read also about pm.Dirichlet, but it seems it can’t handle the kind of constrains my modeling needs.
In short, I’m not sure If/how I can join y1 and y2 to fit on my data, or said differently, model my data considering the two processes with their specific parameters.
Sorry if the question is too obvious, and thanks for your advices !
Best