!date

import matplotlib.pylab as plt
%matplotlib inline

import pymc3 as pm # fork of pymc to allow importing pymc2 as well, when necessary

# I prefer to use pm.*, but this keeps consistency with the SO question
from pymc3 import Model, Uniform, Normal, Poisson, Metropolis, traceplot 
from pymc3 import sample


# Simulate data
import scipy.stats

totalRates = scipy.stats.norm(loc=120, scale=20).rvs(size=10000)
totalCounts = scipy.stats.poisson.rvs(mu=totalRates) 
successRate = 0.1*totalRates 
successCounts = scipy.stats.poisson.rvs(mu=successRate) 

with Model() as success_model: 
    total_lambda_tau= Uniform('total_lambda_tau', lower=0, upper=100000)
    total_lambda_mu = Uniform('total_lambda_mu', lower=0, upper=1000000)
    total_lambda = Normal('total_lambda', mu=total_lambda_mu, tau=total_lambda_tau)
    total = Poisson('total', mu=total_lambda, observed=totalCounts) 

    loss_lambda_factor = Uniform('loss_lambda_factor', lower=0, upper=1)
    error = Poisson('error', mu=total_lambda*loss_lambda_factor, observed=successCounts) 

with success_model: 
    step =  Metropolis() 
    success_samples = sample(20000, step) #, start)
    
plt.figure(figsize=(10, 10)) 
_ = traceplot(success_samples)

plt.figure(figsize=(10, 10)) 
_ = traceplot(success_samples[10000:])

with Model() as success_model: 
    total_lambda_tau= Uniform('total_lambda_tau', lower=0, upper=100)
    total_lambda_mu = Uniform('total_lambda_mu', lower=0, upper=1000)
    total_lambda = Normal('total_lambda', mu=total_lambda_mu, tau=total_lambda_tau)
    total = Poisson('total', mu=total_lambda, observed=totalCounts) 

    loss_lambda_factor = Uniform('loss_lambda_factor', lower=0, upper=1)
    error = Poisson('error', mu=total_lambda*loss_lambda_factor, observed=successCounts) 

with success_model: 
    step =  Metropolis() 
    success_samples = sample(20000, step) #, start)
    
plt.figure(figsize=(10, 10)) 
_ = traceplot(success_samples[10000:])

with success_model: 
    start = pm.find_MAP()
    step =  pm.CompoundStep([pm.Metropolis([total_lambda_mu]),
                             pm.Metropolis([total_lambda_tau]),
                             pm.Metropolis([total_lambda]),
                             pm.Metropolis([loss_lambda_factor]),
                             ]) 
    success_samples = sample(20000, step, start)

plt.figure(figsize=(10, 10)) 
_ = traceplot(success_samples[10000:])

with success_model: 
    step =  pm.NUTS() 
    success_samples = sample(20000, step)
plt.figure(figsize=(10, 10)) 
_ = traceplot(success_samples[10000:])

import theano.tensor as T
with Model() as success_model: 
    loss_lambda_rate = pm.Flat('loss_lambda_rate')
    error = Poisson('error', mu=totalCounts*T.exp(loss_lambda_rate), observed=successCounts)

with success_model: 
    step =  Metropolis() 
    success_samples = sample(20000, step) #, start)
    
plt.figure(figsize=(10, 10)) 
_ = traceplot(success_samples[10000:])

import numpy as np
plt.hist(np.exp(success_samples['loss_lambda_rate'][10000:]));