MAT8886 a short word on profile likelihood
[This article was first published on Freakonometrics - Tag - R-english, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Profile likelihood is an interesting theory to visualize and compute confidence interval for estimators (see e.g. Venzon & Moolgavkar (1988)). As we will use is, we will plot
But more generally, it is possible to consider
> base1=read.table( + "http://freakonometrics.free.fr/danish-univariate.txt", + header=TRUE) > library(evir) > X=base1$Loss.in.DKM > u=5
The function to draw the profile likelihood for the tail index parameter is then
> Y=X[X>u]-u > loglikelihood=function(xi,beta){ + sum(log(dgpd(Y,xi,mu=0,beta))) } > XIV=(1:300)/100;L=rep(NA,300) > for(i in 1:300){ + XI=XIV[i] + profilelikelihood=function(beta){ + -loglikelihood(XI,beta) } + L[i]=-optim(par=1,fn=profilelikelihood)$value } > plot(XIV,L,type="l")
It is possible to use it that profile likelihood function to derive a confidence interval,
> PL=function(XI){ + profilelikelihood=function(beta){ + -loglikelihood(XI,beta) } + return(optim(par=1,fn=profilelikelihood)$value)} > (OPT=optimize(f=PL,interval=c(0,3))) $minimum [1] 0.6315989 $objective [1] 754.1115 > up=OPT$objective > abline(h=-up) > abline(h=-up-qchisq(p=.95,df=1),col="red") > I=which(L>=-up-qchisq(p=.95,df=1)) > lines(XIV[I],rep(-up-qchisq(p=.95,df=1),length(I)), + lwd=5,col="red") > abline(v=range(XIV[I]),lty=2,col="red")
This is done with the following code
> library(ismev) > gpd.profxi(gpd.fit(X,5),xlow=0,xup=3)
To leave a comment for the author, please follow the link and comment on their blog: Freakonometrics - Tag - R-english.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.