区間打ち切りデータの最尤推定のシミュレーション：JuliaとR

Julia です。

@time using Optim
@time using Distributions
@time using StatsFuns

function make_dat_gamma0(n,a,b)
t = rand(Gamma(a,b),n)
E_len = rand(Exponential(),n)
u = rand(n)
E_R = t - E_len .* (1.0.-u)
E_L = t + E_len .* u
E_L, E_R
end

function make_dat_gamma(n,a,b)
    E_R = zeros(n)
    E_L = zeros(n)
    for i in 1:n
    t = rand(Gamma(a,b))
    E = rand(Exponential())
    u = rand()
    E_R[i] = max(0.0,t - E * (1.0-u))
    E_L[i] = t + E * u
    end
    E_L, E_R
end


function llgamma_intcens0(param, L, R)
    a = exp(param[1])
    b = exp(param[2])
    ll = sum(logsubexp.(logcdf.(Gamma(a,b),L),logcdf.(Gamma(a,b),R)))
    return -ll
end


function llgamma_intcens(param, L, R)
    a = exp(param[1])
    b = exp(param[2])
    n = size(R,1)
    ll = 0
    for i in 1:n
        ll += logsubexp(logcdf(Gamma(a,b),L[i]),logcdf(Gamma(a,b),R[i]))
    end
    return -ll
end


@time make_dat_gamma0(100,2.0,1.0)
@time make_dat_gamma(100,2.0,1.0)
@time llgamma_intcens0([0.0,0.0],E_L,E_R)
@time llgamma_intcens([0.0,0.0],E_L,E_R)

# ドット . を使って書くより for で回したほうが速いみたい.


function sim_gamma(n,a,b,iter)
    out = zeros(iter,2)
    for i in 1:iter
    E_L, E_R = make_dat_gamma(n,a,b)
    opt = optimize(vars -> llgamma_intcens(vars, E_L, E_R), zeros(2),BFGS(),Optim.Options(f_reltol=1.0e-8))
    parameters = exp.(Optim.minimizer(opt))
    out[i,:] = parameters
    end
    return out
end


@time simout = sim_gamma(100,2.0,0.2,1000)
#  5.546045 seconds (341.61 k allocations: 17.523 MiB)

R です。

simGam <- function(N,a,b,lambda=1){
  inc <- rgamma(N,a,b)
  u <- rexp(N,lambda) #exposed
  r_e <- runif(N)
  E_L <- inc + r_e*u
  E_R <- inc - (1-r_e)*u
  data.frame(E_L=E_L,E_R=pmax(0,E_R))
}

ll <- function(par,E_L,E_R){
  a <- exp(par[1])
  b <- exp(par[2])
  -sum(log(pgamma(E_L,a,b)-pgamma(E_R,a,b)))
}

simfunc <- function(n,a,b,iter){
  out <- matrix(0,iter,2)
  for(i in 1:iter){
    dat <- simGam(n,a,b)
    opt <- optim(c(0,0),ll,E_L=dat$E_L,E_R=dat$E_R,method = "BFGS",
                 control=list(reltol=1e-8))
    out[i,] <- exp(opt$par) 
  }
  return(out)
}

system.time({
out <- simfunc(n=100, a=2, b=5.0,iter=1000)
})
#   user  system elapsed 
#  4.568   0.024   4.674