案例R语言F79MA 2018-19: Assessed Project 1 Report
当前位置:以往案例 > >案例R语言F79MA 2018-19: Assessed Project 1 Report
2019-03-15

F79MA 2018-19: Assessed Project 1 Report

Introduction and Summary

This project is aim to look at parameter estimation for an inverse gamma model, which can be considered as a survival distribution. Also, we will show that computer simulation can be a very useful tool to characterize the performance of statistical estimators.

Analysis

The probability density function of inverse distribution is

image.png

where is a shape parameter, is a scale parameter, and is the gamma function.

Let be a random sample from an inverse gamma distribution with shape parameter and unknown scale parameter .

Task 1

As is i.i.d, the joint density function is

image.png

Then

image.png

Differentiating with respect to , we get the likelihood equation

image.png

Differentiating again, we get

image.png

Which is negative for all , and any n.

Solving equation (1.3) for , we get

image.png

image.png

因部分公式不兼容请查看doc文档

sol-F79MA.docx

r code:

rm(list=ls(all=TRUE))
library('invgamma')

# task 3
alpha<-10 beta<-1 #true value ntime<-1000 nmax<-50 betahat<-rep(0,ntime) vac<-rep(0,nmax) MSE<-rep(0,nmax) for (n in 1:nmax){ # repeat times for (t in 1:ntime){ Y<-rinvgamma(n,alpha,rate=1) betahat[t]<-alpha/mean(1/Y) } vac[n]<-var(betahat) MSE[n]<-1/ntime*sum(betahat-beta)^2 } n<-c(1:50) plot(n,MSE,main='MSE vs n') # task 5 crbound<-beta^2/alpha/n^2 plot(MSE-crbound) abline(h=0,col='red') # task 6 betas<-c(0.3,0.5,0.8,1,1.2) MSEMME<-matrix(0,nrow=5,ncol=3) MSEMLE<-matrix(0,nrow=5,ncol=3) ns<-c(20,40,60) for (i in 1:5){ beta<-betas[i] for (j in 1:3){ n<-ns[j] mme<-rep(0,100) mle<-rep(0,100) for (t in 1:100){ sp<- rinvgamma(n,alpha,beta) mme[t]<-mean(sp)*((mean(sp))^2/var(sp)+1) mle[t]<-alpha/mean(1/sp) } MSEMME[i,j]<-1/100*sum((mme-beta)^2) MSEMLE[i,j]<-1/100*sum((mle-beta)^2) } } round(MSEMME,4) round(MSEMLE,4) # task 7 betas<-c(0.3,0.5,0.8,1,1.2) AVsMLE<-matrix(0,nrow=5,ncol=3) ns<-c(20,40,60) for (i in 1:5){ beta<-betas[i] for (j in 1:3){ n<-ns[j] mle<-rep(0,100) for (t in 1:100){ sp<- rinvgamma(n,alpha,beta) mle[t]<-alpha/mean(1/sp) } AVsMLE[i,j]<-sd(mle) } } round(AVsMLE,4)

在线提交订单