floating-point - 数値精度が原因でスタンの基本的な例を再現できない

Question

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しかし、Stan はそれが正定値ではないことを訴えています。

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];
}