AR(1) logp potentially using wrong variance

[ljmartin]

Hi all,
I posted a funny output from fitting AR(1) models in the discourse - and before I go ahead and start using the fix I thought Id check here for feedback.
The issue and example code are described here: https://discourse.pymc.io/t/diagnosing-ar-1-model-whose-results-are-way-off-course/4947

Currently the AR(1) model assumes the first observation comes from a normal distribution with precision equivalent to the noise (i.e. innovation) term. This is fine if the data was initiated at zero with no previous data. In most cases, an experimenter won't observe the initiation point and so their AR(1) data can start anywhere within the dispersion of that AR(1) process.

Perhaps a better specification then is to use the variance of an AR(1) process to calculate the likelihood of the first observation. To be honest I'm using this as reference for the variance of an AR(1) process: https://stats.stackexchange.com/questions/103405/prove-expression-for-variance-ar1

Basically it all boils down to changing this line: https://github.com/pymc-devs/pymc3/blob/60cca231e59bf882507021570fb925d1812eb095/pymc3/distributions/timeseries.py#L76

to something like this

var_ar1 = 1 / ((1-k**2)*tau_e)
sd_ar1 = tt.sqrt(var_ar1)
boundary = pm.Normal.dist(0., sigma=sd_ar1).logp

This solved the issue I described in the discourse and brings the highest posterior density of the model fit on that data back in line with other (non bayesian) techniques for estimating the SEM of an AR(1).
thanks for your time!