如何组织最终的 table 来模拟具有不同样本量的 r 中的最大似然?

how to organize the final table to simulate maximum likelihood in r with different sample sizes?

我使用 R 程序进行了模拟,但我无法 organize/make R 中的最终 table。

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, c(2, 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,1, mean),3) # bias of all the applied iterations 
row<- c(L, bs); row

这将 运行 MLE 100 次。我想计算不同样本大小 (n=20,50,100) 和不同参数集的 MLE c(c(mu= 100, sigma=3), c(mu=80 , sigma=4))

我想要两件事,第一件事是如何 运行 代码来计算不同样本大小和不同参数集的 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, c(2, rep));
  rownames(bias)= c("bias_mu", "bias_sigma")
  #Simulation
  set.seed(1)
  #Loop
  for(i in 1:rep){
    myx <- rnorm(rep, mu, 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,1, mean),3) # bias of all the applied iterations 
  row<- c(L, bs)
  #Format a label
  lab <- paste0('n= ',n,';',' mu= ',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,mu = 100,sig = 3,rep = 100)
ex3 <- mymle(n = 100,mu = 100,sig = 3,rep = 100)
#Example 2
ex4 <- mymle(n = 20,mu = 80,sig = 4,rep = 100)
ex5 <- mymle(n = 50,mu = 80,sig = 4,rep = 100)
ex6 <- mymle(n = 100,mu = 80,sig = 4,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,数值也会发生变化。