import pymc3 as pm
import arviz as az
import numpy as np
import matplotlib.pyplot as plt
from pymc3.step_methods.arraystep import BlockedStep
class UpdateIndependentNoiseTermStep(BlockedStep):
def __init__(self, epsilon, sigma, model=None):
pm.modelcontext(model)
self.vars = [epsilon]
self.epsilon = epsilon
self.sigma = sigma
pass
def step(self, point: dict):
sigma = self.sigma
if not isinstance(1.0*self.sigma, float):
sigma = np.exp(point[self.sigma.transformed.name])
pass
point[self.epsilon.name] = pm.Normal.dist(mu=0, sigma=sigma).random()
return point
pass
import pymc3 as pm
import arviz as az
import numpy as np
import matplotlib.pyplot as plt
from pymc3.step_methods.arraystep import BlockedStep
from custom import UpdateIndependentNoiseTermStep
data = 100 * np.exp(1) ** pm.Normal.dist(mu=0, sigma=0.2).random(size=2000)
plt.plot(data)

with pm.Model() as model:
alpha = pm.Normal(r"$\alpha$", mu=100, sigma=10)
sigma = pm.HalfStudentT(r"$\sigma$", nu=3, sigma=0.15)
epsilon = pm.Normal(r"$\varepsilon$", mu=0.0, sigma=sigma)
update_fake_step = UpdateIndependentNoiseTermStep(epsilon=epsilon, sigma=sigma)
alpha_eps = pm.Deterministic(r"$\alpha_{\varepsilon}$", alpha * np.exp(1) ** epsilon)
obs = pm.Normal("obs", mu=alpha_eps, sigma=1, observed=data)
step = pm.Metropolis(vars=[alpha, sigma])
trace = pm.sample(4000, cores=4, chains=4, tune=10000, step=[step, update_fake_step], init='advi', return_inferencedata=True)
pass
az.plot_trace(trace)

az.summary(trace)
| mean | sd | hdi_3% | hdi_97% | mcse_mean | mcse_sd | ess_bulk | ess_tail | r_hat | |
|---|---|---|---|---|---|---|---|---|---|
| πΌ | 102.131 | 5.782 | 90.983 | 113.800 | 0.244 | 0.174 | 607.0 | 756.0 | 1.01 |
| π | 0.002 | 0.273 | -0.422 | 0.449 | 0.002 | 0.004 | 15976.0 | 4481.0 | 1.00 |
| π | 0.171 | 0.218 | 0.000 | 0.445 | 0.006 | 0.004 | 1251.0 | 1517.0 | 1.00 |
| πΌπ | 108.864 | 142.697 | 61.872 | 152.625 | 1.751 | 1.238 | 9790.0 | 4724.0 | 1.00 |
Roughly like this.