Lecture7 <- function()
{
 Fig7a(1)    # Example I 
# Fig7a(10)   # Example I 

# Fig7b()      # Example I - sampled from normal importance function
 
# Fig7c()     # Lab exercise 1

 Fig7d(T)    # Lab exercise 2
# Fig7d(F)     # lab exercise 2

}

# =============================================================================
Fig7d <- function(IsInit)
{
 FileName <- "I:\\Lecture7Dat.csv"
 TheData <- matrix(scan(FileName,sep=","),ncol=3,byrow=T)
 print(TheData)
 Catch <- TheData[,2]
 CPUE <- TheData[,3]
 LikeModel(Catch,CPUE,10000,0.205,1.0)

 Threshold <- exp(30.34)
 Nout <- 200
 Ntotal <- 0
 Cumu <- 0
 Ndone <- 0
 AveLike <- 0
 Vals <- matrix(0,nrow=Nout,ncol=4)
 while (Ndone < Nout)
  {
   r <- runif(1,0.15,0.45)
   K <- runif(1,5000,8000)
   if (IsInit == T)
    phi <- runif(1,0.8,1.2)
   else
    phi <- 1 
   outs <- LikeModel(Catch,CPUE,K,r,phi)
   Like <- exp(outs$LogLike)
   Depl <- outs$Depl
   Cumu <- Cumu + Like
   AveLike <- AveLike + Like
   Ntotal <- Ntotal + 1
   while (Cumu > Threshold)
    {
     Cumu <- Cumu - Threshold
     Ndone <- Ndone + 1
     Vals[Ndone,1] <- r
     Vals[Ndone,2] <- K
     Vals[Ndone,3] <- phi
     Vals[Ndone,4] <- Depl
    }
  }
 print(AveLike/Ntotal)
 par(mfrow=c(2,2))
 hist(Vals[,1],xlab="r",main="") 
 hist(Vals[,2],xlab="K",main="") 
 hist(Vals[,3],xlab="Phi",main="") 
 hist(Vals[,4],xlab="Depletion",main="") 
}

LikeModel <- function(Catch,CPUE,K,r,phi)
{
 Ndata <- length(Catch)
 
 Pop <- rep(0,length=length(Catch))
 Pop[1] <- phi*K
 for (I in 2:length(Catch))
   Pop[I] <- max(0.02,Pop[I-1]+r*Pop[I-1]*(1-Pop[I-1]/K)-Catch[I-1])
 
 lnq <- 0
 for (I in 1:length(Catch))
  lnq <- lnq + log(CPUE[I]/Pop[I])
 q <- exp(lnq/length(Catch)) 
 
 SS <- 0
 for (I in 1:length(Catch))
  SS <- SS + (log(CPUE[I])-log(q*Pop[I]))^2
 Sigma <- sqrt(SS/length(Catch))
 
 LogLike <- -1*(Ndata*log(Sigma)+Ndata/2.0)
 
 Outs <- NULL
 Outs$LogLike <- LogLike
 Outs$Depl <- Pop[Ndata]/K
 return(Outs)

}

# =============================================================================
Fig7c <- function()
{
 Lengths <- c(10,20,30,40,50,60,70,80,90,100)
 Obs <- c(0,0,0,5,14,18,20,19,20,20)

 Nout <- 10000  # set to 1000 for class
 
 # This is, of course, the trick!
 Threshold <- exp(-37.6)
 Cumu <- 0
 Ndone <- 0
 ProbMore <- 0
 Vals <- rep(0,length=Nout)
 while (Ndone < Nout)
  {
   Len50 <- runif(1,20,60)
   Delta <- runif(1,0,20)
   Like <- exp(Likelihood(Lengths,Obs,20,Len50,Delta))
   Cumu <- Cumu + Like
   while (Cumu > Threshold)
    {
     Cumu <- Cumu - Threshold
     Ndone <- Ndone + 1
     Vals[Ndone] <- Len50
     Len80 <- Len50 -1*log(1/4)/log(19)*Delta
     if (Len80 > 55) ProbMore <- ProbMore + 1/Nout
    }
  }
 print(ProbMore)
 hist(Vals) 
 
}

# =============================================================================

Fig7b <- function()
{
 Lengths <- c(10,20,30,40,50,60,70,80,90,100)
 Obs <- c(1,4,9,38,98,160,192,196,199,200)
 print(Obs)
 
 L50 <- 40+seq(from=1,to=99)/100*20
 DDD <-  0+seq(from=1,to=99)/100*40
 
 par(mfrow=c(1,1))

 # grid search method
 Plot2 <- matrix(0,ncol=99,nrow=99)
 for (Irow in 1:99)
  for (Icol in 1:99)
   Plot2[Irow,Icol] <- as.double(Likelihood(Lengths,Obs,200,L50[Irow],DDD[Icol]))
 ymax <- max(Plot2)
 Plot2 <- Plot2 -ymax  
 Plot2 <- exp(Plot2)
 print(attributes(Plot2))
 res<-persp(x=L50,y=DDD,z=Plot2,xlab=expression(L[50]),ylab=expression(delta),zlab="Density",zlim=c(0,1.1),r=8,theta=45)
 
 # Now use SIR
 res<-persp(x=L50,y=DDD,z=Plot2,xlab=expression(L[50]),ylab=expression(delta),zlab="Density",zlim=c(0,1.1),r=8,theta=45)
 Nsim <- 10000
 Len50R <- rnorm(Nsim,50,2)
 DeltaR <- rnorm(Nsim,20,3)
 P1 <- dnorm(Len50R,50,2)
 P2 <- dnorm(DeltaR,20,3)
 LogLik <- rep(0,length=Nsim)
 for (II in 1:Nsim)
  LogLik[II] <- exp(Likelihood(Lengths,Obs,200,Len50R[II],DeltaR[II])-ymax)
 Probs <- LogLik/P1/P2
 points(trans3d(Len50R,DeltaR,LogLik,res),pch=16,col="red",cex=0.5) 
 
 # Plot resample only
 N2 <- 1000
 xx <- sample(seq(from=1,to=Nsim),N2,replace=T,prob=Probs)
 persp(x=L50,y=DDD,z=Plot2,xlab=expression(L[50]),ylab=expression(delta),zlab="Density",zlim=c(0,1.1),r=8,theta=45) 
 points(trans3d(Len50R[xx],DeltaR[xx],LogLik[xx],res),pch=16,col="red",cex=0.5) 
 print(max(table(Len50R[xx])))
 


}

# =============================================================================

Fig7a <- function(Mult)
{

 Lengths <- c(10,20,30,40,50,60,70,80,90,100)
 Obs <- c(0,0,0,5,14,18,20,19,20,20)*Mult
 print(Obs)
 
 L50 <- 40+seq(from=1,to=99)/100*20
 DDD <-  0+seq(from=1,to=99)/100*40
 
 par(mfrow=c(1,1))

 # grid search method
 Plot2 <- matrix(0,ncol=99,nrow=99)
 for (Irow in 1:99)
  for (Icol in 1:99)
   Plot2[Irow,Icol] <- as.double(Likelihood(Lengths,Obs,20*Mult,L50[Irow],DDD[Icol]))
 ymax <- max(Plot2)
 Plot2 <- Plot2 -ymax  
 Plot2 <- exp(Plot2)
 print(attributes(Plot2))
 res<-persp(x=L50,y=DDD,z=Plot2,xlab=expression(L[50]),ylab=expression(delta),zlab="Density",zlim=c(0,1.1),r=8,theta=45)
 
 # Now use SIR
 res<-persp(x=L50,y=DDD,z=Plot2,xlab=expression(L[50]),ylab=expression(delta),zlab="Density",zlim=c(0,1.1),r=8,theta=45)
 Nsim <- 10000
 Len50R <- runif(Nsim,40,60)
 DeltaR <- runif(Nsim,0,40)
 LogLik <- rep(0,length=Nsim)
 for (II in 1:Nsim)
  LogLik[II] <- exp(Likelihood(Lengths,Obs,20*Mult,Len50R[II],DeltaR[II])-ymax)
 # Plot ALL points 
 points(trans3d(Len50R,DeltaR,LogLik,res),pch=16,col="red",cex=0.5) 
 
 # Plot resample only
 N2 <- 1000
 xx <- sample(seq(from=1,to=Nsim),N2,replace=T,prob=LogLik)
 persp(x=L50,y=DDD,z=Plot2,xlab=expression(L[50]),ylab=expression(delta),zlab="Density",zlim=c(0,1.1),r=8,theta=45) 
 points(trans3d(Len50R[xx],DeltaR[xx],LogLik[xx],res),pch=16,col="red",cex=0.5) 
 print(max(table(Len50R[xx])))

 hist(Len50R[xx])
 hist(DeltaR[xx])
 AA

}

# ==================================================================================

Likelihood <- function(Lengths,Obs,NN,Len50,Delta)
{
 Pred <- 1/(1+exp(-log(19)*(Lengths-Len50)/Delta))
 Log_like <- Obs[1]*log(Pred[1])+(NN-Obs[1])*log(1-Pred[1])
 Log_like <- Obs*log(Pred)+(NN-Obs)*log(1-Pred)
 LogL <- sum(Log_like)
 if (is.na(LogL)) LogL <- -1000
 return(LogL)

}

# ==================================================================================

Lecture7()