In [13]:
%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
In [2]:
y = shared(np.array([15, 10, 16, 11, 9, 11, 10, 18], dtype=np.float32))
In [3]:
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.
In [15]:
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
In [16]:
traceplot(trace[1000:]);