import numpy as np
import pymc3 as pm
Simulate data from Gaussian with mean=2 and SD=5 and try to recover these parameters with PyMC3.
x = np.random.normal(2, 5, size=1000)
with pm.Model() as model:
mu = pm.Normal('mu', 0, sd=1e4)
sigma = pm.HalfCauchy('sigma', 25)
likelihood = pm.Normal('likelihood', mu, sd=sigma, observed=x)
Applied log-transform to sigma and added transformed sigma_log to model.
with model:
tr = pm.sample(1000)
Assigned NUTS to mu Assigned NUTS to sigma_log [-----------------100%-----------------] 1000 of 1000 complete in 0.7 sec
pm.summary(tr)
mu: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------- 2.174 0.203 0.008 [1.851, 2.509] Posterior quantiles: 2.5 25 50 75 97.5 |--------------|==============|==============|--------------| 1.836 2.069 2.184 2.301 2.497 sigma_log: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------- 1.609 0.042 0.001 [1.563, 1.651] Posterior quantiles: 2.5 25 50 75 97.5 |--------------|==============|==============|--------------| 1.565 1.592 1.609 1.623 1.655 sigma: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------- 5.003 0.280 0.009 [4.774, 5.212] Posterior quantiles: 2.5 25 50 75 97.5 |--------------|==============|==============|--------------| 4.782 4.912 4.996 5.069 5.231