#rm(list=ls())
par(mfrow=c(1,1), mar=c(3.5,3.5,0.5,0.5), mgp=c(2, 1, 0))
par.ori <- par(no.readonly=TRUE)

dat04 <- read.table("http://www.leg.ufpr.br/~paulojus/dados/data04.txt", head=TRUE)
nd <- data.frame(x = seq(0.5, 3.5, length=201))

## reta única
plot(y ~ x, data = dat04, ylim=c(0, 55))
m1 <- lm(y ~ x, data = dat04)
abline(m1)
nd$m1 <- predict(m1, newdata = nd)
with(nd, lines(m1 ~ x, col = 4))

## Ponto de corte em 2.7
xc <- 2.7
dat04 <- transform(dat04,
                   ind = as.numeric(x > xc))
nd <- transform(nd,
                ind = as.numeric(x >= xc))

## retas paralelas
(m6 <- lm(y ~ x + ind, data = dat04))
nd$m6 <- predict(m6, newdata = nd)
plot(y ~ x, data = dat04, ylim=c(0, 55))
with(subset(nd, x < xc), lines(m6 ~ x, col = 4))
with(subset(nd, x > xc), lines(m6 ~ x, col = 4))

## retas com mesmo intercepto
(m7 <- lm(y ~ x + I(x*ind), data = dat04))
nd$m7 <- predict(m7, newdata = nd)
plot(y ~ x, data = dat04, ylim=c(0, 55))
with(subset(nd, x < xc), lines(m7 ~ x, col = 4))
with(subset(nd, x > xc), lines(m7 ~ x, col = 4))

plot(y ~ x, data = dat04, ylim=c(-10, 55), xlim = c(0, 3.5))
abline(coef(m7)[1], coef(m7)[2], lty=3)
abline(coef(m7)[1], sum(coef(m7)[2:3]), lty=3)
with(subset(nd, x < xc), lines(m7 ~ x, col = 4))
with(subset(nd, x > xc), lines(m7 ~ x, col = 4))

## retas "livres"
(m8 <- lm(y ~ x + ind + I(ind*x), data = dat04))
nd$m8 <- predict(m8, newdata = nd)
plot(y ~ x, data = dat04, ylim=c(0, 55))
with(subset(nd, x < xc), lines(m8 ~ x, col = 4))
with(subset(nd, x > xc), lines(m8 ~ x, col = 4))

## retas "conectadas"
(m9 <- lm(y ~ x + I(ind*(x-xc)), data = dat04))
nd$m9 <- predict(m9, newdata = nd)
plot(y ~ I(x), data = dat04, ylim=c(0, 55))
with(nd, lines(m9 ~ x, col = 4))


##
## "Estimando" o ponto de corte
##
xc.seq <- seq(1, 3.0, by = 0.01)
sigma.f <- function(corte, dados){
    modelo <-  lm(y ~ x + I(ind*(x-corte)), data = dados)
    return(summary(modelo)$sigma)
}
sigma.seq <- sapply(xc.seq, sigma.f, dados = dat04)
plot(xc.seq, sigma.seq, type="l")
(xc.est <- xc.seq[which.min(sigma.seq)])

(m10 <- lm(y ~ x + I((x > xc.est)*(x-xc.est)), data = dat04))
nd$m10 <- predict(m10, newdata = nd)
plot(y ~ x, data = dat04, ylim=c(0, 55))
with(nd, lines(m10 ~ x, col = 4))


## 1as Derivadas analíticas (com coeficientes) e numéricas 
plot(y ~ x, data = dat04, ylim=c(0, 55))
with(nd, lines(m1 ~ x, col = 4))
coef(m1)
(dx <- with(nd, unique(diff(x)))[1])
with(nd, diff(m1)/diff(x+dx/2))
unique(with(nd, round(diff(m1)/diff(x+dx/2), dig=10)))


plot(y ~ x, data = dat04, ylim=c(0, 55))
with(nd, lines(m9 ~ x, col = 4))
coef(m9)
c(coef(m9)[2], sum(coef(m9)[2:3]))
with(nd, diff(m9)/diff(x+dx/2))
unique(with(nd, round(diff(m9)/diff(x+dx/2), dig=10)))


plot(y ~ x, data = dat04, ylim=c(0, 55))
with(nd, lines(m10 ~ x, col = 4))
coef(m10)
c(coef(m10)[2], sum(coef(m10)[2:3]))
with(nd, diff(m10)/diff(x+dx/2))
unique(with(nd, round(diff(m10)/diff(x+dx/2), dig=10)))

plot(y ~ x, data = dat04, ylim=c(0, 55))
with(nd, points(m10 ~ x, pch=19, cex=0.25))

(dx <- with(nd, unique(diff(x)))[1])
nd <- transform(nd, m10d1 = c(NA,diff(m10)/diff(x+dx/2)))

plot(m10d1 ~ x, data = nd, type = "n")
lines(m10d1 ~ x, data = subset(nd, x <= xc.est), type = "l")
lines(m10d1 ~ x, data = subset(nd, x > xc.est+dx), type = "l")

##
## 2 pontos de corte (3 segmentos)
##
xc1 <- 1.4
dat04 <- transform(dat04,
                   ind1 = as.numeric(x > xc1))
nd <- transform(nd,
                ind1 = as.numeric(x >= xc1))

## retas "livres"
(m11 <- lm(y ~ x + ind1 + I(ind1*x) + ind + I(ind*x), data = dat04))
nd$m11 <- predict(m11, newdata = nd)
plot(y ~ x, data = dat04, ylim=c(0, 55))
with(subset(nd, x < xc1), lines(m11 ~ x, col = 4))
with(subset(nd, x > xc1 & x < xc ), lines(m11 ~ x, col = 4))
with(subset(nd, x > xc), lines(m11 ~ x, col = 4))


## 3 retas "conectadas"
(m12 <- lm(y ~ x + I(ind1*(x-xc1)) + I(ind*(x-xc)), data = dat04))
nd$m12 <- predict(m12, newdata = nd)
plot(y ~ x, data = dat04, ylim=c(0, 55))
with(nd, lines(m12 ~ x, col = 4))


## 1as Derivadas
plot(y ~ x, data = dat04, ylim=c(0, 55))
with(nd, lines(m12 ~ x, col = 4))

coef(m12)
c(coef(m12)[2], sum(coef(m12)[2:3]), sum(coef(m12)[2:4]))
#with(nd, diff(m12)/diff(x+dx/2))
unique(with(nd, round(diff(m12)/diff(x+dx/2), dig=10)))

nd <- transform(nd, m12d1 = c(NA,diff(m12)/diff(x)))

plot(m12d1 ~ x, data = nd, type = "n")
lines(m12d1 ~ x, data = subset(nd, x <= xc1), type = "l")
lines(m12d1 ~ x, data = subset(nd, x > xc1+dx & x <= xc), type = "l")
lines(m12d1 ~ x, data = subset(nd, x > xc+dx), type = "l")


## Funções "base" I(ind1*(x-xc1)) e I(ind*(x-xc))
with(nd, matplot(x, cbind(I(ind1*(x-xc1)),I(ind*(x-xc))), type = "l"))