Refine pm.timeseries

[junpenglao]

Currently, the timeseries docstrings are not indexed on the website under API Reference > Distributions > Timeseries. The example notebook is also quite scattered.

We should unify all the timeseries notebooks into a single notebook, and potentially refactor some code and tests:

• We only have test for EulerMaruyama, AR currently, more tests are needed.
  • The tests for AR() and AR1 should be further refined
• The general AR logp could be improved (without using theano.scan)
• Unifying multiple timeseries docs.
• Timeseries docstring (link to doc API session)
• Implement random method (see discourse)

[cshenton]

Hey, I've starting giving this a stab. The running implementation I've got works fine as a distribution, but I get strange theano errors when I try to estimate it's parameters as part of a model. I'm running into the following issues:

• find_MAP ValueError: Optimization error: max, logp or dlogp at max have non-finite values.
• sample ValueError: Bad initial energy: inf. The model might be misspecified.
• and variably at other times theano complaining about using subtensors where symbolic tensors were expected when just compiling the model, though I can't replicate those at this time.

In general I'm wondering how pymc deals with distributions over vectors? The distribution certainly generates data sensibly.

Here is the full implementation:

import numpy as np
import pymc3 as pm
import theano.tensor as tt

from pymc3.distributions import distribution

class AR1(distribution.Continuous):
    """An Autoregressive distribution (order 1) with gaussian noise.

    Args:
        decay (float): The period to period decay parameter.
        scale (float): The scale param of the underlying normal distribution.
    """
    def __init__(self, decay, scale, *args, **kwargs):
        super(AR1, self).__init__(*args, **kwargs)
        self.decay = tt.as_tensor_variable(decay)
        self.scale = tt.as_tensor_variable(decay)
        self.mode = tt.as_tensor_variable(0.)
        self.disc_add = np.frompyfunc(lambda x1, x2: x1*decay + x2, 2, 1)
        self.base = pm.Normal.dist(mu=0., sd=self.scale)

    def default(self):
        return np.zeros(self.shape)

    def logp(self, x):
        """Calculative the log probability of the sequence x."""
        x_init = x[0]
        x_now = x[1:]
        x_prev = x[:-1]

        # Log probability of the first observation
        logp_initial = self.base.logp(x_init)

        # Log probability of the subsequent sequence
        x_zeroed = x_now - self.decay * x_prev
        logp_sequence = self.base.logp(x_zeroed)

        logp_total = logp_initial + tt.sum(logp_sequence)
        return logp_total

    def random(self, size):
        """Generates draws of sequences from the distribution.

        Achieves this by drawing sequences of noise from the base distribution
        and taking the discounted cumulative sum of the noise over the length
        of the sequence.

        Args:
            size (int or list of ints): The sample size to draw. Note
                that the last dimension will be taken to be the sequence
                length.

        Returns:
            numpy.ndarray: An array of dimensions shape * length.
        """
        x_base = self.base.random(size=size)
        # Type casting necessary since we accumulate an autogenerated ufunc
        result = self.disc_add.accumulate(
            x_base, axis=-1, dtype=np.object).astype(np.float)
        return result

if __name__ == '__main__':
    dist = AR1.dist(decay=0.9, scale=1.0)
    sample_data = dist.random(1000)

    with pm.Model() as model:
        decay = pm.Normal('decay', mu=0.0, sd=0.3)
        scale = pm.HalfNormal('scale', sd=1.0)
        series = AR1('series', decay=decay, scale=scale, observed=sample_data)

[cshenton]

Is the issue potentially the use of a numpy ufunc in random? I know other distributions don't output a theano tensor from random, so I wasn't too worried about maintaining a pure theano computational graph (since it would be broken by casting to a numpy array anyway).

[cshenton]

Never mind, assigned self.scale = tt.as_tensor_variable(decay), doh.

[junpenglao]

Also, we need more test for timeseries.

[springcoil]

#2938 I opened a separate issue for this, just because I think it's well carved out work.

[twiecki]

I think this was fixed or at least improved recently, reopen if not.