[R] Understanding R code

kfcnhl zhengchenji18 at hotmail.com
Thu Aug 20 22:31:32 CEST 2009


What is 
1. par.ests <- optimfit$par
2. fisher <- hessb(negloglik, par.ests, maxvalue=maxima);
3. varcov <- solve(fisher);
4. par.ses <- sqrt(diag(varcov));

Thanks a lot,

fit.GEV <- function(maxima)
{
sigma0 <- sqrt((6. * var(maxima))/pi)
mu0 <- mean(maxima) - 0.57722 * sigma0
xi0 <- 0.1
theta <- c(xi0, mu0, sigma0)

#10/5/2007: removed assign() for maxima.nl
#10/5/2007: passed additional parameter to avoid using assign()
negloglik <- function(theta, maxvalue)
{
-sum(dGEV(maxvalue,theta[1],theta[2],abs(theta[3]),logvalue=TRUE)); 
}

#10/5/2007: passed additional 4th parameter which will be passed to
negloglik
optimfit <- optim(theta, fn=negloglik, gr=NULL, maxvalue=maxima,
method="BFGS");
par.ests <- optimfit$par;
if(optimfit$convergence == 0) #if code is 0 we are OK; 
converged <- TRUE
else
converged <- FALSE;
#replace 3rd parameter estimate with its absolute value:
par.ests[3] <- abs(par.ests[3]);
#10/5/2007: added final parameter which must be passed to negloglik
fisher <- hessb(negloglik, par.ests, maxvalue=maxima);
varcov <- solve(fisher);

par.ses <- sqrt(diag(varcov));
out <- list(par.ests = par.ests, par.ses = par.ses, varcov = varcov,
converged = 
#10/5/2007: passed additional 2nd parameter to -negloglik()
converged, llmax = -negloglik(par.ests,maxima));
#06/30/2008: the order of the names is incorrect; note that the parameter
vector theta 
#(passed to the optimization function) has the order c(xi0, mu0, sigma0);
the order of 
#the return parameters should be identical;
#hence we need to make the following correction to maintain the same order
#names(out$par.ests) <- c("xi", "sigma", "mu");
#names(out$par.ses) <- c("xi", "sigma", "mu");
names(out$par.ests) <- c("xi", "mu", "sigma");
names(out$par.ses) <- c("xi", "mu", "sigma");
out;
}
-- 
View this message in context: http://www.nabble.com/Understanding-R-code-tp25069337p25069337.html
Sent from the R help mailing list archive at Nabble.com.




More information about the R-help mailing list