Skip to content

AR1 log likelihood #2284

@eph

Description

@eph

I think the log likelihood calculation for the AR1 model is wrong.

class AR1(distribution.Continuous):

    ...

    def logp(self, x):
        k = self.k
        tau_e = self.tau_e

        x_im1 = x[:-1]
        x_i = x[1:]
        boundary = Normal.dist(0., tau=tau_e).logp

        innov_like = Normal.dist(k * x_im1, tau=tau_e).logp(x_i)
        return boundary(x[0]) + tt.sum(innov_like) + boundary(x[-1])

The last term in the log likelihood should not be there.

import pymc3
from scipy.stats import norm

with pymc3.Model() as t:
    my_ar1 = pymc3.AR1('y',0,1,shape=2)

print('pymc3 log lik     :', t['y'].logp({'y':np.zeros(2,)}) )
print('scipy log lik(n=2):', norm.logpdf([0,0]).sum() )
print('scipy log lik(n=3):', norm.logpdf([0,0,0]).sum() )

yields

pymc3 log lik     : -2.756815599614018
scipy log lik(n=2): -1.83787706641
scipy log lik(n=3): -2.75681559961

I will try to write some tests following what's in test_distributions.py. BTW, this is fantastic software -- thanks for all your efforts.

Metadata

Metadata

Assignees

No one assigned

    Labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions