library(MASS)
library(stats)
source('formule.R')
source('mat_modele.R')
source("triangulaire.R")

# étude des plans : graphique + optimalité
P1 <- read.table('Unif_15_1.dat',header = T)
P2 <- read.table('Unif_15_2.dat',header = T)
P3 <- read.table('Unif_15_3.dat',header = T)
P4 <- read.table('Unif_15_4.dat',header = T)

det <- NULL

P1 <- mat_modele('quadratique',P1)
P2 <- mat_modele('quadratique',P2)
P3 <- mat_modele('quadratique',P3)
P4 <- mat_modele('quadratique',P4)

det <- c(det,det(t(P1)%*%P1))
det <- c(det,det(t(P2)%*%P2))
det <- c(det,det(t(P3)%*%P3))
det <- c(det,det(t(P4)%*%P4))

barplot(det,names.arg = c('P1','P2','P3','P4'),main ='det(t(X)X)')

A_opt <- NULL
A_opt <- c(A_opt,sum(diag(solve(t(P1)%*%P1,diag(rep(1,10))))))
A_opt <- c(A_opt,sum(diag(solve(t(P2)%*%P2,diag(rep(1,10))))))
A_opt <- c(A_opt,sum(diag(solve(t(P3)%*%P3,diag(rep(1,10))))))
A_opt <- c(A_opt,sum(diag(solve(t(P4)%*%P4,diag(rep(1,10))))))

barplot(A_opt,names.arg = c('P1','P2','P3','P4'),main ='Tr(inv(t(X)X))')


#calcul de la surface de réponse
ff <- formule('quadratique')
mod1 <- lm(ff,read.table('Unif_15_2.dat',header = T))
summary(mod1)
res <- studres(mod1)
plot(res~mod1$fitted.values)


##propagation des incertitudes`

#utilisation de la fonction predict
#cas de tirages uniformes des facteurs
PORO <- runif(1000,-1,1)
MK3 <- runif(1000,-1,1)
KR <- runif(1000,-1,1)

F <- data.frame(PORO,MK3,KR)
CP1 <- predict.lm(mod1,F,interval = "none",level = 0.95, type = "response")
CP2 <- predict.lm(mod1,F,interval = "confidence",level = 0.95, type = "response")
V2 <- ((CP2[,3]-CP1)/1.96)^2 #obtention des variances à l'aide de predict
CP2 <-mvrnorm(1,CP1,diag(V2))
CP3 <- predict.lm(mod1,F,interval = "prediction",level = 0.95, type = "response")
V3 <- ((CP3[,3]-CP1)/1.96)^2
CP3 <-mvrnorm(1,CP1,diag(V3))

summary(CP1)
summary(CP2)
summary(CP3)

layout(matrix(1:3,3,1))
plot(density(CP1),main=c("propagation des incertitudes","modèle linéaire ajusté"),xlim=c(3500000,6500000))
plot(density(CP2),main=c("propagation des incertitudes","modèle linéaire avec intervalle de confiance"),xlim=c(3500000,6500000))
plot(density(CP3),main=c("propagation des incertitudes","modèle linéaire avec intervalle de prévision"),xlim=c(3500000,6500000))


plot(ecdf(CP1),main=c("propagation des incertitudes","modèle linéaire ajusté"),xlim=c(3500000,6500000),xaxp=c(3500000,6500000,12))
abline(a=0.1,b=0)
abline(a=0.9,b=0)
plot(ecdf(CP2),main=c("propagation des incertitudes","modèle linéaire avec intervalle de confiance"),xlim=c(3500000,6500000),xaxp=c(3500000,6500000,12))
abline(a=0.1,b=0)
abline(a=0.9,b=0)
plot(ecdf(CP3),main=c("propagation des incertitudes","modèle linéaire avec intervalle de prévision"),xlim=c(3500000,6500000),xaxp=c(3500000,6500000,12))
abline(a=0.1,b=0)
abline(a=0.9,b=0)

#utilisation de predict
#tirages des facteurs selon leur loi
PORO <- triangulaire(1000)
MK3 <- triangulaire(1000)
KR <- runif(1000,-1,1)

F <- data.frame(PORO,MK3,KR)
CP1 <- predict.lm(mod1,F,interval = "none",level = 0.95, type = "response")
CP2 <- predict.lm(mod1,F,interval = "confidence",level = 0.95, type = "response")
V2 <- ((CP2[,3]-CP1)/1.96)^2
CP2 <-mvrnorm(1,CP1,diag(V2))
CP3 <- predict.lm(mod1,F,interval = "prediction",level = 0.95, type = "response")
V3 <- ((CP3[,3]-CP1)/1.96)^2
CP3 <-mvrnorm(1,CP1,diag(V3))

summary(CP1)
summary(CP2)
summary(CP3)

layout(matrix(1:3,3,1))
plot(density(CP1),main=c("propagation des incertitudes","modèle linéaire ajusté"),xlim=c(3500000,6000000))
plot(density(CP2),main=c("propagation des incertitudes","modèle linéaire avec intervalle de confiance"),xlim=c(3500000,6000000))
plot(density(CP3),main=c("propagation des incertitudes","modèle linéaire avec intervalle de prévision"),xlim=c(3500000,6000000))


plot(ecdf(CP1),main=c("propagation des incertitudes","modèle linéaire ajusté"),xlim=c(3500000,6000000),xaxp=c(3500000,6000000,10))
abline(a=0.1,b=0)
abline(a=0.9,b=0)
plot(ecdf(CP2),main=c("propagation des incertitudes","modèle linéaire avec intervalle de confiance"),xlim=c(3500000,6000000),xaxp=c(3500000,6000000,10))
abline(a=0.1,b=0)
abline(a=0.9,b=0)
plot(ecdf(CP3),main=c("propagation des incertitudes","modèle linéaire avec intervalle de prévision"),xlim=c(3500000,6000000),xaxp=c(3500000,6000000,10))
abline(a=0.1,b=0)
abline(a=0.9,b=0)

#propagation d'incertitudes 
#à l'aide de simulations de Monte Carlo
mse<-var(mod1$residuals)*(mod1$rank+mod1$df-1)/mod1$df
cov_beta<-vcov(mod1)
beta<-mod1$coefficients

PORO <- triangulaire(1000)
MK3 <- triangulaire(1000)
KR <- runif(1000,-1,1)
CP<-rep(1,1000) # variable inutile, définie uniquement pour utiliser la fonction mat_modele
F <- data.frame(PORO,MK3,KR,CP)
FF <- mat_modele('quadratique',F)


#propagation de F seul
Y_hat<- FF%*%beta

#propagation avec intervalle de confiance
beta_simu<-mvrnorm(1000,mu=beta,Sigma=cov_beta)
Y_conf <- rep(0,1000)
for (i in 1:1000)
{
Y_conf[i] <- FF[i,]%*%beta_simu[i,]
}

#propagation avec intervalle de prévision
Y_prev <- rep(0,1000)
beta_simu <- mvrnorm(1000,mu=beta,Sigma=cov_beta)
eps_simu <- rnorm(1000,0,sqrt(mse))
for (i in 1:1000)
{
Y_prev[i] <- FF[i,]%*%beta_simu[i,]+ eps_simu[i]
}

#intervalle de prévision et simulation des résidus par bootstrap
Y_prev_boot <- rep(0,1000)
beta_simu <- mvrnorm(1000,mu=beta,Sigma=cov_beta)
eps<-sample(mod1$residuals,1000,replace=T)
for (i in 1:1000)
{
Y_prev_boot[i] <- FF[i,]%*%beta_simu[i,]+eps[i]
}


layout(matrix(1:3,3,1))
plot(density(Y_hat),main=c("propagation des incertitudes","modèle linéaire ajusté"),xlim=c(3500000,6000000))
plot(density(Y_conf),main=c("propagation des incertitudes","modèle linéaire avec intervalle de confiance"),xlim=c(3500000,6000000))
plot(density(Y_prev),main=c("propagation des incertitudes","modèle linéaire avec intervalle de prévision"),xlim=c(3500000,6000000))

layout(matrix(1:2,2,1))
plot(density(Y_prev),main=c("propagation des incertitudes","modèle linéaire avec intervalle de prévision"),xlim=c(3500000,6000000))
plot(density(Y_prev_boot),main=c("propagation des incertitudes","modèle linéaire avec intervalle de prévision par bootstrap"),xlim=c(3500000,6000000))


#calcul de la surface de réponse avec un D_opt
ff <- formule('quadratique')
mod2 <- lm(ff,read.table('15pts_3.dat',header = T))
summary(mod2)
res <- studres(mod2)
plot(res~mod1$fitted.values)


CP4 <- predict.lm(mod2,F,interval = "prediction",level = 0.95, type = "response")
V5 <- ((CP4[,3]-CP1)/1.96)^2
CP4 <-mvrnorm(1,CP1,diag(V5))


layout(matrix(1:2,2,1))
plot(density(CP3),main=c("propagation des incertitudes","plan uniforme"),xlim=c(3500000,6000000))
plot(density(CP4),main=c("propagation des incertitudes","plan D-optimal"),xlim=c(3500000,6000000))