Максимальное правдоподобие для многопараметрического нелинейного нормального распределения

Question

Я пытаюсь написать код для определения g, sigma и лямбда, используя максимальное значение для следующих Function Equation. Мне нужно было ограничить интервалы, чтобы гарантировать, что уравнение не берет журнал отрицательного числа, однако оно приходит с ошибкой «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") 

Answer 1

Вам нужны ограничения 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