如何组织最终表以模拟不同样本量的r中的最大似然？

问题描述

我使用R程序进行模拟，但是无法在R中组织/制作决赛桌。

rep=100
n=20

# define starting values 
mu<-100
sig<-3
theta=c(mu,sig) # store starting values

#Tables
#********

LHE=array(0,c(2,rep));
rownames(LHE)= c("MLE_mu","MLE_sigma")

bias= array(0,rep));
rownames(bias)= c("bias_mu","bias_sigma")

#Simulation {FOR LOOP}
#***********************
set.seed(1)

for(i in 1:rep){
  myx <- rnorm(100,100,3)  

  loglikenorm<-function(x,myx) # always use x to hold parameter values
  {
    mu<-x[1]
    sig<-x[2]
    n<-length(myx)
    loglike<- -n*log(sig*sqrt(2*pi))- sum((1/(2*sig^2))*(myx-mu)^2) # note
    # use of sum
    loglike<- -loglike
  }
  
  result<-nlm(loglikenorm,theta,myx=myx,hessian=TRUE,print.level=1) #ML estimation using nlm
  mle<-result$estimate #extract and store mles
  LHE[,i]=  c(mle[1],mle[2])
  bias[,i]= c(mle[1]-theta[1],mle[2]-theta[2])
} # end for i

L <-round(apply(LHE,1,mean),3) # MLE  of all the applied iterations  

bs <-round(apply(bias,3) # bias of all the applied iterations 
row<- c(L,bs); row

这将运行MLE 100次。我想针对不同的样本量（n = 20,50,100）和不同的参数集c(c(mu= 100,sigma=3),c(mu=80,sigma=4))

计算MLE

我想要两件事，第一件事是如何运行代码以针对不同的样本量和不同的参数集计算MLE。第二个如何使用R程序组织输出（如附件的图像）。

The final table from R

任何帮助将不胜感激。

解决方法

我建议您构建一个用于MLE估计的函数，然后将其应用于不同的参数设置。这里使用您的代码的功能：

#Function
mymle <- function(n,mu,sig,rep)
{
  #Set reps
  rep=rep
  # define starting values 
  mu<-mu
  sig<-sig
  theta=c(mu,sig) # store starting values
  #Tables
  LHE=array(0,c(2,rep));
  rownames(LHE)= c("MLE_mu","MLE_sigma")
  #Bias
  bias= array(0,rep));
  rownames(bias)= c("bias_mu","bias_sigma")
  #Simulation
  set.seed(1)
  #Loop
  for(i in 1:rep){
    myx <- rnorm(rep,sig)  
    
    loglikenorm<-function(x,myx) # always use x to hold parameter values
    {
      mu<-x[1]
      sig<-x[2]
      n<-length(myx)
      loglike<- -n*log(sig*sqrt(2*pi))- sum((1/(2*sig^2))*(myx-mu)^2) # note
      # use of sum
      loglike<- -loglike
    }
    
    result<-nlm(loglikenorm,theta,myx=myx,hessian=TRUE,print.level=1) #ML estimation using nlm
    mle<-result$estimate #extract and store mles
    LHE[,i]=  c(mle[1],mle[2])
    bias[,i]= c(mle[1]-theta[1],mle[2]-theta[2])
  } # end for i
  #Format results
  L <-round(apply(LHE,1,mean),3) # MLE  of all the applied iterations  
  bs <-round(apply(bias,3) # bias of all the applied iterations 
  row<- c(L,bs)
  #Format a label
  lab <- paste0('n= ',n,';',' mu= ',' sig= ',sig)
  row2 <- c(lab,row)
  row2 <- as.data.frame(t(row2))
  return(row2)
}

现在我们申请不同的参数设置：

#Example 1
ex1 <- mymle(n = 20,mu = 100,sig = 3,rep = 100)
ex2 <- mymle(n = 50,rep = 100)
ex3 <- mymle(n = 100,rep = 100)
#Example 2
ex4 <- mymle(n = 20,mu = 80,sig = 4,rep = 100)
ex5 <- mymle(n = 50,rep = 100)
ex6 <- mymle(n = 100,rep = 100)

最后，我们将所有结果绑定为接近您想要的输出：

#Bind all
df <- rbind(ex1,ex2,ex3,ex4,ex5,ex6)

输出：

                       V1 MLE_mu MLE_sigma bias_mu bias_sigma
1  n= 20; mu= 100; sig= 3  99.98     3.015   -0.02      0.015
2  n= 50; mu= 100; sig= 3  99.98     3.015   -0.02      0.015
3 n= 100; mu= 100; sig= 3  99.98     3.015   -0.02      0.015
4   n= 20; mu= 80; sig= 4 79.974      4.02  -0.026       0.02
5   n= 50; mu= 80; sig= 4 79.974      4.02  -0.026       0.02
6  n= 100; mu= 80; sig= 4 79.974      4.02  -0.026       0.02

我分别测试了循环中参数的不同组合，结果是相等的。就像基于大数原理的提醒一样，当我们添加更多rep时，结果将趋于真实值，因此，如果将rep更改为1000，则值将改变。

mle r r