pints icon indicating copy to clipboard operation
pints copied to clipboard

NUTS / transform / Goodwin failing

Open ben18785 opened this issue 3 years ago • 0 comments

I often seem to be getting an error when trying to run NUTS on the Goodwin-Oscillator example. The below code reproduces this. The failure only occurs when I use transformations but this doesn't seem to affect other (non-gradient-based) samplers.

import pints
import pints.toy
import numpy as np

model = pints.toy.GoodwinOscillatorModel()

# Create data
real_parameters = model.suggested_parameters()
times = np.linspace(0, 200, 1000)
values = model.simulate(real_parameters, times)
noise1 = 0.0001
noise2 = 0.003
noise3 = 0.05
noisy_values = np.array(values, copy=True)
noisy_values[:, 0] += np.random.normal(0, noise1, len(times))
noisy_values[:, 1] += np.random.normal(0, noise2, len(times))
noisy_values[:, 2] += np.random.normal(0, noise3, len(times))

# Create an object with links to the model and time series
problem = pints.MultiOutputProblem(model, times, values)

# Create a log posterior
log_prior = pints.MultivariateGaussianLogPrior(real_parameters,
                       np.diag([0.1, 0.1, 0.1, 0.1, 0.1]))
log_likelihood = pints.GaussianKnownSigmaLogLikelihood(problem,
                       [noise1, noise2, noise3])
log_posterior = pints.LogPosterior(log_likelihood, log_prior)

# Starting point that often gives errors
xs = [np.array([1.99890288, 4.01490882, 0.07433401, 0.11583922, 0.09271098])]

# Fails only when I use a transformation
transform = pints.LogTransformation(n_parameters=len(xs[0]))
mcmc = pints.MCMCController(log_posterior, 1, xs,
                            method=pints.NoUTurnMCMC,
                            transformation=transform)
mcmc.run()

When I print the sampled values as sampling progresses, I see that they diverge leading to infs quite quickly. Here's an example set of the first two iterations before failure occurs:

[1.99890288 4.01490882 0.07433401 0.11583922 0.09271098] [2.95701368e-131 8.84698508e+059 inf inf inf]

The failure ultimately happens when log_prior.evaluateS1(x) is done. But, to me, this doesn't seem like the problem itself necessarily, because the sampler has already gone infinite at that point.

This seems like a really hard one to figure out as could be any combination of issues with NUTS, transformations, Goodwin oscillator or the priors!

ben18785 avatar Aug 26 '21 10:08 ben18785