%matplotlib inline
from pymc3 import Normal, sample, Model, traceplot, plots, Metropolis, Potential, variational, HalfCauchy,summary
from theano import scan, shared
from scipy import optimize
import theano.tensor as tt
import matplotlib.pyplot as plt
import numpy as np
y = shared(np.array([15, 10, 16, 11, 9, 11, 10, 18], dtype=np.float32))
with Model() as arma_model:
sigma = HalfCauchy('sigma', 5)
theta = Normal('theta', 0, sd=2)
phi = Normal('phi', 0, sd=2)
mu = Normal('mu', 0, sd=10)
err0 = y[0] - (mu + phi*mu)
def calc_next(last_y, this_y, err, mu, phi, theta):
nu_t = mu + phi*last_y + theta*err
return this_y - nu_t
err, _ = scan(fn=calc_next,
sequences=dict(input=y, taps=[-1,0]),
outputs_info=[err0],
non_sequences=[mu, phi, theta])
like = Potential('like', Normal.dist(0, sd=sigma).logp(err))
Applied log-transform to sigma and added transformed sigma_log to model.
with arma_model:
trace = sample(2000)
Assigned NUTS to sigma_log Assigned NUTS to theta Assigned NUTS to phi Assigned NUTS to mu [-----------------100%-----------------] 2000 of 2000 complete in 35.4 sec
traceplot(trace[1000:]);