Expectation of number of right truncation; numerical integration

Expectation_numtrunc.R

#intensity function
lambda <- function(x, shape){
  shape*x^(shape-1)
}

#number of truncated observation
counttrunc <- function(tau, WS, weibull_param,
                       shape=1, maxit=10000){
  Et = 0
  count = 0L
  for (it in 1:maxit) {
    z = runif(1) #next infection
    Et =  (-log(z) + Et^shape)^(1/shape) #Weibull process
    if(Et>tau){break}
      u = runif(1)
      y = rweibull(1, weibull_param$shape, weibull_param$scale)
      RS = Et+y+WS*u
      count = count + as.integer(RS>tau)
  }
  return(count)
}

#Monte-Carlo integration
mcint = function(tau, WS, weibull_param, sigma=1, iter=10000){
  tab = table(ressim)
  intv = numeric(iter)
  for(i in 1:iter){
    u = runif(1)
    int1 = integrate(function(x){
      lambda(x,sigma)*pweibull(tau-WS*u-x, weibull_param$shape, weibull_param$scale, lower.tail = FALSE)
      }, 
      0, tau)
    intv[i] = int1$value  
  }
  mean(intv)
}

WS = 2 #width of interval censored
tau = 10 #width of observation
sigma = 0.9 #parameter of the infection process
weibull_param <- list(shape=1.5, scale=2)
ressim = replicate(10000, counttrunc(tau, WS, weibull_param, sigma))
mu = mcint(tau, WS, weibull_param, sigma)
print(mean(ressim))
print(mu)

tab = table(ressim)
prob = tab/sum(tab)
xran = range(as.integer(names(tab)))
plot(prob)
points(xran[1]:xran[2],
       dpois(xran[1]:xran[2],mu), type = "b")