!date

import numpy as np, pymc as pm, seaborn as sns
%matplotlib inline
np.random.seed(12345) # for reproducibility

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()

n = 10
values = pm.rnegative_binomial(mu=2, alpha=10, size=n) + pm.rpoisson(mu=1, size=n)
m = pm.MCMC(model(values))

m.use_step_method(pm.AdaptiveMetropolis, [m.mu_A, m.alpha_A, m.mu_B_minus_A])
m.use_step_method(pm.DiscreteMetropolis, m.A)

m.sample(100000, 50000, 50)
print
pm.Matplot.plot(m)