0 
library(mvtnorm)

set.seed(14)
n=10000 
sigmatrue<-1
 rhotrue<-0.3
  b1=0.05
  b0=0
   y<-arima.sim(model=list(ar=c(0.3)),n=10000 ,sd=sigmatrue)#kataskevi 
    #xronoseiras

 x=rep(0,n)

for(i in 1:n){
  x[i]=i
}
for(t in 1:n)
{   
  y[t]=y[t]+b0+b1*x[t]
}
est=arima(y,order=c(1,0,0),xreg=x,include.mean=TRUE,method="ML",kappa=1e+06)
cens<-rep(0, n)
c=(9/10)*(n*b1+b0)
for (i in 1:n) {
  if(y[i]>c){
    y[i]<-c
    cens[i]<-1
  }
}

  ll<-function(p){
    sigma=matrix(c(p[2]^2/(1-p[3]^2), p[2]^2*p[3]/(1-p[3]^2),p[2]^2*p[3]/(1-p[3]^2),p[2]^2/(1-p[3]^2)),ncol=2,nrow=2,byrow=TRUE)
    likelihood<-rep(0,n)
    for(t in 2 :n){
      if(cens[t]==0 & cens[t-1]==0){
        likelihood[t]<-dnorm(((y[t]-(p[1]+p[4]*t)-p[3]*(y[t-1]-(p[1]+p[4]*(t-1)))/p[2]) )/p[2])
      }
      else if(cens[t]==0 & cens[t-1]==1){
        likelihood[t]<-(1/(1-pnorm((c-(p[1]+p[4]*t)*sqrt(1-p[3]^2)/p[2]))*sqrt(1-p[3]^2)/p[2]*dnorm(((y[t]-(p[1]+p[4]*t)*sqrt(1-p[3]^2))/p[2])*(1-pnorm(((c-(p[1]+p[4]*(t))-p[3]*(y[t]-(p[1]+p[4]*(t-1)))/p[2])))))))
      }
      else if(cens[t]==1 & cens[t-1]==0){
        likelihood[t]<-1-pnorm(((c-(p[1]+p[4]*t)-p[3]*(y[t-1]-(p[1]+p[4]*(t-1)))/p[2])))
      }
      else
      {
        likelihood[t]<-(((pmvnorm(lower=c, upper=Inf , mean=c(p[1]+p[4]*(t-1),p[1]+p[4]*t),sigma=sigma))/(1-pnorm((c-(p[1]+p[4]*(t-1))*sqrt(1-p[3]^2)/p[2])))))
      }
    }


    f0=(sqrt(1-p[3])/p[2]*dnorm(((y[1]-p[1]-p[4])*sqrt(1-p[3]^2))/p[2]))

    likelihood[1]=f0
    #Ta prosthesa
    if (any(likelihood==0)){
      likelihood[likelihood==0] = 0.000001 #poly mikros arithmos
    }

    if (any(likelihood==Inf)){
      likelihood[likelihood==Inf] = 1 #poly megalos h 1, an milame gia pi8anothta
    }

    if (any(is.nan(likelihood))){
      likelihood[is.nan(likelihood)] = 0.000001
    }

    minusloglike=-sum(log(likelihood))
    #l1=list(Minusloglike=minusloglike,Loglikelihood=log(likelihood))
    return(minusloglike)
  }
  fit<-optim(c(0,1,0.3,0.05),ll,method="L-BFGS-B",lower=c(-Inf,0.001,-0.999,-Inf),upper = c(Inf,Inf,0.999,Inf),hessian=TRUE)
  fisher.info<-solve(fit$hessian)
  fisher.info
  prop.sigma<-sqrt(diag(fisher.info))
  sigmas<-diag(prop.sigma)
  upper<-fit$par+1.96*sigmas
  lower<-fit$par-1.96*sigmas
  interval<-data.frame(value=fit$par, lower=diag(lower),upper=diag(upper))
  interval

我运行此代码（它用于带有协变量的审查一阶自回归过程，我有 4 个 x(t) 的案例，x(t-1) 要么是审查的，要么是未审查的，我不希望可能性接近零并且inf).我得到错误

if (any(likelihood == Inf)) { : 需要 TRUE/FALSE 的缺失值调用自：fn(par, ...)

该程序适用于 n=100 但当 n 大于 100 时出现此错误。我认为这个错误会导致对四个参数（b1、rho、sigma、b0）的错误估计。有人知道我能做什么吗？

谢谢您的帮助。

4

0 回答 0