Оценка максимального правдоподобия, основанная на Ньютоне-Рафсоне и методе моментов

Question

Я проводил небольшое симуляционное исследование, чтобы асимптотически исследовать свойства метода моментов и оценок максимального правдоподобия.

Метод оценки моментов легко получить (он задан второй строкой моего кода), но для mle мне пришлось написать алгоритм Newton-Raphson algoritmh (который отлично работает для одного образца). Mle должен использовать метод оценки моментов в качестве начальной точки (a0), поскольку таким образом он обладает некоторыми оптимальными статистическими свойствами. Здесь

x<-rbeta(500,0.5,3) 
mom<-3*mean(x)/(1-mean(x)) 

mlea<-function(x,a0,numstp=100,eps=0.001){ 
    n=length(x) 
    numfin=numstp 
    ic=0 
    istop=0 
    while(istop==0){ 
     ic=ic+1 
     lprime=n/a0+n/(a0+1)+n/(a0+2)+sum(log(x)) 
     ldprime=-n/a0^2-n/(a0+1)^2-n/(a0+2)^2 
     a1=a0-(lprime/ldprime) 
     check=abs((a1-a0)/a0) 
     if(check<eps){istop=1} 
     a0=a1 
    } 
    list(a1=a1,check=check,realnumstps=ic) 
} 

Это работает для одного образца, но могу ли я получить эти оценки для 1000 образцов? Как я мог легко обобщить этот процесс? Моя основная трудность связана с тем, что mle нуждается в маме в качестве отправной точки, и оба они должны быть вычислены из одного и того же образца.

Заранее спасибо.

Answer 1

Я думаю, вы хотите это сделать?

n<-100 
replicate(n, { 
    x<-rbeta(500,0.5,3) 
    mom<-3*mean(x)/(1-mean(x)) 
    mlea(x, mom) 
}) 

Теперь вектор чисел не будет возвращен, потому что ваша функция mlea возвращает список. Скажем, значение из этого списка вы действительно заботитесь о a1, то вы могли бы сделать

n<-100 
replicate(n, { 
    x<-rbeta(500,0.5,3) 
    mom<-3*mean(x)/(1-mean(x)) 
    mlea(x, mom)$a1 
}) 

уведомление я называю $ a1 в конце вызова функции.

Итак, вот что происходит здесь, это репликация вытащит 500 новых наблюдений из вашего бета-распределения для каждой итерации в реплике (которая будет повторяться n раз), а затем вычислить маму на основе этого x, а затем дать результат mlea

Мои результаты?

replicate(3, { 
      x<-rbeta(500,0.5,3) 
      mom<-3*mean(x)/(1-mean(x)) 
      list(a1=mlea(x, mom)$a1, mom=mom) 
     }) 
# [,1]  [,2]  [,3]  
#a1 0.494497 0.522325 0.5153832 
#mom 0.4955767 0.5083678 0.5206997 

где каждый столбец здесь находится смотровая