Lecture2 <- function()
{
# Part1()
 Ex1()

}

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

Part1a <- function(Propn,Plot=F)
{
 TheData <- matrix(scan("C:\\Courses\\FISH507G_10\\Lectures\\Lect3a.dat",skip=1),byrow=T,ncol=2)
 TheData <- TheData[,2]

 Blow1 <- seq(from=1,to=40,by=1)
 Blow2 <- seq(from=1,to=60,by=1)
 Likes <- matrix(0,nrow=length(Blow1),ncol=length(Blow2))
 for (ii in 1:length(Blow1))
  for (jj in 1:length(Blow2))
   {
   LikeOut <- Propn*dpois(TheData,Blow1[ii])+(1-Propn)*dpois(TheData,Blow2[jj])
   Likes[ii,jj] <- -1*sum(log(LikeOut))
  }
 if (Plot==T)
  {
   par(mfrow=c(1,1))
   persp(x=Blow1,y=Blow2,z=Likes)
  } 
 Ibest <- which(Likes==min(Likes),arr.ind=T)

 cat(Propn,Blow1[Ibest[1]],Blow2[Ibest[2]],min(Likes),"\n")
 return(min(Likes))

}

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

Part1 <- function()
{
 Part1a(0.5,Plot=T)

 Propns <- seq(from=0.05,to=0.5,by=0.025)
 Values <- rep(0,length(Propns))
 for (ii in 1:length(Propns))
  Values[ii] <- Part1a(Propns[ii],Plot=F)
  
 par(mfrow=c(2,2))
 plot(Propns,Values,xlab="Fraction",ylab="Negative log-likelihood",type='b',pch=16,lty=1,lwd=2)
 Ibest <- which(Values==min(Values),arr.ind=T)
 print(cbind(Propns,Values))
 print(Propns[Ibest])
}


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

Ex1a <- function(TheData,Sigma,Plot=F)
{
 Linf <- seq(from=50,to=150,by=1)
 Kappa <- seq(from=0.05,to=0.5,by=0.01)
 Likes <- matrix(0,nrow=length(Linf),ncol=length(Kappa))
 for (ii in 1:length(Linf))
  for (jj in 1:length(Kappa))
   {
   LikeOut <- dnorm(TheData[,3],Linf[ii]*(1-exp(-Kappa[jj]*(TheData[,2]))),Sigma)
   Likes[ii,jj] <- -1*sum(log(LikeOut))
  }
 Ibest <- which(Likes==min(Likes),arr.ind=T)

 cat(Sigma,Linf[Ibest[1]],Kappa[Ibest[2]],min(Likes),"\n")

 if (Plot==T)
  {
   par(mfrow=c(2,2))
#   persp(x=Linf,y=Kappa,z=Likes)
   filled.contour(x=Linf,y=Kappa,z=Likes,nlevels=100,color = terrain.colors)
   plot(TheData[,2],TheData[,3],xlab="Age",ylab="Length",ylim=c(0,max(TheData[,3])))
   ages <- seq(from=0,to=100)
   lens <- Linf[Ibest[1]]*(1-exp(-Kappa[Ibest[2]]*ages))
   lines(ages,lens,lty=1,lwd=2)
  } 

 return(min(Likes))

}

# ------------------------------------------------------------------------------------------

Ex1 <- function()
{
 TheData <- matrix(scan("C:\\Courses\\Fish507G_10\\Lectures\\Lect3b.dat",skip=1),byrow=T,ncol=3)
 TheData <- TheData[TheData[,1]==1,]
 
 Ex1a(TheData,9,Plot=T)
}

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

Lecture2()