Почему optim () всегда дает мне один и тот же результат? - PullRequest
Новая
       91

Почему optim () всегда дает мне один и тот же результат?

0 голосов
/ 10 ноября 2019

Я пытаюсь использовать начальную загрузку для вычисления доверительного интервала MLE параметра v. Но на каждой итерации 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

1 Ответ

1 голос
/ 10 ноября 2019 
> 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)) действительна. Он выполняется, но ничего не хранит, и, поскольку мы не на верхнем уровне, он даже не печатает.

Вам также нужно убедиться, что конец вашей первой строки не приведет к завершениюоператор (добавив знак минус в конец первой строки) или сделайте что-нибудь еще, чтобы разделить это по-другому.

...