I am relatively new to pymc, so this problem might have a very simple solution that I am not seeing:
For my current project I am trying to infer some physical parameters from a set of time-resolved measurements. Essentially the output of the measurement is a histogram as the one shown in the figure, where a number of counts is binned into time-segments:
You can see that a) the data follows an almost-exponential distribution and b) the relative error is larger at lower count numbers (obviously, I guess).
I have a physical model that I have translated into pytensor, so that I can simulate the same decay from a set of parameters. Now I would like to use pymc to infer the set of parameters that best fit the data.
However, if I just use
Y_obs = pm.Normal('Y_obs', mu = simulation, sigma = sigma, observed = data)
the lower counts at later times are essentially ignored, because they don’t change the overall likelihood too much anymore. Similarly I have tried
Y_obs_log = pm.Normal('Y_obs_log', mu = at.log(simulation), sigma = sigma, observed = np.log(data))
but then the relative change in sigma with time seems to become an issue. Lastly, I have also tried
Y_obs_poiss = pm.Poisson('Y_obs_poiss', mu = simulation/simulation.sum(), observed = data/data.sum() )
w = data Y_obs_poiss2 = pm.Potential('Y_obs_poiss2 ', w*pm.logp(rv=pm.Poisson.dist(simulation/simulation.sum()), value=data/data.sum()))
but neither have resulted in anything useful so far.
I would appreciate any input of how to combine the continuous model I have with the discrete, binned data in a ‘proper’ way to infer the parameters. Maybe there is something obvious I am missing here.
Thank you very much already!