EulerMaruyama testval

Questions
1

I am trying to estimate the parameters of a SDE using the EulerMaruyama class in pymc3. The SDE is defined as:

dp_t = s p_t (1-p_t) + \sqrt{p_t(1-p_t)/N}dW

where p\in[0,1].
When I encode this in Pymc3 using;

pt = EulerMaruyama("pt", dt, wright_fisher_sde, (s,), shape=(T,), testval=0.1*tt.ones(T))

the model fails to initialize sampling, complaining about SamplingError: Initial evaluation of model at starting point failed!. Inspecting the initial values of the variables, it seems that the initial value of pt (an array with shape (T,) ) includes negative values as well as values above 1, which makes the variance part of the SDE negative, which of course breaks everything.
However, I am giving it a testval value that should be valid. However, the sampler does not take it into account.

How can I solve this?

thanks for any help

2

Would you be able to share a self-contained example that reproduces this? I’d be happy to take a look!

3

Sure, thanks a lot. Here’s a small example.

import arviz as az
import pymc3 as pm
import numpy as np
import theano.tensor as tt
from pymc3.distributions.timeseries import EulerMaruyama

# generate some data
s = .1
N = 100
dt = 1e-1
T = 50


# time series
p = 0.1
p_t = []

# simulate
for i in range(T):
    p += dt * s * p * (1-p) + np.sqrt(dt*p*(1-p)/N) * np.random.normal()
    p_t.append(p)

p_t = np.array(p_t)

# z_t noisy observation
z_t = p_t + np.random.normal(size=p_t.size) * 5e-3

Now, we implement the model:

def wright_fisher_sde(p,s):
    return  s*p*(1-p),np.sqrt((p*(1-p)/N))

with pm.Model() as model:
    s = pm.Normal("s",0,1)
    pt = EulerMaruyama("pt", dt, wright_fisher_sde, (s,), shape=(T,), testval=0.1*tt.ones(T))

    zh = pm.Normal("zh", mu=pt, sigma=5e-3, observed=z_t)

    trace = pm.sample(1000,return_inferencedata=True,cores=1,chains=2)

This fails with error: SamplingError: Initial evaluation of model at starting point failed!.

Thanks again for any help you can provide

4

I was able to reproduce your problem and would have expected it to work fine the way you wrote it. I’m not sure why testval was getting ignored. Either way, you can use the start keyword argument for pm.sample to properly initialize the sampler. Try this:

with pm.Model() as model:
    s  = pm.Normal("s",0,1)
    pt = EulerMaruyama("pt", dt, wright_fisher_sde, (s,), shape=(T,), testval=0.1*np.ones(T))

    zh = pm.Normal("zh", mu=pt, sigma=5e-3, observed=z_t)
    
    start = {
    's' :0.5,
    'pt': 0.1*np.ones(T)
    }

    trace = pm.sample(1000, return_inferencedata=True,cores=1,chains=2, start=start)
1 Like
5

Indeed it does solve the problem! I forgot about the start keyword.

Thanks a lot!

6

Sampling will use the testval (soon to be deprecated in favor of initval, FYI), but during default initialization (i.e., init="jitter+adapt_diag") the initial value is jittered. This can cause problems when the initial value is near invalid or otherwise problematic values. To avoid the jitter, use another initialization. I often use "adapt_diag" in these cases.

trace = pm.sample(1000,
                  init="adapt_diag",
                  return_inferencedata=True,
                  cores=1,
                  chains=2)
2 Likes
7

that explains it!

1 Like