下面是程序：
library(emplik)
library(stats)
require(MASS)

p=5
t=10
L=1
lam=5
nt=rep(0,L)
ww=rep(0,L)
ind=0
Z=rep(0,1000)
q1=rep(0,L)
q2=rep(0,L)

for(l in 1:L){
  nt[l]=rpois(1,lam*t)
  x1 <- rexp(nt[l]*p)
  x <- matrix(x1,nt[l],p)
  mu=rep(1,p)
  theta=lam*mu

  g1<-function(z,theta,agama){
    temtheta<-matrix(theta*t/nt[l],nt[l],p)
    deno<-diag(c(1/(1+(z-temtheta)%*%agama)))
    nomi<-z-temtheta
    temp1<-deno%*%nomi
    temp<-colSums(temp1)/nt[l]
    temp
  }

  g2<-function(z,thate,agama){
    temtheta<-matrix(theta*t/nt[l],nt[l],p)
    deno2<-diag(c(-1/(1+(z-temtheta)%*%agama)))
    nomi2<-matrix(1,nt[l],p)
    temp2<-deno2%*%nomi2%*%agama
    tempp<-colSums(temp2)*t/nt[l]
  }

  cal<-function(z,theta,agama){
    sol1<-(g1(z,theta,agama))^2
    sol2<-(g2(z,theta,agama))^2
    return(sum(sol1)+sum(sol2))
  }

  Fin<-function(z,theta,agama){
    temtheta<-matrix(theta*t/nt[l],nt[l],p)
    temp<-1+(z-temtheta)%*%agama
    fi<-sum(log(temp))
    fi
  }

  "%^%" <- function(x, n)
    with(eigen(x), vectors %*% (values^n * t(vectors)))

  agama<-nlm(cal,z=x,theta=theta,rep(0.01,p),iterlim=500)$estimate

  ww[l]<-2*Fin(x,theta,agama)

  for(i in 1:1000){
    ome=0
    sig=0
    for(j in 1:nt[l]){
      ome=ome+x[j,]%*%t(x[j,])
    }
    omega=ome/nt

    for(j in 1:nt[l]){
      sig=sig+(x[j,]-colMeans(x))%*%t(x[j,]-colMeans(x))
    }
    sigma=sig/nt
    sigma <- sigma%^% (-0.5)
    ma=sigma%*%omega%*%(sigma)
    z<-mvrnorm(1,rep(0,p),Sigma<-ma)
    Z[i]<-t(z)%*%(z)
  }



  q1[l]<-quantile(Z,probs=0.05)
  q2[l]<-quantile(Z,probs=0.95)

  if(ww[l]<q2[l] & ww[l]>q1[l]){
    ind=ind+1
  }
  ind
}

covp=ind/L


程序内容是检验一个复合泊松过程的分布，每一个过程都是一个向量，等于有nt个x，每个x都是p维的向量。
L就是循环次数，当L取1的时候程序ok可以运行，超过1报错了。报错内容如下：
错误于eigen(x) : infinite or missing values in 'x'
此外: 警告信息：
1: In ome/nt : 长的对象长度不是短的对象长度的整倍数
2: In sig/nt : 长的对象长度不是短的对象长度的整倍数

十万火急！万分感谢！！！