# Observed deterministic variables in MCMC

I need to model a measurement of an "exponential decay" i.e. I have a histogram of counts $$Y$$ over an array of (time-) intervalls. I want to use MCMC to infer parameters ($$A_1,\lambda_1,A_2,\lambda_2,\lambda_{bgd}$$) of a model $$Y = A_1\mathrm{Poisson}(\lambda_1)+A_2\mathrm{Poisson}(\lambda_2)+\mathrm{Poisson}(\lambda_{bgd})$$ that describes a decay with two two poisson variables, two associated amplitudes and a poisson background.

Currently I am working with pymc3, where it appears that observed variables (e.g. $$Y$$) cannot be deterministic variables at the same time. As I understand the problem, no likelihood can be associated to a determinitisc variable for sampling (no likelihood for the sum/product of the stochastic variables $$A_i,\lambda_i$$).

So far I could think of two approaches to adress the problem:

• Try to find an explicit expression for the likelihood of parameters $$A_1,...,\lambda_{bgd}$$ given the data. Similar to here (under lighthouse problem) http://www.awebb.info/blog/observing_functions
• Model $$Y$$ itself as a stochastic variable (e.g. normal) with some fixed variance and mean at $$A_1\mathrm{Poisson}(\lambda_1)+A_2\mathrm{Poisson}(\lambda_2)+\mathrm{Poisson}(\lambda_{bgd})$$ as suggested here: https://discourse.pymc.io/t/observed-deterministic/483

The first seems complicated, the latter seems quite hackish an I do not know if it would be rigorous.

Can you point me towards a way to solve this problem of observed varibale is deterministic variable? (I do not mean this as a purely technical question for pymc, but more conceptually how to approach such a problem.)

My original comment misunderstood the nature of your problem. Indeed, it is correct that PyMC cannot model systems where the observed variable is deterministic. That was a bit of a surprise to me, but the article that you link provides a good reason why that is the case.

Now, for your problem: Given that we cannot directly observe the sum-of-products $$Y$$ as you have given it, a good solution would be to find an expression for $$Y$$ as a single random variable parametrised by ($$A_1,\lambda_1,A_2,\lambda_2,\lambda_{bgd}$$). We could then plug this into PyMC as originally intended (but now the observed variable would not be deterministic).

Fortunately, that's actually pretty easy: multiplying a Poisson random variable $$Poisson(\lambda)$$ by a constant $$A$$ results in the Poisson random variable $$Poisson(\lambda A)$$. Similarly, adding two Poisson random variables with lambdas $$\lambda_1$$ and $$\lambda_2$$ results in the Poisson random variable $$Poisson(\lambda_1 + \lambda_2)$$. So that means your random variable $$Y$$ can simply be modelled as:

$$Y \sim Poisson(A_1 \lambda_1 + A_2 \lambda_2 + \lambda_3)$$

Which of course brings us to the bad news - if we were to estimate the lambda of this random variable $$Y$$ (either using PyMC or, now that the distribution is much simpler, a more direct method), there is clearly no way that we can isolate the three separate influences of the underlying Poisson distributions.

In other words, your experiment does not give us enough information to surmise the underlying distributions because any linear combination of Poisson random variables is itself just a Poisson random variable.

# Edit

multiplying a Poisson random variable $$Poisson(\lambda)$$ by a constant $$A$$ results in the Poisson random variable $$Poisson(\lambda A)$$

Oops, that's clearly wrong. Per this answer, the distribution of $$Z = \alpha X$$ (where $$X \sim Poisson(\lambda)$$) actually has the PMF:

$$P(Z=z) = \frac{\lambda^{z / \alpha}}{(z / \alpha)!}e^{-\lambda}$$

when $$z / \alpha$$ is a non-negative integer and $$0$$ otherwise.

Which suggests to me that this problem is even more impossible if $$A_1$$ and $$A_2$$ are continuous. If they are continuous, the likelihood will be zero for all samples (since the probability that $$z / \alpha$$ is an integer is zero for continous $$\alpha$$). If they are discrete, then it may be possible to solve this problem but you will probably need to implement your own custom distribution in PyMC, using the PMF given above.

This problem becomes trivial if $$A_1 = A_2 = 1$$.