#Dessin ortho de Gram Schmidt
GSplot<-function(){
library(rgl)
#Matrice des vecteurs de depart
M<-cbind(x=c(1,0,1,0,0.5),y=c(0,0,1,0,0.5),z=c(0,0,0,0,1))

plot3d(M,type="l",xlim=c(-1,1),ylim=c(-1,1),
zlim=c(0,1),col="dark grey")
arrow<-cbind(x=c(0.9,0.9,1,1-0.05/sqrt(2),1-.15/sqrt(2),1,0.45,0.5,0.5),
             y=c(-0.05,0.05,0,1-0.15/sqrt(2),1-0.05/sqrt(2),1,0.5,0.45,0.5),
             z=c(0,0,0,0,0,0,0.95,0.95,1))
triangles3d(arrow,col="dark grey")
readline()


Mchap<-cbind(x=c(1,1,0.5,0.5,0.5,0),y=c(1,0,0.5,0,0.5,0.5),z=c(0,0,1,0,1,0))
segments3d(Mchap[1:2,],col="blue")
readline()

pi0<-cbind(x=c(0,1),y=c(0.01,0.01),z=c(0,0))
segments3d(pi0)
arrow<-cbind(x=c(0.9,0.9,1),y=c(-0.06,0.04,0),z=c(0,0,0))
triangles3d(arrow)
text3d(x=0.95,
	 y=-0.1,
       z=0,"p0")
readline()



Minvpi0<-cbind(x=c(0,-1),y=c(0,0),z=c(0,0))
segments3d(Minvpi0,col="red")
arrow<-cbind(x=c(-1,-0.9,-0.9),y=c(0,-0.05,0.05),z=c(0,0,0))
triangles3d(arrow,col="red")
readline()

Adnegpi0vec1<-cbind(x=c(-1,0,0,1),
                      y=c(0,1,1,1),
                      z=c(0,0,0,0))
segments3d(Adnegpi0vec1,col="red")
readline()

pi1<-cbind(x=c(0,0),y=c(0,1),z=c(0,0))
segments3d(pi1)
arrow<-cbind(x=c(-0.05,0.05,0),y=c(0.95,0.95,1),z=c(0,0,0))
triangles3d(arrow)
text3d(x=-0.05,
	 y=0.95,
       z=0,"p1")
readline()

segments3d(Mchap[3:4,],col="blue")
segments3d(Mchap[5:6,],col="blue")
readline()


Minvpi1<-cbind(x=c(0,0),y=c(0,-0.5),z=c(0,0))
segments3d(Minvpi1,col="red")
Adnegpi0negpi1<-cbind(x=c(0,0,-0.5,-0.5,-0.5,0,0,-0.5),
                      y=c(0,-0.5,0,-0.5,-0.5,-0.5,0,-0.5),
                      z=c(0,0,0,0,0,0,0,0))
segments3d(Adnegpi0negpi1,col="red")
arrow<-cbind(x=c(0,-0.05,0.05,-0.5,-0.4,-0.4),
             y=c(-0.5,-0.4,-0.4,0,-0.05,0.05),
             z=c(0,0,0,0,0,0))
triangles3d(arrow,col="red")
readline()



pi2<-cbind(x=c(0,0),y=c(0,0),z=c(0,1))
ad<-cbind(x=c(-0.5,0,0,0.5),y=c(-0.5,0,0,0.5),z=c(0,1,1,1))
segments3d(ad,col="red")
segments3d(pi2)
arrow<-cbind(x=c(0,0,0),y=c(-0.05,0.05,0),z=c(0.9,0.9,1))
triangles3d(arrow)
text3d(x=-0.05,
	 y=0.05,
       z=0.95,"p2")
readline()

}

regpol<-function(){

###########################################################################"
#Exemple de régression polynomiale
#y=2+2x+3*x^2
nind<-30
x<-runif(nind,0,10)
y<-2+2*x+3*x^2+rnorm(nind,0,5)

####CALCUL PARAM
X<-cbind(rep(1,nind),x,x^2)

XtX<-t(X)%*%X
XtXinv<-solve(XtX)
achap<-XtXinv%*%t(X)%*%y

###PERMET DE REPRESENTER UNE GRILLE
### CREATION MATRICE en COL J et LIG I CORRESPOND Z
### GRACE A outer
x1<-seq(0,10,length=100)
x2<-x1^2
reg<-function(x1,x2,achap){
	res<-achap[1]+achap[2]*x1+achap[3]*x2
	return(res)
}
z<-outer(x1,x2,reg,achap)

###REPRESENTATION GRILLE + DONNEES + EQM
persp3d(x1,x2,z,
xlab="x1",ylab="x2",zlab="y",
col="red", alpha=0.7, aspect=c(1,1,0.5),box=F)
grid3d(c("x", "y+", "z"))
points3d(x,x^2,y,size=5)

#plot3d(x,x^2,y,type="s",radius=1,col="red",box=F,
#xlim=c(min(x),max(x)),
#ylim=c(min(x^2),max(x^2)),zlim=c(min(y),max(y)))
xabs<-seq(0,10,0.1)
ymod<-achap[1]+achap[2]*xabs+achap[3]*xabs^2
lines3d(xabs,xabs^2,ymod,col=2)

##REPRESENTATION DES DONNEES DANS PLAN XY (X1,Z en fait)
lines3d(xabs,0,ymod,col="blue")
points3d(x,0,y,radius=1,col="blue",size=5)

###TRACE DES EQM
for(i in 1:nind){
	yvrai<-c(x[i],x[i]^2,y[i])
	yest<-c(x[i],x[i]^2,reg(x[i],x[i]^2,achap))
	segments3d(rbind(yvrai,yest))
	yvrai<-c(x[i],0,y[i])
	yest<-c(x[i],0,reg(x[i],x[i]^2,achap))
	segments3d(rbind(yvrai,yest),col="blue")
}

}

regpol2<-function(){
#EXEMPLE DE REGRESSION A 2 VARIABLES

x1<-runif(50,0,20)
x2<-runif(50,0,20)
y<-2+3*x1-4*x2+rnorm(50,0,10)

plot3d(x1,x2,y,type="s",radius=1)


X<-cbind(rep(1,50),x1,x2)

XtX<-t(X)%*%X
XtXinv<-solve(XtX)
achap<-XtXinv%*%t(X)%*%y

reg<-function(x1,x2,achap){
	res<-achap[1]+achap[2]*x1+achap[3]*x2
	return(res)
}
z<-outer(1:20,1:20,reg,achap)

persp3d(1:20,1:20,z,
xlab="x1",ylab="x2",zlab="y",
col="red", alpha=0.7, aspect=c(1,1,0.5),box=F)
grid3d(c("x", "y+", "z"))
points3d(x1,x2,y,size=5)

for(i in 1:50){
	yvrai<-c(x1[i],x2[i],y[i])
	yest<-c(x1[i],x2[i],reg(x1[i],x2[i],achap))
	segments3d(rbind(yvrai,yest))
}


}

#REGRESSION 3 VARIABLES
reg3V<-function(){
x1<-runif(100,0,20)
x2<-runif(100,0,20)
y<-2+3*x1-4*x2+2*x1*x2+rnorm(50,0,100)

plot3d(x1,x2,y,type="s",radius=5)


X<-cbind(rep(1,100),x1,x2,x1*x2)

XtX<-t(X)%*%X
XtXinv<-solve(XtX)
achap<-XtXinv%*%t(X)%*%y

reg<-function(x1,x2,achap){
	res<-achap[1]+achap[2]*x1+achap[3]*x2+achap[4]*x1*x2
	return(res)
}
z<-outer(0:20,0:20,reg,achap)

persp3d(0:20,0:20,z,
xlab="x1",ylab="x2",zlab="y",
col="red", alpha=0.7, aspect=c(1,1,0.5),box=F)
grid3d(c("x", "y+", "z"))
points3d(x1,x2,y,size=5)

for(i in 1:100){
	yvrai<-c(x1[i],x2[i],y[i])
	yest<-c(x1[i],x2[i],reg(x1[i],x2[i],achap))
	segments3d(rbind(yvrai,yest))
}
rgl.viewpoint(0,-60,zoom=2/3)

cat("appuyez sur une touche", fill=T)
rep<-scan("",character(),1)

start <- proc.time()[3]
while ((i <- 36*(proc.time()[3]-start)) < 360) {
rgl.viewpoint(i,-60,zoom=2/3); 
}
}