##########################################################################################
## Descarga de bateria NHS ###############################################################
##########################################################################################

## Lendo a base de dados
load("aq2.RData")
aq2 <- as.data.frame(aq2)
names(aq2) <- c("Bat", "y", "x")
table(aq2$Bat)
aq2$tempo <- c(1:2001,1:1166, 1:2225, 1:2069, 1:2276)
## Graficos descritivos
require(lattice)

xyplot(y ~ tempo , group = Bat, type="l", data = aq2, auto.key=TRUE)

## Vou tirar o final de cada serie 5% das observações no final de cada série
n <- length(serie[[1]]$y)
retirada <- round(n*0.04,0)
efe <- n - retirada
efetiva <- serie[[1]][1:efe,]
plot(log(efetiva$y) ~  efetiva$tempo)
plot(efetiva$y ~  log(efetiva$tempo))

## Ajustando um modelo não-linear
tt = nls(y ~ A1/(1+exp(-(log(tempo) - Xo)/S1)) + (A1 - A2)/(1 + exp( -(log(tempo) - X02)/S2)), start = list(A1 = 12.7, Xo = 2, S1 = 0.05, A2 = 10.5, X02 = 7, S2=0.5), data= efetiva)

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

require(rpanel)
source("rp.nlregression.R")
load("efetiva.RData")

#save.image("efetiva.RData")
#efetiva$x <- log(efetiva$tempo)

model <- y~12.8/(1+exp(-(x-th2)/th3))+(th4-12.8)/(1+exp(-(x-th5)/th1))
start <- list(th1=c(init=0.5, from=0.05, to=3),
              th2=c(init=2, from=1, to=3),
              th3=c(init=0.05, from=0.005, to=0.5),
              th4=c(init=10, from=5, to=15),
              th5=c(init=8, from=6, to=10))

rm("ajuste")
rp.nlr(model=model, data=efetiva, start=start)

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

model <- y~th1/(1+exp((x-th2)/th3))
start <- list(th1=c(init=12.8, from=10, to=15),
              th2=c(init=2.4, from=2, to=3),
              th3=c(init=0.01, from=0.001, to=0.1))

rm("ajuste")
rp.nlr(model=model, data=efetiva, start=start)

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

model <- y~th1/(1+exp((x-th2)/th3))
start <- list(th1=c(init=10.8, from=8, to=12),
              th2=c(init=8, from=6, to=10),
              th3=c(init=0.6, from=0.1, to=2))

rm("ajuste")
rp.nlr(model=model, data=efetiva, start=start)

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

model <- y~(th3+2)/(1+exp((x-th1)/th2))+(th3)/(1+exp((x-th4)/th5))
start <- list(th1=c(init=2.6, from=2, to=3),
              th2=c(init=0.1, from=0.01, to=0.5),
              th3=c(init=10, from=2, to=18),
              th4=c(init=9.2, from=6, to=12),
              th5=c(init=0.62, from=0.2, to=1))

rm("ajuste")
rp.nlr(model=model, data=efetiva, start=start)

n0 <- nls(y~(th0-th3)/(1+exp((x-th1)/th2))+(th3)/(1+exp((x-th4)/th5)),
          data=efetiva,
          start=list(th0=12.8, th1=2.6, th2=0.1, th3=10.8, th4=9.2, th5=0.62))
coef(n0)

plot(y~x, efetiva, ylim=c(0,13), xlim=c(0,12))
with(as.list(coef(n0)),
             curve((th0-th3)/(1+exp((x-th1)/th2))+(th3)/(1+exp((x-th4)/th5)),
                   add=TRUE, col=2))

qqnorm(residuals(n0))
acf(residuals(n0))
pacf(residuals(n0))

nrow(efetiva)

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

efetiva2 <- efetiva[c(1:3, seq(4, nrow(efetiva), by=10)),]
plot(y~x, efetiva2, ylim=c(0,13), xlim=c(0,12))

n1 <- nls(y~(th0-th3)/(1+exp((x-th1)/th2))+(th3)/(1+exp((x-th4)/th5)),
          data=efetiva2,
          start=list(th0=12.8, th1=2.6, th2=0.1, th3=10.8, th4=9.2, th5=0.62))
coef(n1)

plot(y~x, efetiva2, ylim=c(0,13), xlim=c(0,12))
with(as.list(coef(n1)),
             curve((th0-th3)/(1+exp((x-th1)/th2))+(th3)/(1+exp((x-th4)/th5)),
                   add=TRUE, col=2))

qqnorm(residuals(n1))
acf(residuals(n1))
pacf(residuals(n1))