STAN, PyStan RuntimeError: Initialization failed

Question

I try to use STAN for inference of an SIR-model expressed as a three-dimensional SDE: The code for the model compiles without problems.

toy_model = """
functions{
    vector mu(vector x, real alpha, real beta) {
        vector[3] m;
        m[1] = -alpha*x[1]*x[2];
        m[2] = alpha*x[1]*x[2]-beta*x[2];
        m[3] = beta*x[2];
        return m;
    }

    matrix sigma(vector x, real alpha, real beta) {
        matrix[3, 3] sig;
        sig[1, 1] = sqrt(alpha*x[1]*x[2])/10;
        sig[1, 2] = 0;
        sig[1, 3] = 0;
        sig[2, 1] = -sqrt(alpha*x[1]*x[2])/10;
        sig[2, 2] = 0;
        sig[2 ,3] = sqrt(beta*x[2])/10;
        sig[3, 1] = 0;
        sig[3, 2] = 0;
        sig[3, 3] = -sqrt(beta*x[2])/10;
        return sig;
    }
}
data{
    int<lower=0> N; // number of timesteps
    real<lower=0> dt; // stepsize
    array[N] vector[3] x; // vector process
}
parameters {
    real<lower=0, upper=1> alpha;
    real<lower=0, upper=1> beta;
}
model {
for (n in 2:N) {
        x[n] ~ multi_normal(x[n-1]+dt*mu(x[n-1], alpha, beta), sqrt(dt)*sigma(x[n-1], alpha, beta));
    }
}
"""

But when I try to sample from the model with

post = stan.build(toy_model, data=input_data)

I get the following error message:

Sampling:   0%Sampling: Initialization failed.
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
/home/vinc777/masterthesis/two_variant_model/PyStan/try.ipynb Cell 6' in <module>
----> 1 fit = post.sample(num_chains=4, num_samples=1000)

File ~/.local/lib/python3.8/site-packages/stan/model.py:89, in Model.sample(self, num_chains, **kwargs)
     61 def sample(self, *, num_chains=4, **kwargs) -> stan.fit.Fit:
     62     """Draw samples from the model.
     63 
     64     Parameters in ``kwargs`` will be passed to the default sample function.
   (...)
     87 
     88     """
---> 89     return self.hmc_nuts_diag_e_adapt(num_chains=num_chains, **kwargs)

File ~/.local/lib/python3.8/site-packages/stan/model.py:108, in Model.hmc_nuts_diag_e_adapt(self, num_chains, **kwargs)
     92 """Draw samples from the model using ``stan::services::sample::hmc_nuts_diag_e_adapt``.
     93 
     94 Parameters in ``kwargs`` will be passed to the (Python wrapper of)
   (...)
    105 
    106 """
    107 function = "stan::services::sample::hmc_nuts_diag_e_adapt"
--> 108 return self._create_fit(function=function, num_chains=num_chains, **kwargs)

File ~/.local/lib/python3.8/site-packages/stan/model.py:311, in Model._create_fit(self, function, num_chains, **kwargs)
    308     return fit
    310 try:
--> 311     return asyncio.run(go())
    312 except KeyboardInterrupt:
    313     return

File ~/.local/lib/python3.8/site-packages/nest_asyncio.py:38, in _patch_asyncio.<locals>.run(main, debug)
     36 task = asyncio.ensure_future(main)
     37 try:
---> 38     return loop.run_until_complete(task)
     39 finally:
     40     if not task.done():

File ~/.local/lib/python3.8/site-packages/nest_asyncio.py:81, in _patch_loop.<locals>.run_until_complete(self, future)
     78 if not f.done():
     79     raise RuntimeError(
     80         'Event loop stopped before Future completed.')
---> 81 return f.result()

File /usr/lib/python3.8/asyncio/futures.py:178, in Future.result(self)
    176 self.__log_traceback = False
    177 if self._exception is not None:
--> 178     raise self._exception
    179 return self._result

File /usr/lib/python3.8/asyncio/tasks.py:280, in Task.__step(***failed resolving arguments***)
    276 try:
    277     if exc is None:
    278         # We use the `send` method directly, because coroutines
    279         # don't have `__iter__` and `__next__` methods.
--> 280         result = coro.send(None)
    281     else:
    282         result = coro.throw(exc)

File ~/.local/lib/python3.8/site-packages/stan/model.py:235, in Model._create_fit.<locals>.go()
    233         sampling_output.clear()
    234         sampling_output.write_line("<info>Sampling:</info> <error>Initialization failed.</error>")
--> 235         raise RuntimeError("Initialization failed.")
    236     raise RuntimeError(message)
    238 resp = await client.get(f"/{fit_name}")

RuntimeError: Initialization failed.

I have no idea where this error comes from and how to solve it. I already double-checked whether the input data matches the model, but I seems that it does. Hopefully, someone has an idea or already had the same problem.

Update: I managed to get a smaller example with the same problem I reduced the dimension and complexity of the SDE a bit. Using a normal distribution to sample from, with mean and standard deviation provided as vectors

second_model = """
functions {
    vector mu(vector y, real alpha, real beta, real gamma) {
        vector[2] m;
        m = y*(beta-alpha-gamma-beta);
        return m;
    }
    vector Sigma(vector y, real sigma) {
        vector[2] sig;
        sig[1] = sigma*(1-y[1]*y[2])/10;
        sig[2] = sigma*(1-y[1]*y[2])/10;
        return sig;
    }
}
data{
    int<lower=0> N;
    real<lower=0> dt;
    array[N] vector[2] y;
}
parameters{
    real<lower=0, upper=1> alpha;
    real<lower=0, upper=1> beta;
    real<lower=0, upper=1> gamma;
    real<lower=0, upper=1> sigma;
}
model{
for (n in 2:N) {
        y[n] ~ normal(y[n-1]+mu(y[n-1], alpha, beta, gamma)*dt, Sigma(y[n-1], sigma));
    }
}
""" 

This works totally fine and I can compile and fit the model. But if I now want to use a multivariate-distribution to draw samples fro, it does not work anymore. I specified the coavariance matrix as diagonal matrix such that mathematically it should be no difference to above, where I sample a vector from a vector of normal distributions

third_model = """
functions {
    vector mu(vector y, real alpha, real beta, real gamma) {
        vector[2] m;
        m = y*(beta-alpha-gamma-beta);
        return m;
    }
    matrix Sigma(vector y, real sigma) {
        matrix[2, 2] sig;
        sig[1, 1] = sigma*(1-y[1]*y[2])/10;
        sig[2, 2] = sigma*(1-y[1]*y[2])/10;
        return sig;
    }
}
data{
    int<lower=0> N;
    real<lower=0> dt;
    array[N] vector[2] y;
}
parameters{
    real<lower=0, upper=1> alpha;
    real<lower=0, upper=1> beta;
    real<lower=0, upper=1> gamma;
    real<lower=0, upper=1> sigma;
}
model{
for (n in 2:N) {
        y[n] ~ multi_normal(y[n-1]+mu(y[n-1], alpha, beta, gamma)*dt, Sigma(y[n-1], sigma));
    }
}
""" 

But this gives the same RuntimeError: Initialization failed. problem as above. So I guess there is a problem with the use of multi_normal, but I don't know which. In general the multi_normalshould be supported in STAN since they use it in their example of Gaussian processes in the documentation.

Answer 1

A covariance matrix for a multivariate normal distribution must be positive semi-definite and symmetric. In the context of a stan model, I expect that the covariance matrix must also be invertible, and hence must also be positive definite as well.

Looking at the definition of the sigma(vector x, real alpha, real beta) function, it looks like the covariance that you defined is neither positive-definite nor symmetric. I suspect that the initialization might be failing when first attempting to invert the covariance matrix, resulting in the observed Initialization failed error.