Cacher is the code snippet organizer for pro developers

import numpy as np
import pymc3 as mc
import scipy as sci
import theano.tensor as th

np.random.seed(13)

n = 10
tau_scale = 2
tau0 = sci.stats.expon.rvs() * tau_scale
mu0 = np.random.randn(n) / np.sqrt(tau0)
x0 = mu0 + np.random.randn(n)

with mc.Model() as model1:
    tau = mc.Exponential('tau', lam=tau_scale)
    mu = mc.Normal('mu', tau=tau, shape=(n,))
    mc.Normal('x', mu=mu, observed=x0)

with mc.Model() as model2:
    tau = mc.Exponential('tau', lam=tau_scale)
    mu_z = mc.Normal('mu_z', shape=(n,))
    mu = mc.Deterministic('mu', mu_z / th.sqrt(tau))
    mc.Normal('x', mu=mu, observed=x0)
    
stan_model = '''
data {
  int<lower=0> n; 
  real<lower=0> tau_scale; 
  real x0[n];
} 
parameters {
  real<lower=0> tau; 
  vector[n] mu;
} 
model {
  tau ~ exponential(tau_scale);
  mu ~ normal(0, 1/tau);
  for (i in 1:n) x0[i] ~ normal(mu[i], 1); 
}
'''

import pystan

stan_data = {'n': n, 'tau_scale': tau_scale, 'x0': x0}

fit = pystan.stan(model_code=stan_model, data=stan_data,
                  iter=1000, chains=4)
print(fit)

with model1:
    advi_fit = mc.variational.advi(n=100000)

with model2:
    advi_fit2 = mc.variational.advi(n=100000)    

def infer(model):
    with model:
        map_ = mc.find_MAP(fmin=sci.optimize.fmin_l_bfgs_b)
        step = mc.NUTS(scaling=map_)
        trace = mc.sample(100, step=step, start=map_, progressbar=False)
        step = mc.NUTS(scaling=trace[-1])
        return mc.sample(11000, step=step, start=trace[-1], progressbar=False)

trace1 = infer(model1)
trace2 = infer(model2)

with model2:
    trace3 = mc.sample(100000, step=mc.Metropolis(), progressbar=False,
                       start=mc.find_MAP(fmin=sci.optimize.fmin_l_bfgs_b))

samples_tau1 = trace1['tau'][1000:]
samples_tau2 = trace2['tau'][1000:]
samples_tau3 = trace3['tau'][10000:]

print()
print( 'pymc3 version: ' + mc.__version__)
print()
print( 'Model 1 NUTS tau')
print( 'Mean: {0:3.1f}'.format(samples_tau1.mean()))
print( 'Standard Deviation: {0:3.1f}'.format(samples_tau1.std()))
print( 'Median {0:3.1f}'.format(np.percentile(samples_tau1, 50)))
print()
print( 'Model 2 NUTS tau')
print( 'Mean: {0:3.1f}'.format(samples_tau2.mean()))
print( 'Standard Deviation: {0:3.1f}'.format(samples_tau2.std()))
print( 'Median {0:3.1f}'.format(np.percentile(samples_tau2, 50)))
print()
print( 'Model 2 Metropolis tau')
print( 'Mean: {0:3.1f}'.format(samples_tau3.mean()))
print( 'Standard Deviation: {0:3.1f}'.format(samples_tau3.std()))
print( 'Median {0:3.1f}'.format(np.percentile(samples_tau3, 50)))
print()
print( 'Model 1 ADVI tau')
print( 'Mean: {0:3.1f}'.format(np.exp(advi_fit[0]['tau_log'])))
print( 'Standard Deviation: {0:3.1f}'.format(np.exp(advi_fit[1]['tau_log'])))
print()
print( 'Model 2 ADVI tau')
print( 'Mean: {0:3.1f}'.format(np.exp(advi_fit2[0]['tau_log'])))
print( 'Standard Deviation: {0:3.1f}'.format(np.exp(advi_fit2[1]['tau_log'])))