# NOT RUN {
library(mvtnorm)
library(MASS)
data(abdat)
param <- log(c(r= 0.42,K=9400,Binit=3400,sigma=0.05))
bestmod <- nlm(f=negLL,p=param,funk=simpspm,indat=abdat,
hessian=TRUE,logobs=log(abdat$cpue),
typsize=magnitude(param),iterlim=1000)
optpar <- bestmod$estimate
vcov <- solve(bestmod$hessian) # solve inverts matrices
columns <- c("r","K","Binit","sigma")
N <- 1000 # the contours improve as N increases; try 5000
mvnpar <- matrix(exp(rmvnorm(N,mean=optpar,sigma=vcov)),
nrow=N,ncol=4,dimnames=list(1:N,columns))
xv <- mvnpar[,"K"]
yv <- mvnpar[,"r"]
# plotprep(width=6,height=5,newdev = FALSE)
plot(xv,yv,type="p") # use default 0.5 and 0.9 contours.
addcontours(xv,yv,range(xv),range(yv),lwd=2,col=2)
points(mean(xv),mean(yv),pch=16,cex=1.5,col=2)
# }
Run the code above in your browser using DataLab