#==========================================================================================
# Aula 19 da disciplina ce223 (24/05/2011)
# Mais sobre a optim()
#                                                               Professor Walmes M. Zeviani
#                                                                     www.leg.ufpr.br/ce223
#==========================================================================================

#------------------------------------------------------------------------------------------
# ajuste de modelos de regressão

# potássio liberado acumulado para o esterco de codorna ----

klib <- c(51.03, 57.76, 26.60, 60.65, 87.07, 64.67,
          91.28, 105.22, 72.74, 81.88, 97.62, 90.14,
          89.88, 113.22, 90.91, 115.39, 112.63, 87.51,
          104.69, 120.58, 114.32, 130.07, 117.65, 111.69,
          128.54, 126.88, 127.00, 134.17, 149.66, 118.25,
          132.67, 154.48, 129.11, 151.83, 147.66, 127.30)

# tempo em que foram feitas as coletas ---------------------

tempo <- rep(c(15, 30, 45, 60, 75, 90,
               120, 150, 180, 210, 240, 270), each=3)

# data.frame com todos os dados ----------------------------

liber <- data.frame(tempo, klib)
str(liber)

#------------------------------------------------------------------------------------------
# gráfico dos dados com adição de uma curva estimada visualmente

plot(klib~tempo, data=liber)
A <- 120; V <- 20; D <- 0.05
curve(A*x/(V+x)+D*x, add=TRUE)

#------------------------------------------------------------------------------------------
# escrevendo a função de verossimilhança segundo o modelo de regressão

fun.ver <- function(parms, k, t){
  ## params 1, 2, 3 e 4 correspondem a A, V, D, e s2
  ## k são valores observados de potássio, t é o tempo
  ## função retorna a log-verossimilhança
  n <- length(k)
  ft <- parms[1]*t/(parms[2]+t)+parms[3]*t
  ll <- -(n/2)*log(2*pi)-(n/2)*log(parms[4])-(0.5/parms[4])*sum((k-ft)^2)
  return(ll)
}

fun.ver(c(100, 20, 0.05, 1), k=liber$klib, t=liber$tempo)

#------------------------------------------------------------------------------------------
# usando a optim() com a função de verossimilhança, chutes da curva e dados

A <- 117; V <- 23; D <- 0.12; s2 <- 11^2
op <- optim(c(A, V, D, s2), fun.ver, k=liber$klib, t=liber$tempo,
            control=list(fnscale=-1))
str(op)

A <- op$par[1]; V <- op$par[2]; D <- op$par[3];
plot(klib~tempo, data=liber)
curve(A*x/(V+x)+D*x, add=TRUE, col=2)

#------------------------------------------------------------------------------------------
# minimizando a soma de quadrados dos resíduos

fun.sqr <- function(parms, k, t){
  ## params 1, 2, e 3 correspondem a A, V, D
  ## k são valores observados de potássio, t é o tempo
  ## função retorna a soma de quadrados dos resíduos
  ft <- parms[1]*t/(parms[2]+t)+parms[3]*t
  sqr <- sum((k-ft)^2)
  return(sqr)
}

op <- optim(c(120, 23, 0.12), fun.sqr, k=liber$klib, t=liber$tempo)
str(op)
op$par

#------------------------------------------------------------------------------------------
# usando a função nls() nonlinear leasts squares do R para minimizar a sqr

n0 <- nls(klib~A*tempo/(V+tempo)+D*tempo, data=liber, start=c(A=120, V=20, D=0.05))
summary(n0)
coef(n0)

#------------------------------------------------------------------------------------------
# compare as estimativas dos 3 procedimentos, elas devem ser iguais

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