## code illustrating section 7, on the problem with non-linear
## dynamics + section 8 on statistics

rick <- function(r=3,n0=-1,n.t=40,e=rep(0,n.t)) {
## get Ricker (log scale), given random effects.
  n <- rep(n0,n.t)
  for (i in 2:n.t) n[i] <- r + n[i-1] - exp(n[i-1]) + e[i-1]
  n
}

rick.ll <- function(r,y) {
## deterministic Ricker log likelihood
  n <- rick(r,n.t=length(y))
  sum(dpois(y,exp(n),log=TRUE))
}

rick.fey <- function(r,e,y) {
## stochastic Ricker log(f(e,y))
  n <- rick(r,n.t=length(y),e=e)
  fy.e <- sum(dpois(y,exp(n),log=TRUE)) ## f(y|e)
  fe <- sum(dnorm(e,sd=.1,log=TRUE)) ## f(e)
  fy.e + fe
}

## illustrate the problem with non-linear dynamics
ps <- FALSE
if (ps) postscript("c:/lnotes/stat4biodyn/notes/ricker1.eps",height=4)
par(mfrow=c(1,2),mar=c(6,6,1,1))
plot(1:40,exp(rick(2,n0=log(.5))),type="l",ylab="population",xlab="time",cex.lab=1.8,cex.axis=1.6)
lines(1:40,exp(rick(2,n0=log(.49))),lty=2)
text(38,.7,"r=2",cex=2)
plot(1:40,exp(rick(3,n0=log(.5))),type="l",ylab="population",xlab="time",cex.lab=1.8,cex.axis=1.6)
lines(1:40,exp(rick(3,n0=log(.49))),lty=2)
text(38,6.8,"r=3",cex=2)
if (ps) dev.off()

## simulate some data...
set.seed(0)
ll3 <- ll2  <- r <- seq(1,4,length=1000)
n <- rick(3);y <- rpois(n,exp(n))

## look at transects through log-likeihood...
for (i in 1:length(r)) ll3[i] <- rick.ll(r[i],y)
n <- rick(2);y <- rpois(n,exp(n))
for (i in 1:length(r)) ll2[i] <- rick.ll(r[i],y)
ps <- FALSE
if (ps) postscript("c:/lnotes/stat4biodyn/notes/ricker2.eps",height=4)
par(mfrow=c(1,2),mar=c(6,6,1,1))
par(mfrow=c(1,2))
plot(r,ll2,type="l",cex.lab=1.8,cex.axis=1.6,ylab="log lik")
ind <- ll2==max(ll2);points(r[ind],ll2[ind],pch=20)
text(1.5,-450,expression(paste(r[true],"=2")),cex=2)
plot(r,ll3,type="l",cex.lab=1.8,cex.axis=1.6,ylab="log lik")
ind <- ll3==max(ll3);points(r[ind],ll3[ind],pch=20)
text(1.5,-700,expression(paste(r[true],"=3")),cex=2)
if (ps) dev.off()

## stochastic Ricker model...
ps <- FALSE
if (ps) postscript("c:/lnotes/stat4biodyn/notes/ricker3.eps",height=4)
par(mfrow=c(1,2),mar=c(6,6,3,1))

## simulate data with r=2....
e <- rnorm(40)*.1
r <- 2
n <- rick(r)
y <- rpois(n,exp(n))

## transect through log(f(e,y)), as e[1] varies
fey <- e1 <- seq(-.2,.2,length=1000)
for (i in 1:length(e1)) {
 e[1] <- e1[i]
 fey[i] <- rick.fey(r,e,y)
}
plot(e1,fey,type="l",ylab="log{f(e,y)}",xlab=expression(e[1]),
     cex.lab=1.8,cex.axis=1.6,main=expression(paste(r[true],"=2")),cex.main=1.6)

## simulate data with r=3....
r <-3
n <- rick(r)
y <- rpois(n,exp(n))

## corresponding transect through log(f(e,y)), as e[1] varies
fey <- e1 <- seq(-.2,.2,length=1000)
for (i in 1:length(e1)) {
 e[1] <- e1[i]
 fey[i] <- rick.fey(r,e,y)
}
plot(e1,fey,type="l",ylab="log{f(e,y)}",xlab=expression(e[1]),
     cex.lab=1.8,cex.axis=1.6,main=expression(paste(r[true],"=3")),cex.main=1.6)
if (ps) dev.off()