Learning Time-Dependent Models

Time-Dependent Parameters

In addition to learning static parameters, QInfer can be used to learn the values of parameters that change stochastically as a function of time. In this case, the model parameter vector \(\vec{x}\) is interpreted as a time-dependent vector \(\vec{x}(t)\) representing the state of an underlying process. The resulting statistical problem is often referred to as state-space estimation. By using an appropriate resampling algorithm, such as the Liu-West algorithm [LW01], state-space and static parameter estimation can be combined such that a subset of the components of \(\vec{x}\) are allowed to vary with time.

QInfer represents state-space filtering by the use of the Simulatable.update_timestep() method, which samples how a model parameter vector is updated as a function of time. In particular, this method is used by SMCUpdater to draw samples from the distribution \(f\)

\[\vec{x}(t_{\ell+1}) \sim f(\vec{x}(t_{\ell}, e(t_{\ell})),\]

where \(t_{\ell}\) is the time at which the experiment \(e_{\ell}\) is measured, and where \(t_{\ell+1}\) is the time step immediately following \(t_{\ell}\). As this distribution is in general dependent on the experiment being performed, update_timestep() is vectorized in a manner similar to likelihood() (see Designing and Using Models for details). That is, given a tensor \(X_{i,j}\) of model parameter vectors and a vector \(e_k\) of experiments, update_timestep() returns a tensor \(X_{i,j,k}'\) of sampled model parameters at the next time step.

Random Walk Models

As an example, RandomWalkModel implements update_timestep() by taking as an input a Distribution specifying steps \(\Delta \vec{x} = \vec{x}(t + \delta t) - \vec{x}(t)\). An instance of RandomWalkModel decorates another model in a similar fashion to other derived models. For instance, the following code declares a precession model in which the unknown frequency \(\omega\) changes by a normal random variable with mean 0 and standard deviation 0.005 after each measurement.

>>> from qinfer import (
...     SimplePrecessionModel, RandomWalkModel, NormalDistribution
... )
>>> model = RandomWalkModel(
...     underlying_model=SimplePrecessionModel(),
...     step_distribution=NormalDistribution(0, 0.005 ** 2)
... )

We can then draw simulated trajectories for the true and estimated value of \(\omega\) using a minor modification to the updater loop discussed in Sequential Monte Carlo.

model = RandomWalkModel(
    # Note that we set a minimum frequency of negative
    # infinity to prevent problems if the random walk
    # causes omega to cross zero.
    underlying_model=SimplePrecessionModel(min_freq=-np.inf),
    step_distribution=NormalDistribution(0, 0.005 ** 2)
)
prior = UniformDistribution([0, 1])
updater = SMCUpdater(model, 2000, prior)

expparams = np.empty((1, ), dtype=model.expparams_dtype)

true_trajectory = []
est_trajectory = []

true_params = prior.sample()

for idx_exp in range(400):
    # We don't want to grow the evolution time to be arbitrarily
    # large, so we'll instead choose a random time up to some
    # maximum.
    expparams[0] = np.random.random() * 10 * np.pi
    datum = model.simulate_experiment(true_params, expparams)
    updater.update(datum, expparams)

    # We index by [:, :, 0] to pull off the index corresponding
    # to experiment parameters.
    true_params = model.update_timestep(true_params, expparams)[:, :, 0]

    true_trajectory.append(true_params[0])
    est_trajectory.append(updater.est_mean())

plt.plot(true_trajectory, label='True')
plt.plot(est_trajectory, label='Est.')
plt.legend()
plt.xlabel('# of Measurements')
plt.ylabel(r'$\omega$')

plt.show()

(Source code, svg, pdf, hires.png, png)

../_images/timedep-1.svg

Specifying Custom Time-Step Updates

The RandomWalkModel example above is somewhat unrealistic, however, in that the step distribution is independent of the evolution time. For a more reasonable noise process, we would expect that \(\mathbb{V}[\omega(t + \Delta t) - \omega(t)] \propto \Delta t\). We can subclass SimplePrecessionModel to add this behavior with a custom update_timestep() implementation.

class DiffusivePrecessionModel(SimplePrecessionModel):
    diffusion_rate = 0.0005 # We'll multiply this by
                             # sqrt(time) below.

    def update_timestep(self, modelparams, expparams):
        step_std_dev = self.diffusion_rate * np.sqrt(expparams)
        steps = step_std_dev * np.random.randn(
            # We want the shape of the steps in omega
            # to match the input model parameter and experiment
            # parameter shapes.
            # The axis of length 1 represents that this model
            # has only one model parameter (omega).
            modelparams.shape[0], 1, expparams.shape[0]
        )
        # Finally, we add a new axis to the input model parameters
        # to match the experiment parameters.
        return modelparams[:, :, np.newaxis] + steps

model = DiffusivePrecessionModel()
prior = UniformDistribution([0, 1])
updater = SMCUpdater(model, 2000, prior)

expparams = np.empty((1, ), dtype=model.expparams_dtype)

true_trajectory = []
est_trajectory = []

true_params = prior.sample()

for idx_exp in range(400):
    expparams[0] = np.random.random() * 10 * np.pi
    datum = model.simulate_experiment(true_params, expparams)
    updater.update(datum, expparams)

    true_params = model.update_timestep(true_params, expparams)[:, :, 0]

    true_trajectory.append(true_params[0])
    est_trajectory.append(updater.est_mean())

plt.plot(true_trajectory, label='True')
plt.plot(est_trajectory, label='Est.')
plt.legend()
plt.xlabel('# of Measurements')
plt.ylabel(r'$\omega$')

plt.show()

(Source code, svg, pdf, hires.png, png)

../_images/timedep-2.svg

Learning Walk Parameters

In the above examples, the diffusion distribution was treated as exactly known by the model. We can also parameterize this distribution, adding its parameters to model to be learned as well. GaussianRandomWalkModel is a built in model similar to RandomWalkModel. It is more restrictive in the sense that it is limited to gaussian time-step updates, but more general in that it has the ability to automatically append a parameterization of the gaussian time-step distribution, either diagonal or dense, to the underlying model.

For example suppose that we have a coin whose bias is taking a random walk in time with an unknown diffusion constant. To avoid exiting the allowable space of biases, \([0,1]\), we transform to inverse-logit space before taking a gaussian step, and transform back to the probability interval after each step.

import numpy as np
from scipy.special import expit, logit
from qinfer import (
   CoinModel, BinomialModel, GaussianRandomWalkModel,
   UniformDistribution, SMCUpdater
)

# Put a random walk on top of a binomial coin model
model = GaussianRandomWalkModel(
   BinomialModel(CoinModel()),
   model_transformation=(logit, expit)
)

# Generate some data with a true diffusion 0.05
true_sigma_p = 0.05
Nbin = 10
p = expit(logit(0.5) + np.cumsum(true_sigma_p * np.random.normal(size=300)))
data = np.random.binomial(Nbin, 1-p)

# Analyse the data
prior = UniformDistribution([[0.2,0.8],[0,0.1]])
u = SMCUpdater(model, 10000, prior)
ests, stds = np.empty((data.size+1, 2)), np.empty((data.size+1, 2))
ts = np.arange(ests.shape[0])
ests[0,:] = u.est_mean()
for idx in range(data.size):
    expparam = np.array([Nbin]).astype(model.expparams_dtype)
    u.update(np.array([data[idx]]), expparam)
    ests[idx+1,:] = u.est_mean()
    stds[idx+1,:] = np.sqrt(np.diag(u.est_covariance_mtx()))

plt.figure(figsize=(10,10))
plt.subplot(2,1,1)
u.plot_posterior_marginal(1)
plt.title('Diffusion parameter posterior')

plt.subplot(2,1,2)
plt.plot(ts, ests[:,0], label='estimated')
plt.fill_between(ts, ests[:,0]-stds[:,0], ests[:,0]+stds[:,0],
alpha=0.2, antialiased=True)
plt.plot(ts[1:], p, '--', label='actual')
plt.legend()
plt.title('Coin bias vs. time')
plt.show()

(Source code, svg, pdf, hires.png, png)

../_images/timedep-3.svg

As a second example, consider a 5-sided die for which the 3rd, 4th and 5th sides are taking a correlated gaussian random walk, and the other two sides are constant. We can attempt to learn the six parameters of the cholesky factorization of the random walk covariance matrix as we track the drift of the die probabilities.

import numpy as np
from qinfer.utils import to_simplex, from_simplex, sample_multinomial
from qinfer import (
   NDieModel, MultinomialModel, GaussianRandomWalkModel,
   UniformDistribution, ConstrainedSumDistribution, SMCUpdater, ProductDistribution
)

# Put a random walk on top of a multinomial die model
randomwalk_idxs = [2,3,4] # only these sides of the die are taking a walk
model = GaussianRandomWalkModel(
    MultinomialModel(NDieModel(5)),
    model_transformation=(from_simplex, to_simplex),
    diagonal=False,
    random_walk_idxs = randomwalk_idxs
)

# Generate some data with some true covariance matrix
true_cov = 0.1 * np.random.random(size=(3,3))
true_cov = np.dot(true_cov, true_cov.T)
Nmult = 40
ps = from_simplex(np.array([[0.1,0.2,0.2,0.4,.1]] * 200))
ps[:, randomwalk_idxs] += np.random.multivariate_normal(np.zeros(3), true_cov, size=200).cumsum(axis=0)
ps = to_simplex(ps)
expparam = np.array([(0,Nmult)],dtype=model.expparams_dtype)
data = sample_multinomial(Nmult, ps.T).T

# Analyse the data
prior = ProductDistribution(
    ConstrainedSumDistribution(UniformDistribution([[0,1]] * 5)),
    UniformDistribution([[0,0.2]] * 6)
)
u = SMCUpdater(model, 10000, prior)
ests, stds = np.empty((data.shape[0]+1, model.n_modelparams)), np.empty((data.shape[0]+1, model.n_modelparams))
ts = np.arange(ests.shape[0])
ests[0,:] = u.est_mean()
for idx in range(data.shape[0]):
    expparam = np.array([(0,Nmult)],dtype=model.expparams_dtype)
    outcome = np.array([(data[idx],)], dtype=model.domain(expparam)[0].dtype)
    u.update(outcome, expparam)
    ests[idx+1,:] = u.est_mean()
    stds[idx+1,:] = np.sqrt(np.diag(u.est_covariance_mtx()))

true_chol = np.linalg.cholesky(true_cov)
k = 1
plt.figure(figsize=(10,10))
for idx, coord in enumerate(zip(*model._srw_tri_idxs)):
    i, j = coord
    plt.subplot(3,3,i*3 + j + 1)
    u.plot_posterior_marginal(5 + idx)
    plt.axvline(true_chol[i,j],color='b')
plt.show()

plt.figure(figsize=(12,10))
color=iter(plt.cm.Vega10(range(5)))
for idx in range(5):
    c=next(color)
    plt.plot(ps[:, idx], '--', label='$p_{}$ actual'.format(idx), color=c)
    plt.plot(ests[:,idx], label='$p_{}$ estimated'.format(idx), color=c)
    plt.fill_between(range(len(ests)), ests[:,idx]-stds[:,idx], ests[:,idx]+stds[:,idx],alpha=0.2, color=c, antialiased=True)
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.show()

(Source code)