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.
