R kode for moment- og ML-estimatene i gammafordelingen

 

 

# Til illustrasjon bruker vi dataene fra eksempel 4.27 og 7.14 om levetidene (i uker) for mus som er eksponert for 240 rad gamma str?ling.

# Vi leser f?rst inn dataene (jf. eksempel 7.14):

uker=c(152, 115, 109, 94, 88, 137, 152, 77, 160, 165, 125, 40, 128, 123, 136, 101, 62, 153, 83, 69)

 

 

# Vi finner momentestimatene (jf? eksempel 7.14):

n=length(uker)

beta=(sum(uker^2)/n-mean(uker)^2)/mean(uker)

alfa=mean(uker)/beta

 

 

# Vi lager en funksjon som finner log-likelihooden ganget med minus ?n:

negloglik=function(tetta,x)

{

? a=tetta[1]

? b=tetta[2]

? nn=length(x)

? nn*a*log(b)+nn*lgamma(a)-(a-1)*sum(log(x))+sum(x)/b

}

 

 

#Vi minimerer minus log-likelihooden numerisk (som er det same som ? maksimere log-likelihooden).

# Vi bruker da momentestimatene som startverdier for den numeriske optimeringen:

mle=optim(c(alfa,beta),negloglik,x=uker)$par

 

 

# Vi vil finne bootstrap estimatene for standardfeilene til momentestimatene og ML-estimatene.

 

# Vi bruker ML-estimatene som parameterverdier i bootstrapmetoden.

a.mle=mle[1]

b.mle=mle[2]

 

# Vi genererer B = 1000 bootstrap utvalg av st?rrelse n = 20 fra gammafordelingen med formparameter a.mle og skalaparameter b.mle og bestemmer momentestimatene og ML-estimatene for hvert av utvalgene:

B=1000

a.boot.moment=b.boot.moment=a.boot.mle=b.boot.mle=rep(0,B)

for (i in 1:B)

{

     uker.boot=rgamma(n,shape=a.mle,scale=b.mle)

???? b.boot.moment[i]=(sum(uker.boot^2)/n-mean(uker.boot)^2)/mean(uker.boot)

???? a.boot.moment[i]=mean(uker.boot)/b.boot.moment[i]

???? mle.boot=optim(c(a.boot.moment[i],b.boot.moment[i]),negloglik,x=uker.boot)$par

???? a.boot.mle[i]=mle.boot[1]

???? b.boot.mle[i]=mle.boot[2]

}

 

# Vi beregner bootstrap estimatene for standardfeilene

sd(a.boot.moment)

sd(b.boot.moment)

sd(a.boot.mle)

sd(b.boot.mle)