Блог

Я пытаюсь написать код для определения g, сигма и лямбда, используя максимальное правдоподобие для следующего функционального уравнения. Мне нужно было ограничить интервалы, чтобы убедиться, что уравнение не принимает логарифм отрицательного числа, однако оно выдает ошибку «L-BFGS-B нужны конечные значения». Что я могу сделать, чтобы исправить это?

 lnQs<- c(0.4452211,  3.2828926,  3.0752400,  2.6613305,  5.8122312,  0.5629881,  0.8445112,  4.4336806,  3.8253957,  0.9336889,  1.5188934,  4.3915304, 2.5368227,  0.3729370,  2.3683679,  2.3555777,  0.5985054,  0.6240360,  0.9462143,  0.5440311,  2.6390102,  0.4921728,  0.2971820,  0.1939826, 0.5621709,  0.7839881,  2.2367834, 10.4101839,  5.8010886,  6.0974008)
Rf<- c(0.0484, 0.0537, 0.0598, 0.0590, 0.0590, 0.0571, 0.0497, 0.0497, 0.0492, 0.0588, 0.0586, 0.0492, 0.0468, 0.0480, 0.0405, 0.0405, 0.0488, 0.0452, 0.0675, 0.0395, 0.0395, 0.0385, 0.0400, 0.0432, 0.0394, 0.0397, 0.0397, 0.0407, 0.0436, 0.0436)
S<- c(0.8,  2.3,  2.2,  3.3,  4.9,  8.7,  0.8,  0.9,  1.8,  2.3,  1.3,  1.9,  5.7,  9.3,  4.9, 18.7,  3.6,  2.4, 15.1, 10.3,  0.8,  2.9, 12.0,  0.8,  9.9,  1.3,  8.9, 12.3,  4.2,  4.2)
T<- c(1.0,  5.0,  5.0,  7.0,  5.1, 21.0, 14.1,  5.0,  5.0, 12.0,  7.0,  3.0,  7.0, 21.0,  7.0, 21.0, 21.0, 21.0, 21.0, 20.0,  5.0, 21.0, 21.0, 21.1, 21.0, 14.0,  12.0, 14.0,  5.0,  5.0)

LL<- function(g,sigma,lambda){
  R=(dnorm(lnQs,(log(1 Rf lambda)-sigma^2/2)*S-log(((1 Rf lambda)/(1 g))^T-1),sigma*S^0.5))
#
    -sum(log(R), log=TRUE)
}
fit<-mle(minuslogl=LL, start=list(g=.05, sigma=.2, lambda=.1), method = "L-BFGS-B",lower=c(0,0,0.0915),upper=c(0.13,Inf,Inf))
summary(fit)
#Criteria required 
#   lambda>-(1 Rf)  - easily done with restriction lambda>0
#   lambda>(g-Rf)   - NOT SURE HOW TO DEAL WITH lowest Rf=0.0385, tried putting upper limit on g and lower limit on lambda for now
#   sigma>0         - easily done with restriction sigma>0
#   Problem that L-BFGS-B needs finite values offit<-mle(minuslogl=LL, start=list(g=.077, sigma=.256, lambda=.110), method = "BFGS")
 

Нужны ли вам пределы Inf для ваших сигма- и лямбда-параметров для MLE? например, следующее работает без ошибок (хотя лямбда-оценка довольно плохая с высокой ошибкой std):

 set.seed(123)
T<-0.5
S<-1
Rf<-2
g <-.1 
sigma <- .5
lambda <- 1
lnQs <- rnorm(100,(log(1 Rf lambda)-sigma^2/2)*S-log(((1 Rf lambda)/(1 g))^T-1),sigma*S^0.5)

# negative ll fn  
LL<- function(g,sigma,lambda){
  R <- dnorm(lnQs,(log(1 Rf lambda)-sigma^2/2)*S-log(((1 Rf lambda)/(1 g))^T-1),sigma*S^0.5)
  -sum(log(R))
}
fit<-mle(minuslogl=LL, start=list(g=.05, sigma=.2, lambda=.1), method = "L-BFGS-B",lower=c(0,0,0.0915),upper=c(0.13,5,5))
summary(fit)
Maximum likelihood estimation

Call:
mle(minuslogl = LL, start = list(g = 0.05, sigma = 0.2, lambda = 0.1), 
    method = "L-BFGS-B", lower = c(0, 0, 0.0915), upper = c(0.13, 
        5, 5))

Coefficients:
       Estimate Std. Error
g       0.10583  16.869598
sigma   0.45412   0.032111
lambda  0.39076 363.867321
 

С предоставленными вами данными LL вычисляется до бесконечности во многих точках с множеством разных значений ваших параметров, например, используя простой поиск по сетке в [0,1] x [0,1] x [0,1], вы можете найти следующие точки, где LL вычисляется до бесконечности (а иногда и NAN), чтовот почему L-BFGS-G терпит неудачу:

 g <- seq(0,1,length=10)
sigma <- seq(0,1,length=10)
lambda <- seq(0,1,length=10)
# grid search
for (i in 1:10) {
  for (j in 1:10) {
    for (k in 1:10) {
      if (LL(g[i],sigma[j],lambda[k]) == Inf) {
        print(paste(g[i],sigma[j],lambda[k]))
      }
    }
  }
}
 

Некоторые точки, в которых LL оценивается как Inf

 [1] "0 0 0"
[1] "0 0 0.111111111111111"
[1] "0 0 0.222222222222222"
[1] "0 0 0.333333333333333"
[1] "0 0 0.444444444444444"
[1] "0 0 0.555555555555556"
[1] "0 0 0.666666666666667"
[1] "0 0 0.777777777777778"
[1] "0 0 0.888888888888889"
[1] "0 0 1"
[1] "0 0.111111111111111 0.111111111111111"
[1] "0 0.111111111111111 0.222222222222222"
[1] "0 0.111111111111111 0.333333333333333"
[1] "0 0.111111111111111 0.444444444444444"
 

Существует простой способ, который может решить проблему:

 LL<- function(g,sigma,lambda){
      R=(dnorm(lnQs,(log(1 Rf lambda)-sigma^2/2)*S-log(((1 Rf lambda)/(1 g))^T-1),sigma*S^0.5))
      #
      ll <- -sum(log(R), log=TRUE)
      ifelse(ll == Inf, 9999999, ll) # return a large enough number if Inf
}

fit<-mle(minuslogl=LL, start=list(g=.05, sigma=.2, lambda=.1), method = "L-BFGS-B",lower=c(0,0,0.0915),upper=c(0.13,Inf,Inf))
summary(fit)


Maximum likelihood estimation

Call:
mle(minuslogl = LL, start = list(g = 0.05, sigma = 0.2, lambda = 0.1), 
    method = "L-BFGS-B", lower = c(0, 0, 0.0915), upper = c(0.13, 
        Inf, Inf))

Coefficients:
        Estimate Std. Error
g      0.1249300 0.07295480
sigma  0.8795551 0.07349924
lambda 0.0915000 0.07459289

-2 log L: 160.1491