library(np)
library(locfit)
### Donnnées #######

## Data A
data(cps71)
head(cps71)


## Data B
load('Alcool.Rdata')
head(Alcool)

## Data C
load('Hormone.Rdata')
head(Hormone)


##### Section I  ################

## I-1-a) 
data(cps71)
head(cps71)

## I-1-b)
X <- cps71[,2]
Y <- cps71[,1]
plot(X,Y)

## I-1-c)
x0 <- seq(min(X), max(X), l=1000)

## I-2) 
mod.kern <- npreg(Y~X, exdat=x0)
rkern <- mod.kern$mean
plot(X,Y)
lines(x0, rkern, col='red')

## I-3) 
mod.lp <- locfit(Y~lp(X))
rlp <- predict(mod.lp, x0)
plot(X,Y)
lines(x0, rlp, col='red')

## I-4)
plot(X,Y)
lines(lowess(X,Y), col='red')
 lines(  lowess(X,Y, f=0.1), col='blue')



##### Section II  ################

## II-1-a)
load('Alcool.Rdata')
X <- Alcool$Alcohol.consumption
Y <- Alcool$Life
n = length(X)

## II-1-b)
x0 <- seq(min(X), max(X), l=1000)



## II-2-a) 
mod.kern <- npreg(Y~X, exdat=x0)
rkern <- mod.kern$mean
mod.lp <- locfit(Y~lp(X))
rlp <- predict(mod.lp, x0)


plot(X,Y)
lines(x0, rkern, col='red')
lines(x0, rlp, col='blue')

## II-3-a) 
J0 <- rep(1:10, l=n)
J1 <- sample(J0)

## II-3-b)
mod <- locfit(Y[J1!=1]~lp(X[J1!=1]))
Yhat <- predict(mod, X[which(J1==1)])
plot(Y[which(J1==1)],Yhat)

## alternative plus détaillée

## ech d'apprentissage
Xtrain <- X[which(J1!=1)]
Ytrain <- Y[which(J1!=1)]
## ech de validation
Xvalid <- X[which(J1==1)]
Yvalid <- Y[which(J1==1)]
## calcul de l'estimateur a partir de l'ech d'apprentissage
mod <- locfit(Ytrain~lp(Xtrain))
Yhat <- predict(mod, Xvalid)
plot(Yvalid, Yhat)


## II-3-c)
Yhat.lp <- vector(l=n)
for(l in 1:10){
	mod <- locfit(Y[J1!=l]~lp(X[J1!=l]))
	Yhat <- predict(mod, X[which(J1==l)])
	Yhat.lp[which(J1==l)] <- Yhat
}

## II-3-d)
plot(Y, Yhat.lp)
abline(0,1)
cor(Y, Yhat.lp)

## II-3-e)
Ecv.lp <- mean((Y-Yhat.lp)^2)


II-4)
r<- function(x) 10^(-4)*(x+40)^3*exp(-(x+40)/20)
plot(X,Y)
lines(x0,rlp, col='red') ## rlp calculé en II-2-b
lines(x0, r(x0), col='blue')

#########   Section III ###########

## III-1)
load('Hormone.Rdata')
X <- Hormone$Dose
Y <- Hormone$Endocrine
n <- length(X)
plot(X,Y)
x0 <- seq(min(X), max(X), l=1000)

## III-2)
mod.kern <- npreg(Y~X, exdat=x0)
rkern <- mod.kern$mean
mod.lp <- locfit(Y~lp(X))
rlp <- predict(mod.lp, x0)



plot(X,Y)
lines(x0, rkern, col='red')
lines(x0, rlp, col='blue')


## III-3) 
J0 <- rep(1:10, l=n)
J1 <- sample(J0)


Yhat.lp <- vector(l=n)
for(l in 1:n){
	mod <- locfit(Y[J1!=l]~lp(X[J1!=l]))
	Yhat.lp[which(J1==l)] <- predict(mod, X[which(J1==l)])
}

## II-3-d)
plot(Y, Yhat.lp)
abline(0,1)
cor(Y, Yhat.lp)

## III-4)

Yhat.lm <- vector(l=n)
for(l in 1:n){
	mod <- lm(Y[J1!=l]~lp(X[J1!=l]))
	co <- coef(mod)
	Yhat.lm[which(J1==l)] <- co[1] + co[2]*X[which(J1==l)]
}
plot(Y, Yhat.lm)
abline(0,1)
cor(Y, Yhat.lm)

Ecv.lp <- mean((Y-Yhat.lp)^2)
Ecv.lm <- mean((Y-Yhat.lm)^2)