Я попробовал следующий код, чтобы найти MLE для четырех параметров: mu, sigma, alpha, lambda:
beta<-0.51#constant value which is needed in the likelihood
Следующие значения являются начальными значениями для mu, sigma, alpha и lambda:
mu=4;sig=0.01;alph=0.2;lam=1.5
composite.trapezoid <- function(f, a, b, n) {
if (is.function(f) == FALSE) {
stop('f must be a function with one variable')}
h <- (b - a) / n
j <- 1:n - 1
xj <- a + j * h
approx <- (h / 2) * (f(a) + 2 * sum(f(xj)) + f(b))
return(approx)}
Выше приведен код для генерации интеграции. И следующие коды - это некоторые функции, которые будут включены в функцию правдоподобия:
f_d1_55<-function(x) {return((0.2*exp(-(log(x)-mu)^2/(2*sig))/(sqrt(2*pi*sig)*x))*(exp(-lam*((55-x)^alph))))}
f_d1_56<-function(x) {return((0.2*exp(-(log(x)-mu)^2/(2*sig))/(sqrt(2*pi*sig)*x))*(exp(-lam*((56-x)^alph))))}
D1_55<-beta*composite.trapezoid(f_d1_55, 1e-200, 55, 120)
D1_56<-beta*composite.trapezoid(f_d1_56, 1e-200, 56, 120)
f_d2_55<-function(x) {return((0.2*exp(-(log(x)-mu)^2/(2*sig))/(sqrt(2*pi*sig)*x))*(exp(-lam*((56-x)^alph))))}
f_d2_56<-function(x) {return((0.2*exp(-(log(x)-mu)^2/(2*sig))/(sqrt(2*pi*sig)*x))*(exp(-lam*((57-x)^alph))))}
D2_55<-beta*((1-beta)*composite.trapezoid(f_d2_55, 1e-15, 55, 100)+composite.trapezoid(f_d2_55, 55, 56, 100))
D2_56<-beta*((1-beta)*composite.trapezoid(f_d2_56, 1e-15, 56, 100)+composite.trapezoid(f_d2_56, 56, 57, 100))
f_i11_55 <- function(x) {return((0.2*exp(-(log(x)-mu)^2/(2*sig))/(sqrt(2*pi*sig)*x))*(exp(-lam*((55-x)^alph))-exp(-lam*((56-x)^alph))))}
f_i11_56 <- function(x) {return((0.2*exp(-(log(x)-mu)^2/(2*sig))/(sqrt(2*pi*sig)*x))*(exp(-lam*((56-x)^alph))-exp(-lam*((57-x)^alph))))}
f_i12_55 <- function(x) {return((0.2*exp(-(log(x)-mu)^2/(2*sig))/(sqrt(2*pi*sig)*x))*(1-exp(-lam*((56-x)^alph))))}
f_i12_56 <- function(x) {return((0.2*exp(-(log(x)-mu)^2/(2*sig))/(sqrt(2*pi*sig)*x))*(1-exp(-lam*((57-x)^alph))))}
I1_55<-(1-beta)*composite.trapezoid(f_i11_55, 1e-200, 55, 20)+composite.trapezoid(f_i12_55, 55, 56, 20)
I1_56<-(1-beta)*composite.trapezoid(f_i11_56, 1e-200, 56, 20)+composite.trapezoid(f_i12_56, 56, 57, 20)
f_i21_55 <- function(x) {return((0.2*exp(-(log(x)-mu)^2/(2*sig))/(sqrt(2*pi*sig)*x))*(exp(-lam*((56-x)^alph))-exp(-lam*((57-x)^alph))))}
f_i21_56 <- function(x) {return((0.2*exp(-(log(x)-mu)^2/(2*sig))/(sqrt(2*pi*sig)*x))*(exp(-lam*((57-x)^alph))-exp(-lam*((58-x)^alph))))}
f_i22_55 <- function(x) {return((0.2*exp(-(log(x)-mu)^2/(2*sig))/(sqrt(2*pi*sig)*x))*(1-exp(-lam*((57-x)^alph))))}
f_i22_56 <- function(x) {return((0.2*exp(-(log(x)-mu)^2/(2*sig))/(sqrt(2*pi*sig)*x))*(1-exp(-lam*((58-x)^alph))))}
I2_55<-(1-beta)*(1-beta)*composite.trapezoid(f_i21_55, 1e-200, 55, 20)+(1-beta)*composite.trapezoid(f_i21_55, 55, 56, 20)+composite.trapezoid(f_i22_55, 56, 57, 20)
I2_56<-(1-beta)*(1-beta)*composite.trapezoid(f_i21_56, 1e-200, 56, 20)+(1-beta)*composite.trapezoid(f_i21_56, 56, 57, 20)+composite.trapezoid(f_i22_56, 57, 58, 20)
Следующие A11 и A12 являются конечными функциями, которые включают в себя вышеуказанные функции. И logL_xray - функция логарифмического правдоподобия.
A11<-3*log(D1_55)+2*log(D1_56)+8*log(D1_57)+1*log(D1_58)+0*log(I1_55)+1*log(I1_56)+1*log(I1_57)+1*log(I1_58)+(1098-3-0)*log(1-D1_55-I1_55)+(1087-2-1)*log(1-D1_56-I1_56)+(917-8-1)*log(1-D1_57-I1_57)+(796-1-1)*log(1-D1_58-I1_58)
A12<-0*log(D2_55)+2*log(D2_56)+1*log(D2_57)+0*log(D2_58)+4*log(I2_55)+1*log(I2_56)+2*log(I2_57)+2*log(I2_58)+(1006-0-4)*log(1-D2_55-I2_55)+(998-2-1)*log(1-D2_56-I2_56)+(849-1-2)*log(1-D2_57-I2_57)+(729-0-2)*log(1-D2_58-I2_58)
logL_xray<-function(t,mu,sig,alph,lam) sum(A11,A12)
optim(par=c(mu=4,sig=0.01,alph=0.2,lam=1.5),fn=logL_xray)
И получил вывод примерно так:
$par
mu sig alph lam
4.00 0.01 0.20 1.50
$value
[1] -155.0813
$counts
function gradient
5 NA
$convergence
[1] 0
, что означает, что функция «optim» не обновляется, а дает те же начальные значения. Как я могу решить эту проблему?