0

我正在尝试重现 Stan 手册 2.11 版第 15.2-15.3 节中的示例：

gpdemo.stan 
/*
 * Stan manual, ch. 15
 */

data {
  int<lower=1> N;
  real x[N];
}


transformed data {
  vector[N] mu;
  cov_matrix[N] Sigma;
  matrix[N, N] L;

  for (i in 1:N)
    mu[i] = 0;

  for (i in 1:N)
    for (j in 1:N)
      Sigma[i, j] = exp(-square(x[i] - x[j])) + i == j ? 0.1 : 0.0;
  L = cholesky_decompose(Sigma);
}


parameters {
  vector[N] z;
}


model {
  z ~ normal(0, 1);
}


generated quantities {
  vector[N] y;
  y = mu + L * z;
}
gpdemo.R 
library('rstan')

x <- -50:50 / 10

data <- list(
  N = length(x),
  x = x
)

sim <- stan("gp-demo.stan", "gp_demo", data = data, iter = 200, chains = 3)

但斯坦抱怨说这Sigma不是肯定的：

Error : Exception thrown at line 23: cholesky_decompose: Matrix m is not positive definite
failed to create the sampler; sampling not done

事实上它不是：

x_dist <- as.matrix(dist(x))
x_kernel <- exp(-x_dist^2)
eigen(x_kernel)$values
# [1] FALSE

但问题似乎是数值精度之一，而不是其他问题：

summary(eigen(x_kernel)$values)
#      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
#  0.000000  0.000000  0.000000  1.000000  0.000049 17.360000 
min(eigen(x_kernel)$values)
# [1] -1.125251e-15

显然这个例子在某人的机器上运行没有数字问题（否则它不会出现在手册中，更不用说在互联网上的博客上转载了）。这是一个非常简单的例子，所以我想知道我是否只是在某个地方犯了错误。

对于它的价值，我能够通过更明智地循环构造来使其工作Sigma：

for (i in 1:N) Sigma[i, i] = 1.1;
for (i in 1:(N-1))
  for (j in (i+1):N) {
    Sigma[i, j] = exp(-square(x[i] - x[j]));
    Sigma[j, i] = Sigma[i, j];
}
4

0 回答 0