Autocorrelation in a model

Questions
1

Is there a practical way to implement autocorrelation in a pymc3 model?

What I want to implement is:

x_0~normal(0,s^2)
x_i~normal( ro*x_i-1 , (1-ro^2)*s^2 ) for i=1,2,3,…n-1

In my current application, n=106

I can’t figure out how to implement this type of model in pymc3. I am hoping there is a simple way.

Thanks,
Wayne Hajas

2

There is AR(n) distribution for this purpose, see eg http://docs.pymc.io/notebooks/AR.html

3

Thanks Junpeng. I took a look at the AR function.

Unfortunately, the autocorrelation is just part of my model. I need to estimate all the x-values as well as the ro and sigma-values. As far as I can tell, the AR function will not help me estimate the x-values.

This model has been implemented in pymc2. I suspect the implementation is also possible in pymc3 - but I am just not seeing it.

Any assistance is appreciated.
Wayne

4

could you post the pymc2 code?

5

I am in the process of cleaning up the pymc2 code. I haven’t fully tested this version but this is what I have right now. I am generating auto-correlated values of LogRecruit. Thanks for taking a look at this.

def Recruitment(self):
    '''Log of Relative annual Recruitment'''

    #Will be subjected to antiLogit transformation
    self.T = Normal('T', 7.78, 0.667)
    #Coefficient for autocorrelation
    self.Ro = Lambda('Ro', lambda t=self.T: expit(t))
    self.Omega = Lambda('Omega', lambda t=self.Ro: np.sqrt(1 - t * t))

    #Standard deviation of recruitment rate
    self.lnSigmaRecruit = Normal('lnSigmaRecruit', 2.5, 1)
    self.tauRecruit = Lambda(
        'tauRecruit', lambda t=self.lnSigmaRecruit: np.exp(-2 * t))

    self.LogRecruit = numpy.empty(self.nyear, dtype=object)
    self.LogRecruit[0] = Normal('LogRecruit_' + str(self.MinYear), 0,
                                self.tauRecruit)
    for y in range(1, self.nyear):
        self.LogRecruit[y] = Normal('LogRecruit_' + str(self.MinYear + y),
                                    self.Ro * self.LogRecruit[y - 1],
                                    self.tauRecruit / self.Omega / self.Omega)
6

What you are descripting is exactly a AR model. For example, the PyMC3 code similar to about would be something like:

with pm.Model():
    T = pm.Normal('T', 7.78, tau=0.667)
    Ro = expit(T)
    Omega = tt.sqrt(1.-Ro**2)
    lnSigmaRecruit = pm.Normal('lnSigmaRecruit', 2.5, tau=1)
    tauRecruit = tt.exp(-2. * lnSigmaRecruit)
    data = pm.AR('y', 
                 Ro, 
                 tau=tauRecruit / Omega / Omega, 
                 init=pm.Normal.dist(0., tau=tauRecruit), 
                 shape=nyear)
1 Like
7

Maybe I am misunderstanding documentation and examples like https://github.com/pymc-devs/pymc3/blob/master/docs/source/notebooks/AR.ipynb

As I understand, the code you generated will let me estimate T and lnSigmaRecruit. From those estimates, I can also estimate Ro, Omega and tauRecruit.

I also need to generate LogRecruit (from my pymc2 code). LogRecruit represents one process. I have another array of nodes representing another process. Then I multiply the two arrays together to represent the combined processes - and then this product is compared against data.

Does pm.AR generate an array of LogRecruit-nodes?

I hope that bit of explanation helps. Let me know if I am misunderstanding pm.AR

Thanks for looking at this,
Wayne

8

Just to be clear:

data = pm.AR('y', 
                 Ro, 
                 tau=tauRecruit / Omega / Omega, 
                 init=pm.Normal.dist(0., tau=tauRecruit), 
                 shape=nyear)

is the LogRecruit as in the pymc2 code.

So yes, you can generate two random arrays similarly with pm.AR and multiply them.

The example above is conditioned on an observed, but pm.AR works similar when you are treating it as a latent Variable following some distribution. See the docstring: http://docs.pymc.io/api/distributions/timeseries.html#pymc3.distributions.timeseries.AR

9

Thanks very much Junpeng!
pm.AR works as you said it would.
Now all my problems are my own fault.

Cheers,
Wayne