如何使用R中的蒙特卡罗模拟来改善威布尔分布的最大似然估计？

如何解决如何使用R中的蒙特卡罗模拟来改善威布尔分布的最大似然估计？

我正在使用两个参数的Weibull分布，生成样本并估计对数似然。 我正在使用的威布尔密度函数是这样的：

我正在使用的Weibull日志可能性是这样的：

我正在使用的Weibull渐变是这样的：

我使用的优化方法是optim中的L-BFGS-B方法。当我在函数optim中使用数值梯度时，与使用解析梯度时相比，结果更好，执行时间也更短，而且这种情况应该不会发生。所以现在我在这里寻求帮助。 ：）

Weibull.mc.ga = function(Nrep=10000,Nobs=100,semente=2000,scale=2,shape=2)
{
  # loglikelihood function
  logLikWeibull = function(theta){
    alpha = theta[1]; beta = theta[2]
    loglik = Nobs*log(beta)-Nobs*beta*log(alpha) + (beta-1)*(sum(log(y))) - (1/(alpha^beta))*(sum(y^beta))
    
    return(loglik)
  }
  # score function (gradient)
  scoreFn = function(theta){
    alpha = theta[1]; beta = theta[2]
    cbind(beta*((alpha^(-beta-1))*sum(y^beta)) - Nobs*(beta/alpha),(Nobs/beta)-Nobs*(log(alpha)) + sum(log(y)) - (log(beta)/((alpha)^beta))*sum(y^beta))
  }
  # begin time count
  tempo.inicio = Sys.time()
  # estimative vectors
  emvalpha = rep(0,Nrep)
  emvbeta = rep(0,Nrep)
  set.seed(semente) # seed
  contadorFalhas = 0 # count failures
  # Monte Carlo Simulation
  i = 1
  while(i <= Nrep){
    # sample
    y <<- rweibull(Nobs,shape=shape,scale=scale)
    alphachute = 1;   betachute = 1 # initial values for parameters
    chute = c(alphachute,betachute)
    # loglikelihood maximization
    ir = optim(par = chute,fn = logLikWeibull,gr = scoreFn,method="L-BFGS-B",lower=c(0.01,0.01),control=list(fnscale=-1,trace = 3))
    alpha = ir$par[1]
    beta =ir$par[2]

    # convergence check
    if(ir$convergence == 0){
      emvalpha[i] = ir$par[1]
      emvbeta[i] = ir$par[2]
      mediaalpha = mean(emvalpha)
      mediabeta = mean(emvbeta)
      varianciaalpha = var(emvalpha)
      varianciabeta = var(emvbeta)
      i = i + 1
    }
    else{
      contadorFalhas = contadorFalhas + 1 #count failures
    }
  } # End Monte Carlo Simulation
  # average estimates,biases and relative biases
  alphamedio = mean(emvalpha) # mean
  betamedio = mean(emvbeta) # mean
  alphavies = alphamedio - scale # bias
  betavies = betamedio - shape # bias
  alphaviesrel = alphavies/scale # relative bias
  betaviesrel = betavies/shape # relative bias
  # mean square errors (EQM means MSE)
  alphaeqm = alphavies^2 + var(emvalpha)
  betaeqm = betavies^2 + var(emvbeta)

  # lower limit of confidence interval for alpha
  Linfalpha = mediaalpha - (1.96)*varianciaalpha
  # upper limit of confidence interval for alpha
  Lsupalpha = mediaalpha + (1.96)*varianciaalpha
  # lower limit of confidence interval for beta
  Linfbeta = mediabeta - (1.96)*varianciabeta
  # upper limit of confidence interval for beta
  Lsupbeta = mediabeta + (1.96)*varianciabeta
  mIntervalodeConf = matrix(c(Linfalpha,Lsupalpha,Linfbeta,Lsupbeta),2,byrow = FALSE)
  rownames(mIntervalodeConf) = c("inferior limit","superior limit")
  colnames(mIntervalodeConf) = c("Interval alpha","Interval beta")
  # matrix of results
  mResultados = matrix(c(alphamedio,betamedio,alphavies,betavies,alphaviesrel,betaviesrel,alphaeqm,betaeqm),4,byrow=TRUE)
  rownames(mResultados) = c("mean","bias","rel. bias","MSE")
  colnames(mResultados) = c("emv alpha","emv beta")
  # calculating execution time
  tempo.fim = Sys.time()
  tempo.exec = tempo.fim - tempo.inicio
  # list of all results
  finalresults = list(Nobs=Nobs,Nrep=Nrep,seed=semente,alpha=scale,beta=shape,failures=contadorFalhas,results=mResultados,interval=mIntervalodeConf,horario = tempo.inicio,timeexec = tempo.exec)
  return(finalresults)
}