Rubén Roa-Ureta
2008-May-06 17:34 UTC
[R] Solving for parameter values at likelihood points
ComRades, On Oct 16, 2006, I posted a question regarding how to find the parameter values that made the likelihood function take a value equal to 1/8th of the maximum likelihood value. There was no reply but I found a solution and would like to know if there is a better solution using the function uniroot(). I think with uniroot() I can get better parameter values at these bounds. This is my solution, using the function nearest() from package GenKern. The example concerns the marginal likelihood for the normal variance parameter when the mean is a nuisance parameter (see section 7.3 in Royall, 1997, Statistical evidence. A likelihood paradigm, Chapman & Hall). Rub?n y<-c(-2.250357,10.723614,7.238959,9.179939,5.831603,9.301580,4.019388,9.777280,5.587761,3.600276,9.753381,5.154928,2.776350,9.940350,14.185712,4.281198,10.994180,9.333027,13.082535,2.067826,10.772551,-3.811529,5.446021,-5.333287,1.421544) n <- length(y) s.sq <- (1/(n-1))*sum((y-mean(y))^2) sigma.sq <- seq(1,100,0.1) Lik.M <- sigma.sq^(-(n-1)/2)*exp(-(n-1)*s.sq/(2*sigma.sq)) Rel.Lik.M <- Lik.M/max(Lik.M) plot(sigma.sq,Rel.Lik.M,type="l") abline(h=1/8) arrows(y0=0.4,y1=0.18,x0=0,x1=12) arrows(y0=0.4,y1=0.18,x0=62,x1=50) library(GenKern) cbind(sigma.sq,Rel.Lik.M)[Rel.Lik.M==1] #maximum likelihood estimate max.idx <- nearest(Rel.Lik.M,1,outside=TRUE) # vector index value at the maximum likelihood #finding the parameter value at 1/8th the max. likelihood to the left of the maximum cbind(sigma.sq[1:max.idx],Rel.Lik.M[1:max.idx])[nearest(Rel.Lik.M[1:max.idx],0.125),] #finding the parameter value at 1/8th the max. likelihood to the right of the maximum cbind(sigma.sq[max.idx:length(sigma.sq)],Rel.Lik.M[max.idx:length(sigma.sq)])[nearest(Rel.Lik.M[max.idx:length(sigma.sq)],0.125),]