如题，求助大家一个关于 Stan 建模的问题，下面是具体数据、代码和模型

数据

head(rongelap)

     cX    cY counts time
1 -6050 -3270     75  300
2 -6050 -3165    371  300
3 -5925 -3320   1931  300
4 -5925 -3165   4357  300
5 -5800 -3350   2114  300
6 -5800 -3165   2318  300

cX 和 cY 是横、纵坐标，time 是时间，counts 是时间 time 内的次数。

统计模型

$$Y(x)$$ 和 $$u(x)$$ 是已知的观测值， $x \in \mathbb{R}^2$ 代表二维坐标。

$\log(\lambda(x) )= \alpha + S(x)$
$Y(x) \sim \mathrm{Poisson}(u(x)\lambda(x) )$

其中 $$S(x)$$ 是一个高斯过程，其自协方差函数如下

$\mathsf{Cov}(S(x_i), S(x_j)) = \sigma^2 \exp\big(- \frac{\|x_i - x_j \|_2}{\phi}\big)$

对模型来说， $\sigma$ 和 $\phi$ 是超参数， $\phi$ 是大于 0 的。更多详细的关于此模型的介绍，可以见这里

运行代码

# 加载 cmdstanr 包
library(cmdstanr)
# 设置 cmdstan 安装位置
set_cmdstan_path(path = "/opt/cmdstan/cmdstan-2.31.0/")

rongelap <- readRDS(file = "data/rongelap.rds")

nchains <- 4 # 4 条迭代链
# 准备数据
rongelap_poisson_d <- list(
  N = nrow(rongelap), # 观测记录的条数
  D = 2, # 坐标维度
  X = rongelap[, c("cX", "cY")], # N x 2 矩阵
  counts = rongelap$counts, # N 向量
  time = rongelap$time
)

# 编译模型
mod_rongelap_poisson <- cmdstan_model(
  stan_file = "code/rongelap_gp_poisson.stan",
  compile = TRUE,
  cpp_options = list(stan_threads = TRUE)
)

# 采样拟合模型
fit_rongelap_poisson <- mod_rongelap_poisson$sample(
  data = rongelap_poisson_d, # 观测数据
  # init = inits_data, # 迭代初值
  iter_warmup = 1000, # 每条链预处理迭代次数
  iter_sampling = 4000, # 每条链总迭代次数
  chains = nchains, # 马尔科夫链的数目
  parallel_chains = 4, # 指定 CPU 核心数，可以给每条链分配一个
  threads_per_chain = 1, # 每条链设置一个线程
  show_messages = FALSE, # 不显示迭代的中间过程
  refresh = 0, # 不显示采样的进度
  seed = 20232023 # 设置随机数种子，不要使用 set.seed() 函数
)

fit_rongelap_poisson$diagnostic_summary()

$num_divergent
[1] 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0

$ebfmi
[1] 0.7486386 0.7157464 0.6727842 0.7087742

fit_rongelap_poisson$summary(variables = c("lp__", "sigma", "phi", "alpha"))

# A tibble: 4 × 10
  variable       mean  median      sd     mad      q5     q95  rhat ess_bulk ess_tail
  <chr>         <num>   <num>   <num>   <num>   <num>   <num> <num>    <num>    <num>
1 lp__        3.40e+6 3.40e+6 13.2    14.8    3.40e+6 3.40e+6  1.04     92.5     511.
2 sigma       4.79e-1 4.78e-1  0.0296  0.0292 4.32e-1 5.29e-1  1.08     42.7     144.
3 phi         1.23e+0 1.06e+0  0.676   0.465  5.45e-1 2.49e+0  1.00   8982.     8125.
4 alpha       1.94e+0 1.94e+0  0.0384  0.0360 1.88e+0 2.01e+0  1.10     47.0     105.

ess_bulk 和 ess_bulk 都太少，特别是 alpha ，rhat 明显大于 1。一些诊断信息如下：

fit_rongelap_poisson$cmdstan_diagnose()

Processing csv files: /var/folders/cv/7n2v46gx65bbdtv7ctx7jtmh0000gn/T/RtmptycRum/rongelap_gp_poisson-202303291723-1-134385.csv, /var/folders/cv/7n2v46gx65bbdtv7ctx7jtmh0000gn/T/RtmptycRum/rongelap_gp_poisson-202303291723-2-134385.csv, /var/folders/cv/7n2v46gx65bbdtv7ctx7jtmh0000gn/T/RtmptycRum/rongelap_gp_poisson-202303291723-3-134385.csv, /var/folders/cv/7n2v46gx65bbdtv7ctx7jtmh0000gn/T/RtmptycRum/rongelap_gp_poisson-202303291723-4-134385.csv

Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
No divergent transitions found.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Effective sample size satisfactory.

The following parameters had split R-hat greater than 1.05:
  alpha, sigma, eta[1], eta[2], eta[3], eta[4], eta[5], eta[6], eta[7], eta[8], eta[9], eta[10], eta[11], eta[12], eta[13], eta[14], eta[15], eta[16], eta[18], eta[19], eta[20], eta[21], eta[22], eta[23], eta[25], eta[26], eta[27], eta[28], eta[29], eta[30], eta[31], eta[32], eta[33], eta[34], eta[35], eta[36], eta[37], eta[38], eta[40], eta[41], eta[42], eta[43], eta[45], eta[47], eta[48], eta[49], eta[50], eta[51], eta[52], eta[53], eta[54], eta[56], eta[57], eta[58], eta[59], eta[60], eta[61], eta[62], eta[63], eta[64], eta[65], eta[66], eta[67], eta[68], eta[69], eta[70], eta[71], eta[72], eta[73], eta[74], eta[75], eta[76], eta[77], eta[78], eta[79], eta[80], eta[81], eta[82], eta[83], eta[84], eta[85], eta[86], eta[87], eta[88], eta[89], eta[90], eta[91], eta[92], eta[93], eta[94], eta[95], eta[96], eta[97], eta[98], eta[99], eta[100], eta[101], eta[102], eta[103], eta[104], eta[105], eta[106], eta[107], eta[108], eta[109], eta[110], eta[111], eta[112], eta[113], eta[114], eta[115], eta[116], eta[117], eta[118], eta[119], eta[120], eta[121], eta[122], eta[123], eta[124], eta[125], eta[126], eta[127], eta[128], eta[129], eta[130], eta[131], eta[132], eta[133], eta[134], eta[135], eta[136], eta[137], eta[138], eta[139], eta[140], eta[141], eta[142], eta[143], eta[144], eta[145], eta[146], eta[147], eta[148], eta[149], eta[150], eta[151], eta[152], eta[153], eta[154], eta[155], eta[156], eta[157], K[1,1], K[2,2], K[3,3], K[4,4], K[5,5], K[6,6], K[7,7], K[8,8], K[9,9], K[10,10], K[11,11], K[12,12], K[13,13], K[14,14], K[15,15], K[16,16], K[17,17], K[18,18], K[19,19], K[20,20], K[21,21], K[22,22], K[23,23], K[24,24], K[25,25], K[26,26], K[27,27], K[28,28], K[29,29], K[30,30], K[31,31], K[32,32], K[33,33], K[34,34], K[35,35], K[36,36], K[37,37], K[38,38], K[39,39], K[40,40], K[41,41], K[42,42], K[43,43], K[44,44], K[45,45], K[46,46], K[47,47], K[48,48], K[49,49], K[50,50], K[51,51], K[52,52], K[53,53], K[54,54], K[55,55], K[56,56], K[57,57], K[58,58], K[59,59], K[60,60], K[61,61], K[62,62], K[63,63], K[64,64], K[65,65], K[66,66], K[67,67], K[68,68], K[69,69], K[70,70], K[71,71], K[72,72], K[73,73], K[74,74], K[75,75], K[76,76], K[77,77], K[78,78], K[79,79], K[80,80], K[81,81], K[82,82], K[83,83], K[84,84], K[85,85], K[86,86], K[87,87], K[88,88], K[89,89], K[90,90], K[91,91], K[92,92], K[93,93], K[94,94], K[95,95], K[96,96], K[97,97], K[98,98], K[99,99], K[100,100], K[101,101], K[102,102], K[103,103], K[104,104], K[105,105], K[106,106], K[107,107], K[108,108], K[109,109], K[110,110], K[111,111], K[112,112], K[113,113], K[114,114], K[115,115], K[116,116], K[117,117], K[118,118], K[119,119], K[120,120], K[121,121], K[122,122], K[123,123], K[124,124], K[125,125], K[126,126], K[127,127], K[128,128], K[129,129], K[130,130], K[131,131], K[132,132], K[133,133], K[134,134], K[135,135], K[136,136], K[137,137], K[138,138], K[139,139], K[140,140], K[141,141], K[142,142], K[143,143], K[144,144], K[145,145], K[146,146], K[147,147], K[148,148], K[149,149], K[150,150], K[151,151], K[152,152], K[153,153], K[154,154], K[155,155], K[156,156], K[157,157], L_K[1,1], L_K[2,2], L_K[3,3], L_K[4,4], L_K[5,5], L_K[6,6], L_K[7,7], L_K[8,8], L_K[9,9], L_K[10,10], L_K[11,11], L_K[12,12], L_K[13,13], L_K[14,14], L_K[15,15], L_K[16,16], L_K[17,17], L_K[18,18], L_K[19,19], L_K[20,20], L_K[21,21], L_K[22,22], L_K[23,23], L_K[24,24], L_K[25,25], L_K[26,26], L_K[27,27], L_K[28,28], L_K[29,29], L_K[30,30], L_K[31,31], L_K[32,32], L_K[33,33], L_K[34,34], L_K[35,35], L_K[36,36], L_K[37,37], L_K[38,38], L_K[39,39], L_K[40,40], L_K[41,41], L_K[42,42], L_K[43,43], L_K[44,44], L_K[45,45], L_K[46,46], L_K[47,47], L_K[48,48], L_K[49,49], L_K[50,50], L_K[51,51], L_K[52,52], L_K[53,53], L_K[54,54], L_K[55,55], L_K[56,56], L_K[57,57], L_K[58,58], L_K[59,59], L_K[60,60], L_K[61,61], L_K[62,62], L_K[63,63], L_K[64,64], L_K[65,65], L_K[66,66], L_K[67,67], L_K[68,68], L_K[69,69], L_K[70,70], L_K[71,71], L_K[72,72], L_K[73,73], L_K[74,74], L_K[75,75], L_K[76,76], L_K[77,77], L_K[78,78], L_K[79,79], L_K[80,80], L_K[81,81], L_K[82,82], L_K[83,83], L_K[84,84], L_K[85,85], L_K[86,86], L_K[87,87], L_K[88,88], L_K[89,89], L_K[90,90], L_K[91,91], L_K[92,92], L_K[93,93], L_K[94,94], L_K[95,95], L_K[96,96], L_K[97,97], L_K[98,98], L_K[99,99], L_K[100,100], L_K[101,101], L_K[102,102], L_K[103,103], L_K[104,104], L_K[105,105], L_K[106,106], L_K[107,107], L_K[108,108], L_K[109,109], L_K[110,110], L_K[111,111], L_K[112,112], L_K[113,113], L_K[114,114], L_K[115,115], L_K[116,116], L_K[117,117], L_K[118,118], L_K[119,119], L_K[120,120], L_K[121,121], L_K[122,122], L_K[123,123], L_K[124,124], L_K[125,125], L_K[126,126], L_K[127,127], L_K[128,128], L_K[129,129], L_K[130,130], L_K[131,131], L_K[132,132], L_K[133,133], L_K[134,134], L_K[135,135], L_K[136,136], L_K[137,137], L_K[138,138], L_K[139,139], L_K[140,140], L_K[141,141], L_K[142,142], L_K[143,143], L_K[144,144], L_K[145,145], L_K[146,146], L_K[147,147], L_K[148,148], L_K[149,149], L_K[150,150], L_K[151,151], L_K[152,152], L_K[153,153], L_K[154,154], L_K[155,155], L_K[156,156], L_K[157,157], f[3], f[4], f[5], f[6], f[7], f[8], f[9], f[10], f[11], f[12], f[13], f[14], f[15], f[16], f[17], f[18], f[19], f[20], f[21], f[22], f[23], f[24], f[25], f[26], f[27], f[28], f[29], f[30], f[31], f[32], f[33], f[34], f[35], f[36], f[37], f[38], f[39], f[40], f[41], f[42], f[43], f[44], f[45], f[46], f[47], f[48], f[49], f[50], f[51], f[52], f[53], f[54], f[55], f[56], f[57], f[58], f[59], f[60], f[61], f[62], f[63], f[64], f[65], f[66], f[67], f[68], f[69], f[70], f[71], f[72], f[73], f[74], f[75], f[76], f[77], f[78], f[79], f[80], f[81], f[82], f[83], f[84], f[85], f[87], f[88], f[89], f[90], f[91], f[92], f[93], f[94], f[95], f[96], f[97], f[98], f[99], f[100], f[101], f[102], f[103], f[104], f[105], f[106], f[107], f[108], f[109], f[110], f[112], f[113], f[114], f[115], f[116], f[117], f[118], f[119], f[120], f[121], f[122], f[123], f[124], f[125], f[126], f[127], f[128], f[129], f[130], f[131], f[132], f[133], f[134], f[135], f[136], f[137], f[138], f[139], f[140], f[141], f[142], f[143], f[144], f[145], f[147], f[148], f[149], f[150], f[151], f[152], f[153], f[154], f[155], f[156], f[157]
Such high values indicate incomplete mixing and biased estimation.
You should consider regularizating your model with additional prior information or a more effective parameterization.

Processing complete.

我找了一些材料，但是目前的代码是我已经尽力优化的了，cholesky_decompose 分解做重参数化，代码尽可能向量化，只用 Stan 内置的函数，如指数型自协方差函数 gp_exponential_cov。

目前，运行上述代码约 1 小时，耗时也多，可见采样效率低下。不知道还要如何选择更加合适的先验分布？是不是真实数据本身严重影响先验的设置？

附件

rongelap.rds 数据集在这

rongelap_gp_poisson.stan 文件内容如下：

data {
  int<lower=1> N;
  int<lower=1> D;
  array[N] vector[D] X;
  array[N] int<lower = 0> counts;
  vector[N] time;
}
transformed data {
  real delta = 1e-12; // 防止矩阵奇异保证 cholesky 分解
  vector[N] log_T = log(time);
}
parameters {
  real alpha;
  real<lower=0> sigma;
  real<lower=0> phi;
  vector[N] eta;
}
transformed parameters {
    matrix[N, N] K = gp_exponential_cov(X, sigma, phi) + diag_matrix(rep_vector(delta, N));
    matrix[N, N] L_K = cholesky_decompose(K);
    vector[N] f = L_K * eta; // f 是高斯过程的一个实现
}
model {
  alpha ~ std_normal();
  sigma ~ std_normal();
  phi ~ inv_gamma(5, 5);
  eta ~ std_normal();

  counts ~ poisson_log(log_T + alpha + f);
}

参数 alpha 是标准正态先验，sigma 是非负标准正态先验，phi 是逆伽马先验，目前都是一般推荐的先验设置。~~感觉贝叶斯统计里最麻烦之一就是先验的选择。~~

运行环境

R 包 CmdStanR 版本 0.5.3
CmdStan 版本 2.31.0

参考材料

Fitting a Gaussian process https://mc-stan.org/docs/stan-users-guide/fit-gp.html 模拟结果不错，但是在我的这个实际数据上，不太行。
Stan Geostatistical Model https://rpubs.com/NickClark47/stan_geostatistical
关于 rongelap 数据的介绍见 https://bookdown.org/xiangyun/data-analysis-in-action/analyze-spatial-data.html