为什么 optim() 总是给我相同的结果？

Question

我正在尝试使用自举来计算参数 v 的 MLE 的置信区间。但在每次迭代中，optim() 给了我相同的结果。为什么会这样？

这是我的代码：

# Create log-likelihood function
dofloglik = function(v, x) { n <- length(x) 
                             loglik <- n*lgamma((v+1)/2) - n*lgamma(v/2)
                                       - (n/2)*log(v*pi) - ((v+1)/2)*sum(log(1+x^2/v))
                             return(-loglik) }

# Sample
y = c(0.97, 0.41, -0.13, -0.18, 1.69, 
       0.17, -0.46, -0.65, 2.66, -0.28, 
       2.30, -0.15, 0.24, 0.43, -0.54, 
       -1.55, 2.13, -0.48, 1.24, 1.05)

# Bootstrapping
n = length(y)
sims = 10
mle = numeric(sims)
for (i in 1:sims) { 
  bs = sample(y, n, replace=TRUE) 
  mle[i] = optim(1, dofloglik, x = bs, method = "CG")$par
  print(bs)
}
mle = sort(mle)
mle

输出是：

[1]  1.05  1.69 -0.13  2.66  2.13 -0.28 -0.54 -0.65  1.69  0.43 -0.18  0.17  1.24 -0.13 -1.55  0.41 -0.18  0.24 -0.54  2.13
[1] -0.13 -0.18  2.13 -1.55 -0.65 -0.28  1.24  2.13 -0.48  2.30 -0.28 -0.54 -0.13  0.41  1.24  1.24 -0.54 -0.18  2.66 -0.54
[1]  2.30  0.43 -1.55  1.05  2.13 -0.18  1.24  0.24 -0.13  2.30 -0.48  2.66 -0.13 -1.55 -0.54 -0.13  0.43 -0.13  0.24  2.30
[1]  1.24  1.24 -0.54  2.30 -0.46 -0.65  2.13  2.66 -0.46 -0.13 -0.18  2.66  0.24 -0.48 -1.55 -1.55  1.05  0.24  0.17 -0.46
[1]  0.97  2.66  2.13  0.24 -1.55  0.43  2.13  0.43 -0.18 -1.55 -0.46  1.69 -0.48  2.66 -0.18  2.30  2.66 -1.55 -0.54 -0.54
[1]  0.41  0.43  1.24  1.69 -0.13  0.41  1.69 -0.54 -0.28  2.13 -0.46  2.30  0.17  0.97  1.24  2.30  2.66  0.43  0.43  1.24
[1]  1.24  1.24  1.05  2.66  0.17  2.13  0.17 -0.13  2.30 -0.46 -1.55 -0.28  1.24  2.30 -0.48  0.24 -0.54  0.41 -0.65 -1.55
[1]  2.66 -0.54  1.05 -0.15  0.17 -1.55  0.41  1.69 -1.55 -1.55 -0.28 -0.46 -0.48 -0.13 -0.46  0.43  1.24  0.24 -0.46 -0.28
[1]  1.05  0.24  2.13  0.97  1.69  1.05  2.13 -0.15 -0.48 -1.55  1.05 -0.15  0.43 -0.13 -0.28  0.17  2.66 -0.15  1.24 -0.28
[1] -0.48 -0.18  0.24  2.30 -0.46 -0.54  0.43 -0.54  2.66 -0.48  2.66  0.24  2.13  0.97  1.05 -0.18  2.30 -0.13 -0.46  1.24
> mle = sort(mle)
> mle
[1] 47.16466 47.16466 47.16466 47.16466 47.16466 47.16466 47.16466 47.16466 47.16466 47.16466

Answer 1

> dofloglik = function(v, x) { n <- length(x) 
+                              loglik <- n*lgamma((v+1)/2) - n*lgamma(v/2)
+                                        - (n/2)*log(v*pi) - ((v+1)/2)*sum(log(1+x^2/v))
+                              return(-loglik) }
> 
> dofloglik(1, 1:4)
[1] 2.28946
> dofloglik(1, 2:5)
[1] 2.28946
> dofloglik(1, 3:6)
[1] 2.28946

这可能不是你想要的。

我怀疑不是您想要的部分是您将某些内容分成两行。除了每一行都是一个有效的命令，所以 R 自己执行第一行。

+                              loglik <- n*lgamma((v+1)/2) - n*lgamma(v/2)
+                                        - (n/2)*log(v*pi) - ((v+1)/2)*sum(log(1+x^2/v))

所以loglik <- n*lgamma((v+1)/2) - n*lgamma(v/2) 已运行且有效，因此它被存储到 loglik 中。第二行- (n/2)*log(v*pi) - ((v+1)/2)*sum(log(1+x^2/v)) 有效。它执行但不存储任何内容，而且由于我们不在顶层，它甚至不打印。

您要么需要确保第一行的末尾不会导致完整的语句（通过在第一行的末尾添加减号），要么执行其他操作以以不同的方式将其拆分。