aspmrun.r

text/r-latex aspmrun.r — 1.1 KB

File contents

#filename="aspmrun.r"
# Download this file from http://vr3pc109.rhi.hi.is/fish/fish5108statass/lecture60/aspmrun.r
# Download the other files below from the same directory
# 
# To do a complete assessment give the command:
# source("aspmrun.r")
source("aspminit.r")
source("aspm.r")
source("ssefcn.r")
# Get a better initial value of q...
fit0<-aspm(Fvec,Rvec,totyrs,M,w)
q<-mean(I)/mean(fit0$Bhat)
parameters.0<-c(log(Fvec),log(Rvec),log(q))     # Initial values of all parameters
fm<-nlm(ssefcn,parameters.0,iterlim=500,Y=Y,I=I,M=M,w=w)
parameters<-fm$estimate
Fvec<-exp(parameters[1:totyrs])
Rvec<-exp(parameters[(totyrs+1):(2*totyrs)])
q<-exp(parameters[2*totyrs+1])
fmopt<-aspm(Fvec,Rvec,totyrs,M,w)
Yhat<-fmopt$Yhat
Bhat<-fmopt$Bhat
Ihat<-q*Bhat

time<-1:totyrs
par(mfrow=c(2,3))
plot(yrs,Y,xlab="Year",ylim=c(0,max(Y)))
lines(yrs,Yhat)
plot(yrs,I,xlab="Year",ylim=c(0,max(I)))
lines(yrs,Ihat)
plot(yrs,Rvec,xlab="Year",ylab="Recruitment",ylim=c(0,max(Rvec)))
plot(yrs,Fvec,xlab="Year",ylab="F",ylim=c(0,max(Fvec)))
plot(yrs,Bhat,xlab="Year",ylab="Biomass",ylim=c(0,max(Bhat)))