STAN, PyStan RuntimeError: Initialization failed

Question

我尝试使用 STAN 来推断表示为三维 SDE 的 SIR 模型：模型的代码编译没有问题。

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));
    }
}
"""

但是当我尝试使用

从模型中采样时

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

我收到以下错误消息：

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.

我不知道这个错误是从哪里来的，也不知道如何解决。我已经仔细检查了输入数据是否与模型匹配，但我似乎确实如此。希望有人有想法或已经遇到同样的问题。

更新： 我设法得到了一个具有相同问题的小例子我稍微降低了 SDE 的维度和复杂性。使用正态分布作为样本，均值和标准差作为向量提供

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));
    }
}
"""

这完全没问题，我可以编译和拟合模型。但是，如果我现在想使用多元分布来抽取样本，它就不再起作用了。我将协方差矩阵指定为对角矩阵，这样在数学上它应该与上面没有区别，我从正态分布的向量中采样向量

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));
    }
}
"""

但这会产生与上述相同的 RuntimeError: Initialization failed. 问题。所以我猜测multi_normal的使用有问题，但我不知道是哪个。一般来说，multi_normal 应该在 STAN 中得到支持，因为他们在文档中的高斯过程示例中使用了它。

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，因此也必须是正定的。

查看 sigma(vector x, real alpha, real beta) 函数的定义，您定义的协方差似乎既不是 positive-definite 也不是对称的。我怀疑第一次尝试反转协方差矩阵时初始化可能会失败，导致观察到 Initialization failed 错误。

STAN, PyStan RuntimeError: Initialization failed

STAN, PyStan RuntimeError: Initialization failed

python

stan

pystan