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
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
Then
Differentiating with respect to , we get the likelihood equation
Differentiating again, we get
Which is negative for all , and any n.
Solving equation (1.3) for , we get
因部分公式不兼容请查看doc文档
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)