set.seed(5)
beta0<--8
beta1<-0.03
gamma<-0.0105
#n is sample size
n<-50
N<-2000
beta0hat1<-NULL
beta1hat1<-NULL
gammahat1<-NULL
cp<-NULL
for (i in 1:N)
{
  u<-runif(n)
  x<-rnorm(n)
  c<-rexp(n,1/725)
  t1<-(1/gamma)*log(1-((gamma/exp(beta0+beta1*x))*log(1-u)))
  t<-pmin(t1,c)
  delta<-1*(t1>c)
  cp[i]<-length(delta[delta==1])/n
  delta[delta==1]<-ifelse(rbinom(length(delta[delta==1]),1,0.75),1,2)
  dat=data.frame(t,delta,x)
  dat$interval[delta==2] <- as.character(cut(dat$t[delta==2], breaks=seq(0, 600, 50)))
  labs <- cut(dat$t[delta==2], breaks=seq(0, 600, 50))
  dat$lower[delta==2]<-as.numeric( sub("\((.+),.*", "\1", labs) )
  dat$upper[delta==2]<-as.numeric( sub("[^,]*,([^]]*)\]", "\1", labs) )
  data0<-dat[which(dat$delta==0),]#uncensored data
  data1<-dat[which(dat$delta==1),]#right censored data
  data2<-dat[which(dat$delta==2),]#interval censored data
  library(maxLik)
  ll1<-function(para)
  {
    b0<-para[1]
    b1<-para[2]
    g<-para[3]
    e<-sum((b0+b1*data0$x)+g*data0$t+(1/g)*exp(b0+b1*data0$x)*(1-exp(g*data0$t)))
    r<-sum((1/g)*exp(b0+b1*data1$x)*(1-exp(g*data1$t)))
    i<-sum(log(exp((1/g)*exp(b0+b1*data2$x)*(1-exp(g*data2$lower)))-exp((1/g)*exp(b0+b1*data2$x)*(1-exp(g*data2$upper)))))
    l<-e+r+i
    return(l)
  }
  mle1<-maxLik(logLik=ll1,start=c(para<-c(-8,0.03,0.0105)))
  beta0hat1[i]<-mle1$estimate[1]
  beta1hat1[i]<-mle1$estimate[2]
  gammahat1[i]<-mle1$estimate[3]
}

biasbeta0hat1<-abs(mean(beta0hat1)-beta0)
biasbeta0hat1
biasbeta1hat1<-abs(mean(beta1hat1)-beta1)
biasbeta1hat1
biasgammahat1<-abs(mean(gammahat1)-gamma)
biasgammahat1