Stan拟合联合似然模型采样异常致bridgesampling报错排查

咱们先把问题拆解清楚：你遇到的bridge_sampler报错，本质是Stan模型采样时的发散、低BFMI等问题导致后验样本不可靠，甚至出现大量log_prob计算为无穷的情况，直接让桥抽样无法进行。下面一步步分析原因和解决办法：

1. 先搞定Stan的采样警告（发散+低BFMI） #

你的模型是混合了连续结果（y的正态模型）和离散结果（x1的probit模型）的分层模型，这类模型天生容易出现数值不稳定，尤其是当参数尺度不匹配、先验设置不合理时。

核心优化点： #

• 标准化所有预测变量：这是最容易被忽略但效果显著的步骤。把x、G、X1的每一列都标准化到均值0、方差1（R里用scale()函数就行）。如果预测变量尺度差异太大，会导致参数（比如gamma3）的后验分布极端狭窄或宽泛，HMC采样器根本“走不动”，直接引发发散和低BFMI。
• 调整采样控制参数：你已经把adapt_delta调到0.99了，可以再试调到0.999，同时把max_treedepth拉到18。更高的adapt_delta会让采样器更谨慎地探索参数空间，减少发散；更大的树深度能处理更复杂的后验几何形状。另外，可以适当增加迭代次数，比如iter=60000、warmup=30000，给采样器足够的热身时间。
• 检查先验设置：
  • 你的gamma3用的是N(0,1)先验，如果X1没标准化，这个尺度可能太大，容易让X1*gamma3的绝对值变得极大，导致Phi(X1*gamma3)接近0或1，引发数值问题。可以试试收缩更强的先验，比如N(0, 0.8)或N(0, 0.5)。
  • 检查sigma的InvGamma(4,3)先验：它的均值是1，标准差约0.707，如果你的y样本方差和这个差异很大，会导致sigma的采样异常。可以根据y的样本方差调整先验，比如用InvGamma(2,1)（均值1，标准差1），或者更宽松的先验。

2. 解决bridge_sampler的log_prob无穷问题 #

报错里提到“2744 of the 45000 log_prob() evaluations on the proposal draws produced -Inf/Inf”，这是因为当Phi(X1*gamma3)接近0时，如果对应的x1=1，log(p)就会变成-Inf；反之Phi接近1时x1=0，log(1-p)也会变成-Inf。这说明你的后验样本里有大量gamma3的极端取值，直接导致log概率计算崩溃。

针对性解决： #

• 改用潜变量参数化的probit模型：直接用潜变量替代Phi的计算，能避免数值溢出。修改你的Stan模型代码：

  stanmodelH1 = 'data {
      int<lower=0> n;
      int<lower=0> k;
      vector[n] y;
      int x1[n];
      real<lower=0> a;
      real<lower=0> b;
      matrix[n, k] G;
      matrix[n, k+1] X1;
      matrix[n, 2] x;
  }
  
  parameters {
      vector[2] beta;
      vector[k+1] gamma3;
      vector[k] gamma1;
      real<lower=0> sigma;
      vector[n] z1; // 新增潜变量
  }
  
  model {
      target += inv_gamma_lpdf(sigma | a, b);
      target += normal_lpdf(gamma1 | 0, sqrt(1.1));
      target += normal_lpdf(gamma3 | 0, 1); // 标准化X1后可以调整这个尺度
      target += normal_lpdf(beta | 0, sqrt(1.2));
      // 用潜变量建模probit
      z1 ~ normal(X1*gamma3, 1);
      x1 ~ bernoulli(Phi(z1));
      // 原有的正态模型
      target += normal_lpdf(y | x*beta + G*gamma1, sqrt(sigma));
  }
  '
  

  潜变量z1是连续的，HMC采样器处理起来更稳定，避免了直接计算Phi带来的数值极值问题。

3. 验证采样有效性 #

在运行bridge_sampler之前，一定要先确认Stan采样是正常的：

• 用check_hmc_diagnostics(stanfitmodelH1)检查，确保没有发散、BFMI值在合理范围（通常>0.2）。
• 用pairs(stanfitmodelH1, pars = c("beta", "gamma1", "gamma3", "sigma"))查看参数对的分布，如果出现明显的“漏斗形状”，说明参数尺度还是有问题，得继续调整先验或标准化。
• 查看gamma3的后验分布（plot(stanfitmodelH1, pars = "gamma3")），如果有极端大的取值，说明先验还是不够收缩，或者预测变量没标准化。

只要Stan采样的警告都消失了，bridge_sampler应该就能正常运行了。

内容的提问来源于stack exchange，提问作者geraltofrivia