library(stats4)

# =============================================================================================================================================
Ex1a <- function()
{
 PlotRes <- function(Linf1,Kappa1,Linf2,Kappa2,Sigma,Age,Length,Sex)
  {
   Obs1 <- Length[Sex==1]
   Age2 <- seq(from=0,to=30,by=0.2)
   Pred1 <- Linf1*(1-exp(-Kappa1*Age2))
   Obs2 <- Length[Sex==2]
   Pred2 <- Linf2*(1-exp(-Kappa2*Age2))
   ymax <- max(Pred1,Pred2,Obs1,Obs2)*1.05
   par(mfrow=c(2,2))
   plot(Age[Sex==1],Obs1,pch=16,ylab="Length",xlab="Age",ylim=c(0,ymax))
   lines(Age2,Pred1,lty=1,lwd=2)
   plot(Age[Sex==2],Obs2,pch=16,ylab="Length",xlab="Age",ylim=c(0,ymax))
   lines(Age2,Pred2,lty=1,lwd=2)
  }
  
 # negative log-likelihood
 LL4 <- function(Linf1,Kappa1,Linf2,Kappa2,Sigma,Age,Length,Sex)
 {
  Obs1 <- Length[Sex==1]
  Pred1 <- Linf1*(1-exp(-Kappa1*Age[Sex==1]))
  Obs2 <- Length[Sex==2]
  Pred2 <- Linf2*(1-exp(-Kappa2*Age[Sex==2]))
  Like1 <- dnorm(Obs1,Pred1,Sigma)
  Like2 <- dnorm(Obs2,Pred2,Sigma)
  LikeOut <- -1*(sum(log(Like1))+sum(log(Like2)))
  return(LikeOut)
 }

 # Read the data
 TheData <- matrix(scan("C:\\Courses\\FISH507G_10\\Lectures\\Lect3b.dat",skip=1),byrow=T,ncol=3)
 Sex <- TheData[,1]
 Age <- TheData[,2]
 Length <- TheData[,3]
 
 # Fit the model and plot
 ret4 <- mle(LL4,start=list(Linf1=100,Kappa1=0.2,Linf2=100,Kappa2=0.2,Sigma=10),fixed=list(Age=Age,Length=Length,Sex=Sex),lower=c(0.01,0.01,0.01,0.01,0.01),upper=c(Inf,Inf,Inf,Inf,Inf),method="L-BFGS-B")
 print(summary(ret4))
 ret4 <- coef(ret4)
 PlotRes(ret4[1],ret4[2],ret4[3],ret4[4],ret4[5],Age,Length,Sex)


}
# =============================================================================================================================================
Ex1b <- function()
{

 PlotRes <- function(Linf1,Kappa1,Linf2,Kappa2,Sigma,Age,Length,Sex)
  {
   Obs1 <- Length[Sex==1]
   Age2 <- seq(from=0,to=30,by=0.2)
   Pred1 <- Linf1*(1-exp(-Kappa1*Age2))
   Obs2 <- Length[Sex==2]
   Pred2 <- Linf2*(1-exp(-Kappa2*Age2))
   ymax <- max(Pred1,Pred2,Obs1,Obs2)*1.05
   par(mfrow=c(2,2))
   plot(Age[Sex==1],Obs1,pch=16,ylab="Length",xlab="Age",ylim=c(0,ymax))
   lines(Age2,Pred1,lty=1,lwd=2)
   plot(Age[Sex==2],Obs2,pch=16,ylab="Length",xlab="Age",ylim=c(0,ymax))
   lines(Age2,Pred2,lty=1,lwd=2)
  }

 #negative loglikelihood (note there is only one Linf and kappa)
 LL4a <- function(Linf,Kappa,Sigma,Age,Length,Sex)
 {
  Obs1 <- Length[Sex==1]
  Pred1 <- Linf*(1-exp(-Kappa*Age[Sex==1]))
  Obs2 <- Length[Sex==2]
  Pred2 <- Linf*(1-exp(-Kappa*Age[Sex==2]))
  Like1 <- dnorm(Obs1,Pred1,Sigma)
  Like2 <- dnorm(Obs2,Pred2,Sigma)
  LikeOut <- -1*(sum(log(Like1))+sum(log(Like2)))
  return(LikeOut)
 }

 # Read the data
 TheData <- matrix(scan("C:\\Courses\\FISH507G_10\\Lectures\\Lect3b.dat",skip=1),byrow=T,ncol=3)
 Sex <- TheData[,1]
 Age <- TheData[,2]
 Length <- TheData[,3]
 
 # Fit the model and show the fir
 ret4 <- mle(LL4a,start=list(Linf=100,Kappa=0.2,Sigma=10),fixed=list(Age=Age,Length=Length,Sex=Sex),lower=c(0.01,0.01,0.01,0.01,0.01),upper=c(Inf,Inf,Inf,Inf,Inf),method="L-BFGS-B")
 print(summary(ret4))
 ret4 <- coef(ret4)
 PlotRes(ret4[1],ret4[2],ret4[1],ret4[2],ret4[3],Age,Length,Sex)
}

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