R中的似然最大化

Question

在 R 中，我编写了一个包含两个递归计算的对数似然函数。对数似然函数正常工作（它给出了已知参数值的答案），但是当我尝试使用 optim() 最大化它时，它需要太多时间。如何优化代码？提前感谢您的想法。

这是具有使用 copula 函数的依赖结构的马尔可夫状态切换模型的对数似然函数。

在 for 循环中命名为 g：

在 for 循环中命名为 p：

在代码中命名为f：

一些数据：

u <- cbind(rt(100,10),rt(100,13))

f函数：

f=function(u,p,e1,e2){
     s=diag(2);s[1,2]=p
     ff=dcopula.gauss(cbind(pt(u[,1],e1),pt(u[,2],e2)),Sigma=s)*dt(u[,1],e1)*dt(u[,2],e2)
     return(ff)
  }

对数似然函数：

loglik=function(x){
  p11<-x[1];p12<-x[2];p21<-x[3];p22<-x[4];p31<-x[5];p32<-x[6];r<-x[7];a1<-x[8];a2<-x[9];s<-x[10];b1<-x[11];b2<-x[12];t<-x[13];c1<-x[14];c2<-x[15]
  p1=c(numeric(nrow(u)));p2=c(numeric(nrow(u)));p3=c(numeric(nrow(u)))
  g=c(numeric(nrow(u)))
  p1_0=.3
  p2_0=.3
  g[1]<-(p1_0*f(u,r,a1,a2)[1])+(p2_0*f(u,s,b1,b2)[1])+((1-(p1_0+p2_0))*f(u,t,c1,c2)[1])
  p1[1]<-((p1_0*p11*f(u,r,a1,a2)[1])+(p2_0*p21*f(u,r,a1,a2)[1])+((1-(p1_0+p2_0))*p31*f(u,r,a1,a2)[1]))/g[1]
  p2[1]<-((p1_0*p12*f(u,s,b1,b2)[1])+(p2_0*p22*f(u,s,b1,b2)[1])+((1-(p1_0+p2_0))*p32*f(u,s,b1,b2)[1]))/g[1]
  p3[1]<-((p1_0*(1-(p11+p12))*f(u,t,c1,c2)[1])+(p2_0*(1-(p21+p22))*f(u,t,c1,c2)[1])+((1-(p1_0+p2_0))*(1-(p31+p32))*f(u,t,c1,c2)[1]))/g[1]
  for(i in 2:nrow(u)){
    g[i]<-(p1[i-1]*p11*f(u,r,a1,a2)[i])+(p1[i-1]*p12*f(u,s,b1,b2)[i])+(p1[i-1]*(1-(p11+p12))*f(u,t,c1,c2)[i])+
      (p2[i-1]*p21*f(u,r,a1,a2)[i])+(p2[i-1]*p22*f(u,s,b1,b2)[i])+(p2[i-1]*(1-(p21+p22))*f(u,t,c1,c2)[i])+
      (p3[i-1]*p31*f(u,r,a1,a2)[i])+(p3[i-1]*p32*f(u,s,b1,b2)[i])+(p3[i-1]*(1-(p31+p32))*f(u,t,c1,c2)[i])
    p1[i]<-((p1[i-1]*p11*f(u,r,a1,a2)[i])+(p1[i-1]*p12*f(u,s,b1,b2)[i])+(p1[i-1]*(1-(p11+p12))*f(u,t,c1,c2)[i]))/g[i]
    p2[i]<-((p2[i-1]*p21*f(u,r,a1,a2)[i])+(p2[i-1]*p22*f(u,s,b1,b2)[i])+(p2[i-1]*(1-(p21+p22))*f(u,t,c1,c2)[i]))/g[i]
    p3[i]<-((p3[i-1]*p31*f(u,r,a1,a2)[i])+(p3[i-1]*p32*f(u,s,b1,b2)[i])+(p3[i-1]*(1-(p31+p32))*f(u,t,c1,c2)[i]))/g[i]
  }
  return(-sum(log(g)))
}

优化：

library(QRM)
library(copula)
start=list(0,1,0,0,0,0,1,9,7,-1,10,13,1,6,4)
##
optim(start,loglik,lower=c(rep(0,6),-1,1,1,-1,1,1,-1,1,1),
  upper=c(rep(1,6),1,Inf,Inf,1,Inf,Inf,1,Inf,Inf),
  method="L-BFGS-B") -> fit

Answer 1

这看起来像是 Stack-Overflow 的问题。

我想到的是：

1. 定义一个包含值f(.,.,.,.) 的向量，以避免对同一函数进行k*nrow(u) 评估，只需调用这些感兴趣的条目。

2. 似乎循环可以被矩阵和/或向量乘积代替。但是，如果没有进一步的信息，就不清楚代码在做什么，而且需要很长时间才能从代码中提取这些信息。