PyStan RuntimeError:初始化失败

Question

我尝试使用STAN来推断一个SIR模型，它表示为一个三维SDE:该模型的代码编译时没有问题。

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.

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

更新：I成功地获得了一个较小的示例，与此相同，我降低了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));
    }
}
""" 

这是完全好的工作，我可以编译和符合模型。但是如果我现在想用多元分布来抽取样本，它就不再起作用了。我将coavariance矩阵指定为对角矩阵，因此在数学上，它应该与上面没有任何差别，其中我从正态分布的向量中采样一个向量。

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的使用有一个问题，但我不知道是哪一个。通常，STAN中应该支持multi_normal，因为它们在文档中的高斯过程示例中使用它。

Answer 1

多元正态分布的协方差矩阵必须是正的、半定的、对称的.在stan模型中，我认为协方差矩阵也必须是可逆的，因此也必须是正定的。

看看sigma(vector x, real alpha, real beta)函数的定义，你定义的协方差既不是正定的，也不是对称的。我怀疑，当第一次试图反转协方差矩阵时，初始化可能会失败，从而导致观察到的Initialization failed错误。