# Load all of the previous models fit for this dataset
#
source('prob_2.R')

# Note sockeye.m1 is the Ricker model but without a variance model.
# Here we fit the same model but with a power variance model.
#
sockeye.m3 = gnls(recruits ~ beta1 * spawners * exp(-beta2 * spawners), data = sockeye[-12,], start=list(beta1=2, beta2=0.001), weights=varPower())
recruits_predicted_m3 = predict(sockeye.m3, newdata=data.frame(spawners=spawnersVals))

#postscript("../../WriteUp/Graphics/Chapter6/chap_6_prob_3_data_N_fits.eps", onefile=FALSE, horizontal=FALSE)
plot(sockeye$spawners, sockeye$recruits, xlab='spawners',  ylab='recruits')
lines(spawnersVals, recruits_predicted_m1, col='red')
lines(spawnersVals, recruits_predicted_m2, col='blue')
lines(spawnersVals, recruits_predicted_m3, col='green')
legend('topleft', c('data', 'Ricker model','tbs box-cox model',  'Ricker model (power variance)'), lty=c(1, 1, 1, 1), col=c('black', 'blue', 'red', 'green'))
grid()
#dev.off()

# Lets just make sure the power variance gives normal looking residuals:
#
res = sockeye$recruits[-12] - fitted.values(sockeye.m3) # residuals
sigma = sd(res)

#postscript("../../WriteUp/Graphics/Chapter6/chap_6_prob_3_std_residual_plot.eps", onefile=FALSE, horizontal=FALSE)
standardRes = res/sigma
qqnorm(standardRes, main='')
abline(a=0, b=1)
grid()
#dev.off()

print(shapiro.test(standardRes))