R中如何绕过.Machine$double.xmax上限 解决MH算法重尾采样溢出问题

首先明确：不需要尝试突破双精度浮点数的数值上限，这个1.797693e+308是IEEE 754双精度浮点数的硬上限，R、Python默认数值类型都遵循这个标准，强行用高精度算术包突破上限会让计算速度暴跌几十上百倍，完全没必要，这个溢出问题用对数域计算接受率就能从根源解决。

R端代码修复方案 #

1. 核心修改：对数域计算MH接受概率 #

你当前直接计算后验密度比值的写法是溢出的根源：重尾分布的后验密度值很容易超出双精度数上限变成Inf，导致比值出现NaN。
MH算法的接受判断完全可以在对数尺度完成，二者判断结果完全等价（对数是单调递增函数，不改变不等号方向）：

• 原始接受率：$a = \frac{p(\theta_{propose})}{p(\theta_{current})}$
• 对数接受率：$\log(a) = \log p(\theta_{propose}) - \log p(\theta_{current})$
• 判断规则：如果$\log(a) \geq 0$（即$a\geq1$）直接接受；否则生成均匀随机数$u$，判断$\log(u) < \log(a)$即可，不需要计算原始的$a$值。

你需要先把原来返回原始密度的target_1、target_2改成返回对数密度的版本（比如命名为log_target_1、log_target_2，只需要在你原来的密度计算代码最后加个log()转换即可，如果是自定义密度直接写对数形式数值稳定性更好），然后把采样循环里的接受判断部分替换成如下写法：

# 替换param1更新部分
proposed_param1 <- rnorm(1, current_param1, sd = sqrt(var_param1))
log_a <- log_target_param1(proposed_param1, current_param2) - log_target_param1(current_param1, current_param2)
if (log_a >= 0 || log(runif(1)) < log_a) {
  param_1[i] = proposed_param1 # 接受
} else {
  param_1[i] = current_param1  # 拒绝
}

# 替换param2更新部分
current_param1 = param_1[i]
proposed_param2 <- rnorm(1, current_param2, sd = sqrt(var_param2))
log_b <- log_target_param2(current_param1, proposed_param2) - log_target_param2(current_param1, current_param2)
if (log_b >=0 || log(runif(1)) < log_b) {
  param_2[i] = proposed_param2 # 接受
} else {
  param_2[i] = current_param2  # 拒绝
}

对数转换后，哪怕原始密度值达到$10^{1000}$，取对数后也仅为2300左右，完全落在双精度数的正常计算范围内，从根源上杜绝溢出。

2. 其他需要修复的代码问题 #

你当前代码里还有两个容易引发报错的点：

• 计算建议分布方差的optim调用中，初值写的是c(param1,param2)，但这个变量在MH函数内部没有定义，会错误调用全局环境中的同名变量，应该替换为startval作为初值。
• 海森矩阵求逆时建议给对角元加一个极小的正则项（比如1e-6），避免海森矩阵奇异时算出负的方差，导致正态分布采样出NaN；另外注意optim默认做函数最小化，如果传入的是对数似然/对数后验，需要加control = list(fnscale = -1)切换为最大化，否则海森矩阵符号错误，算出来的方差同样是负数。
  修正后的方差计算代码示例：

  opt_res <- optim(startval, LgLik, x=sample, hessian=T, control = list(fnscale = -1))
  hess_inv <- solve(-opt_res$hessian + diag(1e-6, 2)) # 负海森的逆是渐近协方差阵，加正则避免奇异
  var_param1 <- hess_inv[1,1]
  var_param2 <- hess_inv[2,2]
  

Python端同类问题说明 #

Python默认的浮点数类型同样是IEEE 754双精度float64，数值上限和R完全一致，如果你照搬原始密度比值的写法，一样会遇到Inf、NaN的溢出问题，解决方案和R完全相同：全部转对数域计算接受率即可，不需要用高精度计算库。

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