用于拟合截断正态分布的相同代码不同输出

Question

我是运行同一段代码，具有相同的种子、相同的包版本、相同的 R 版本，在 3 个不同的系统上：1) 我的电脑 2) linux 集群和3) R snippets

packageVersion("truncnorm")
packageVersion("MASS")
set.seed(42)
fit<-NULL
x <- c(0.0916, 0.0084, 0.0442, 0.6254, 0.2021, 0.0135, 0.0259,
       0.1557,0.0191, 0.3575, 0.1843, 0.1792, 0.0476, 0.0765, 
       0.0356, 0.0039, 0.1714, 0.1222, 0.2872, 0.395, 0.3334,
       0.2223, 0.0096, 0.0436, 0.207)
mu0 <- mean(x)
sigma0 <- stats::sd(x)
fit <- MASS::fitdistr(x, densfun = function(xx, mu, sigma) {
    truncnorm::dtruncnorm(xx, a = 0, b = 1, mean = mu, sd = sigma)
}, 
   start = list(mu = mu0, sigma = sigma0), 
   lower = list(mu = -Inf, sigma = 0.05), 
   upper = list(mu = Inf, sigma = Inf))
print(fit)

在我的电脑上，fit 显示为 NULL，而在其他 2 个系统中，该模型确实拟合成功。知道这怎么可能吗？

P.S.: 我系统的问题是

Error in MASS::fitdistr(x, densfun = function(xx, mu, sigma) { : optimization failed

如果我稍微更改数据，例如从数据中删除 0.0084（这是数据中的第二个数字），则模型适合。在所有 3 个系统中给我相同的输出。

这是来自我自己系统的sessionInfo()：

R version 3.6.0 (2019-04-26) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 19043)

Matrix products: default

locale: [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252 LC_MONETARY=English_United States.1252 [4] LC_NUMERIC=C
LC_TIME=English_United States.1252

attached base packages: [1] stats graphics grDevices utils
datasets methods base

other attached packages: [1] truncnorm_1.0-8 MASS_7.3-51.4

loaded via a namespace (and not attached): [1] ks_1.13.2
compiler_3.6.0 Matrix_1.2-17 mclust_5.4.7 tools_3.6.0
simIReff_1.0 [7] mvtnorm_1.1-3 KernSmooth_2.23-15 grid_3.6.0 pracma_2.3.3 lattice_0.20-38

Answer 1

这原来是一个数值 unstable/sensitive 问题。

如果您打开 debug(MASS::fitdistr) 并单步执行，最终您会到达一行

 if (res$convergence > 0L) stop("optimization failed")

如果你此时打印出res的值，你会得到（略有缩写）：

$par
       mu     sigma
-6.411168  1.022651
$value
[1] -21.72969
$counts
function gradient
      81       81
$convergence
[1] 52
$message
[1] "ERROR: ABNORMAL_TERMINATION_IN_LNSRCH"

换句话说，L-BFGS-B 优化器（使用它是因为您指定了边界 - 它非常挑剔）认为存在问题，因此 fitdistr 会抛出错误。据我所知，没有办法告诉 fitdistr“无论如何都给我答案”。

我尝试了很多不同的方法（稍微扰乱了起始条件，即 mu + 1e-3, sigma0 + 1e-3；删除边界，以便 fitdistr 使用更强大的 Nelder-Mead 优化器） .绘制 log10(1e-4 + neg log likelihood)（以便我们可以看到与最小负数 log-likelihood 的微小差异）给出以下图像（下面的代码）：

[红色：Linux，绿色：Windows/convergence失败，蓝色：Windows/perturbed开始，青色：Nelder-Mead]

这些是分布的相应拟合：

如您所见（或者可能看不到！），所有的配合基本上是相同的。如果计算负值 log-likelihood，您会发现它们的差异小于 0.001 个单位 [即可以忽略不计]。（您也可以看出这一点，因为第一张图像中的所有点都位于 log10(difference) = -3 轮廓内。）

所以答案之间的差异并不重要，只是出现错误的烦恼。您可以 (1) 使用 while 循环 + try() 稍微扰动起始值，直到得到答案； (2) 降低界限以允许 Nelder-Mead 工作：(3) 使用 bbmle 或其他一些工具，让您对优化过程有更多 robust/defensive 的了解 ...

nllfun <- function(mu, sigma) {
  -sum(log(dtruncnorm(x, a = 0, b = 1, mean = mu, sd = sigma)))
}
library(emdbook)
library(truncnorm)
p1 <- c(-7.02938981, 1.06779942) ## Linux
p2 <- c(-6.411, 1.022651)  ## Windows (convergence error)
p3 <- c(-6.587645, 1.0359466) ## Windows (perturbed start)
p4 <- c(-5.9937989, 0.9901366) ## Windows (Nelder-Mead/no bounds)
cc <- curve3d(nllfun(x,y), xlim = c(-7.1, -5.98), ylim = c(0.98, 1.07),
              n = c(101, 101), sys3d = "none")

image(cc$x, cc$y, log10(cc$z-min(cc$z) + 1e-4))
contour(cc$x, cc$y, log10(cc$z-min(cc$z) + 1e-4), add = TRUE)
points(p1[1], p1[2], pch = 16, col = 2)
points(p2[1], p2[2], pch = 17, col = 3)
points(p3[1], p3[2], pch = 18, col = 4)
points(p4[1], p4[2], pch = 18, col = 5)

hist(x, freq=FALSE)
curve(dtruncnorm(x, a=0, b=1, mean=p1[1], sd = p1[2]), col = 2, add=TRUE)
curve(dtruncnorm(x, a=0, b=1, mean=p2[1], sd = p2[2]), col = 3, add=TRUE)
curve(dtruncnorm(x, a=0, b=1, mean=p3[1], sd = p3[2]), col = 4, add=TRUE)
curve(dtruncnorm(x, a=0, b=1, mean=p4[1], sd = p4[2]), col = 5, add=TRUE)

用于拟合截断正态分布的相同代码不同输出

Same code different outputs for fitting Truncated Normal Distribution

statistics

r

normal-distribution