PyMC observed data for a sum of random variables

Question

I'm trying to infer models parameters with PyMC. In particular the observed data is modeled as a sum of two different random variables: a negative binomial and a poisson.

In PyMC, an algebraic composition of random variables is described by a "deterministic" object. Is it possible to assign the observed data to this deterministic object?

If not possible, we still know that the PDF of the sum is the convolution the PDF of the components. Is there any trick to compute this convolution efficiently?

Answer 1

It is not possible to make a deterministic node observed in PyMC2, but you can achieve an equivalent model by making one part of your convolution a latent variable. Here is a small example:

def model(values):
    # priors for model parameters
    mu_A = pm.Exponential('mu_A', beta=1, value=1)
    alpha_A = pm.Exponential('alpha_A', beta=1, value=1)
    mu_B_minus_A = pm.Uninformative('mu_B_minus_A', value=1)

    # latent variable for negative binomial
    A = pm.NegativeBinomial('A', mu=mu_A, alpha=alpha_A, value=0)

    # observed variable for conditional poisson
    B = pm.Poisson('B', mu=mu_B_minus_A+A, value=values, observed=True)

    return locals()

Here is a notebook that tests it out. It seems like it will be tough to fit without some additional information on the model parameters. Perhaps there is a clever way to calculate or approximate the convolution of a NB and a Poisson that you could use as a custom observed stochastic instead.