0

我在使用似然函数中的积分函数时遇到了困难。似然函数涉及四个函数，F1R、F1L、F2R、F2L，它们由积分定义。

F1R= 积分(0,R) exp(-(alpha1 v)^tau1 - (alpha2 v)^tau2) ((tau1 (alpha1^tau1)*v^(tau1-1)))

F1L= 积分(0,L) exp(-(alpha1 v)^tau1 - (alpha2 v)^tau2) ((tau1 (alpha1^tau1)*v^(tau1-1)))

F2R= 积分(0,t) exp(-(alpha1 v)^tau1 - (alpha2 v)^tau2) ((tau2 (alpha2^tau2)*v^(tau2-1)))

F2R= 积分(0,t) exp(-(alpha1 v)^tau1 - (alpha2 v)^tau2) ((tau2 (alpha2^tau2)*v^(tau2-1)))

其中 R 和 L 在数据库中给出。

    log_veros <-function(param,x){
  beta01 <- exp(param[1])
  beta11 <- exp(param[2])
  beta12 <- exp(param[3])
  beta13 <- exp(param[4])
  beta02 <- exp(param[5])
  beta21 <- exp(param[6])
  beta22 <- exp(param[7])
  beta23 <- exp(param[8])
  beta1 =c(beta11,beta12,beta13)
  beta2 =c(beta21,beta22,beta23)
  X<- cbind(x$x1,x$x2,x$x3)
  eta1 <- exp(beta01 + X%*%beta1)
  eta2 <- exp(beta02 + X%*%beta2)
  eta <- eta1 + eta2
  Delta0 <- x$Delta0
  Delta1 <- x$Delta1
  Delta2 <- x$Delta2
  F1L <- rep(NA,n)
  F1R <- rep(NA,n)
  F2L <- rep(NA,n)
  F2R <- rep(NA,n)
  for(i in 1: n){
    myf1 <- function(v){
      eta1[i]*exp(-(eta[i])*v)
    } 
    myf2 <- function(v){
      eta2[i]*exp(-(eta[i])*v)
    }
    F1L[i] <-integrate(myf1, lower=0,upper = dados$L[i])$value
    F1R[i] <-integrate(myf1, lower=0,upper = dados$R[i])$value
    F2L[i] <-integrate(myf2, lower=0,upper = dados$L[i])$value
    F2R[i] <-integrate(myf2, lower=0,upper = dados$R[i])$value
   }
  aux<- log((F1R -F1L)^Delta1*(F2R-F2L)^Delta2*(1-(F1L+F2L))^Delta0)
  return(-sum(aux)) 
}

关注数据库的一部分以便更好地理解：

 data
             L          R status Delta0 Delta1 Delta2 x1 x2          x3
1  0.031678365 0.41981090      1      0      1      0  1  1  0.49043448
2  0.015044380 0.33537736      1      0      1      0  0  1  0.07733734
3  0.003636517 0.20048560      1      0      1      0  1  1 -1.78858071
4  0.014671982 0.40071384      1      0      1      0  0  1 -0.76329354
5  0.076556199        Inf      0      1      0      0  0  1 -1.43815999
6  0.000000000 0.03239601      2      0      0      1  0  1  0.84835412
7  0.046961166        Inf      0      1      0      0  1  1 -1.31396966

我认为我integrate错误地使用了该功能，因为它在optimx

library(optimx)
p= c(log(5),log(1.5),log(4),log(1.5),log(4),log(2),log(1.5),log(2.5))

emv <- optimx(par = p, fn =log_veros, x=data,  
              method = c("nlm", "BFGS", "Rcgmin", "nlminb"))
4

0 回答 0