Distribution proposal/enhancement: GaussianRandomWalk for irregularly spaced observations

[drbenvincent]

As far as I understand, the existing GaussianRandomWalk distribution implicitly assumes regularly spaced observations in time because the variance at each timestep is constant.

The proposal would be to either create a new distribution for irregularly spaced observations or to generalise the implementation of the current distribution.

At this point I am not familiar enough with the backend to implement this, but the basic approach would be to:

• provide an additional vector of observation time points
• use the fairly basic maths of the sum of normally distributed variables to adjust for the temporal gap between observations.

So it would still be a cumulative sum of normals, but as far as I can work out, the mu and std would be multiplied by the time delta between observations.

Quick bit of code to illustrate...

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import arviz as az
import aesara.tensor as aet

# create regularly spaced observations from a Gaussian random walk
N = 50
true_mu, true_std = 0.1, 1
t = np.arange(N)
y = (np.random.randn(N) * true_std + true_mu).cumsum()

# throw away some random observations. Those remaining will be irregularly spaced
keep = np.random.choice(np.arange(N), size=20, replace=False)
keep = np.sort(keep)
t = t[keep]
y = y[keep]
plt.plot(t, y, 'o');

def GaussianRandomWalk(name, mu, std, t):
    multiplier = aet.concatenate([[1.], aet.extra_ops.diff(t)])
    Δobs = pm.Normal("Δobs", mu*multiplier, std*multiplier, shape=len(t)) 
    return pm.Deterministic(name, aet.extra_ops.cumsum(Δobs))

with pm.Model() as model:
    mu = pm.Normal('mu', 0, 1)
    std = pm.HalfNormal('std', 1)      
    μ_obs = GaussianRandomWalk('gwr', mu, std, t)
    pm.Normal('y', μ_obs, 0.1, observed=y)
    idata = pm.sample(tune=2000)
    idata.extend(pm.sample_posterior_predictive(idata, var_names=["y"]))

# visual diagnostics
az.plot_posterior(idata, var_names=["mu", "std"], ref_val=[true_mu, true_std]);

# visualise posterior predictive + data
fig, ax = plt.subplots()
az.plot_hdi(t, ppc.posterior_predictive.y, hdi_prob=0.95, fill_kwargs={'alpha':0.25}, smooth=False, ax=ax)
ax.plot(t, y, 'ko');

Which seems to work pretty well.

[twiecki]

CC @canyon289 @mjhajharia

[drbenvincent]

Bonus points to demonstrate/implement how to do posterior prediction with interpolation between observations. I think it's less trivial than a straight regression model?

[canyon289]

This is a nice suggestion! If i'm slow on feedback its because im currently (more like struggling) to implement a 1d fixed step implementation. Let me push that over the line (with lots of help from Ricardo and Michael, many thanks as welll) and I'll be back here.

[drbenvincent]

Bump. Wondered if there were any more thoughts on this idea. In particular, @ricardoV94?

PS. I think my example would be better if the level of observation noise was much higher, and the sigma on the random walk was fixed and low. You'd hope to see some effects of smoothing that way.

[ricardoV94]

This would require a new distribution as the GRW uses a baseclass that is not gaussian specific, and doesn't allow time varying sigmas.

Alternatively we could use imputation for the "missing timepoints" (whenever we bring back imputation for multivariate distributions), but that could be potentially inefficient if you need many "fake" timepoints. The advantage of this is that it works for non gaussian random walks.

If you are so inclined, I suggest pushing something over to pymc-experimental. As you mentioned, you need to write the function that creates a cumulative sum of Normals with the right sigmas given the time delta.

The logp can be obtained automatically from the generative process:

pm.logp(pt.cumsum(pm.Normal.dist(0, sigma=[1, 2, 3, 2])), value=[0, 0, 0, 0, 0]).eval()

Look at Censored for a simple case where we create a fully fledged PyMC distribution from a simple generative process.

[ricardoV94]

If you don't want a fancy Distribution you can use CustomDist:

import pymc as pm

def dist(mu, sigma, dt, size):
  sigma_adj = ...
  raw_innov = pm.Normal.dist(mu, sigma_adj, shape=(size+len(dt)))
  grw = raw_innov.cumsum(-1)
  return grw

with pm.Model() as m:
  ...
  pm.CustomDist("grw", mu, sigma, dt, dist=dist, observed=...)

[ricardoV94]

Here is the original example with a CustomDist:

https://gist.github.com/ricardoV94/4dc16748f99edcc162489352eed8fac7