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()