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.
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.
'''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 = Normal('LogRecruit_' + str(self.MinYear), 0,
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)
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