#####################################
#regression - option maths applis
#EMSE
#L. Carraro (Telecom Saint-Etienne)
#decembre 09
#####################################


###################################################################################
#Première partie : obtention d'un échantillon de points en 2D formant deux familles
###################################################################################

classimu<-function(nsimu,sigma=1/5)
#classimu simule deux familles de centres de gravité, puis des échantillons mélanges de gaussiennes centrées sur ces points.
#pour une simulation avec centres de gravité fixés, voir classimu2
#nsimu est le nombre d'individus pour chaque échantillon
#sigma est l'écart-type des variables simulées
{
#on simule le premier échantillon
#les moyennes des gaussiennes
m1x<-rnorm(10,mean=1)
m1y<-rnorm(10,mean=0)
ind<-sample(1:10,nsimu,replace=T)
#l'échantillon
input1<-matrix(nrow=nsimu,ncol=2)
for (k in 1:nsimu)
	{
		input1[k,1]<-rnorm(n=1,mean=m1x[ind[k]],sd=sigma)
		input1[k,2]<-rnorm(n=1,mean=m1y[ind[k]],sd=sigma)
	}
input1<-cbind(input1,rep(1,times=nsimu))

#on simule le second échantillon
#les moyennes des gaussiennes
m2x<-rnorm(10,mean=0)
m2y<-rnorm(10,mean=1)
ind<-sample(1:10,nsimu,replace=T)
#l'échantillon
input2<-matrix(nrow=nsimu,ncol=2)
for (k in 1:nsimu)
	{
		input2[k,1]<-rnorm(n=1,mean=m2x[ind[k]],sd=sigma)
		input2[k,2]<-rnorm(n=1,mean=m2y[ind[k]],sd=sigma)
	}
input2<-cbind(input2,rep(0,times=nsimu))
#on reporte les résultats dans une liste
classimu<-list(ech=rbind(input1,input2),m1=rbind(m1x,m1y),m2=rbind(m2x,m2y))
}

#exemple
#x<-classimu(100,sigma=1/5)
#xech<-x$ech
#plot(xech[xech[,3]==1,],col="red",xlim=c(min(xech[,1]),max(xech[,1])),ylim=c(min(xech[,2]),max(xech[,2])))
#points(xech[xech[,3]==0,],col="blue")
#points(t(x$m1),pch=23,col="red",bg="red")
#points(t(x$m2),pch=23,col="blue",bg="blue")

classimu2<-function(nsimu,sigma=1/5,m1,m2)
#classimu2 simule des échantillons mélanges de gaussiennes centrées sur les points donnés dans m1, puis dans m2.
#pour une simulation avec centres de gravité tirés au hasard, voir classimu
#nsimu est le nombre d'individus pour chaque échantillon
#sigma est l'écart-type des variables simulées
#m1 et m2 donnent des coordonnées des centres de gravités des deux familles et ont chacun 2 lignes
{
#on simule le premier échantillon
n1<-length(m1[1,])
ind<-sample(1:n1,nsimu,replace=T)
#l'échantillon
input1<-matrix(nrow=nsimu,ncol=2)
for (k in 1:nsimu)
	{
		input1[k,1]<-rnorm(n=1,mean=m1[1,ind[k]],sd=sigma)
		input1[k,2]<-rnorm(n=1,mean=m1[2,ind[k]],sd=sigma)
	}
input1<-cbind(input1,rep(1,times=nsimu))

#on simule le second échantillon
#les moyennes des gaussiennes

n2<-length(m2[1,])
ind<-sample(1:n2,nsimu,replace=T)
#l'échantillon
input2<-matrix(nrow=nsimu,ncol=2)
for (k in 1:nsimu)
	{
		input2[k,1]<-rnorm(n=1,mean=m2[1,ind[k]],sd=sigma)
		input2[k,2]<-rnorm(n=1,mean=m2[2,ind[k]],sd=sigma)
	}
input2<-cbind(input2,rep(0,times=nsimu))
#on reporte les résultats dans une liste
classimu2<-list(ech=rbind(input1,input2),m1=m1,m2=m2,sigma=sigma)
}

#exemple (x est obtenue par classimu : voir exemple précédent)
data.class<-classimu2(100,sigma=0.25,m1=x$m1,m2=x$m2)
ech<-data.class$ech
plot(ech[ech[,3]==1,],col="red",xlim=c(min(ech[,1]),max(ech[,1])),ylim=c(min(ech[,2]),max(ech[,2])))
points(ech[ech[,3]==0,],col="blue")
points(t(data.class$m1),pch=23,col="red",bg="red")
points(t(data.class$m2),pch=23,col="blue",bg="blue")

#une telle simulation est sauvegardée dans le fichier classification.Rdata
save(data.class,file="classification.Rdata")

#############################################################################################
#Deuxième partie : représentation des densités des points simulés et de la frontière optimale
#############################################################################################

load(file="classification.Rdata")
#on représente les centres de gravité
plot(t(data.class$m1),pch=23,col="red",bg="red",xlim=c(min(data.class$ech[,1]),max(data.class$ech[,1])),ylim=c(min(data.class$ech[,2]),max(data.class$ech[,2])),,xlab="x",ylab="y")
points(t(data.class$m2),pch=23,col="blue",bg="blue")
#puis on ajoute les points simulés
ech<-data.class$ech
points(ech[ech[,3]==1,],col="red")
points(ech[ech[,3]==0,],col="blue")

#densités des points simulés
density.class<-function(x,y,sigma,m)
#(x,y) est le point courant
#sigma l'écart-type des gaussiennes (matrice de covariance : sigma^2*Id)
#m la matrice des espérances des gaussiennes
{
density.class<-mean(dnorm(x,mean=m[1,],sd=sigma)*dnorm(y,mean=m[2,],sd=sigma))
}
n=200
x<-seq(-3,4,length.out=n)
y<-seq(-2,3,length.out=n)
plan<-expand.grid(x=x,y=y)
z1<-matrix(nrow=n,ncol=n)
z2<-matrix(nrow=n,ncol=n)
for (i in 1:n)
	{
	for (j in 1:n)
		{
		z1[i,j]<-density.class(x[i],y[j],data.class$sigma,data.class$m1)
		z2[i,j]<-density.class(x[i],y[j],data.class$sigma,data.class$m2)
		}
	}


#on représente les densités des deux populations
ech<-data.class$ech
plot(ech[ech[,3]==1,],col="red",xlim=c(min(ech[,1]),max(ech[,1])),ylim=c(min(ech[,2]),max(ech[,2])),xlab="x",ylab="y")
title(main="densités des deux populations")
points(ech[ech[,3]==0,],col="blue")
points(t(data.class$m1),pch=23,col="red",bg="red")
points(t(data.class$m2),pch=23,col="blue",bg="blue")
contour(x,y,z1,col="red",add=T)
contour(x,y,z2,col="blue",add=T)
indic<-(z1>z2)
contour(x,y,indic,add=T)
legend("topright",cex=0.8,legend=c("densité population 1","densité population 2","frontière optimale"),lty=1,lwd=c(1,1,3),col=c("red","blue","black"))

#représentation 3d interactive
library(rgl)
persp3d(x,y,z1,col="red")
xy<-expand.grid(x,y)
lines3d(xy[,1],xy[,2],z2,col="blue")

#tableau des bons et mauvais classements
#classements optimaux
n.ech<-length(data.class$ech[,1])
dens<-matrix(nrow=n.ech,ncol=2)
class.opt<-vector(length=n.ech)
for (i in 1:n.ech)
	{
	dens[i,1]<-density.class(data.class$ech[i,1],data.class$ech[i,2],data.class$sigma,data.class$m1)
	dens[i,2]<-density.class(data.class$ech[i,1],data.class$ech[i,2],data.class$sigma,data.class$m2)
	class.opt[i]<-(dens[i,1]>dens[i,2])		#TRUE si classe = 1 = red
	}
#comparaison aux véritables classes
error.opt<-(class.opt==data.class$ech[,3])
#% de bien classés
prop.opt<-length(error.opt[error.opt==T])/length(error.opt)

#############################################################################################
#Troisième partie : estimation d'un modèle linéaire et séparation résultante
#############################################################################################
load(file="classification.Rdata")
data<-as.data.frame(data.class$ech)
names(data)<-c("x","y","classe")
modlin<-lm(classe~x+y,data=data)

#classement obtenu avec un modèle linéaire de degré 1
#####################################################
n.ech<-length(data.class$ech[,1])
class.lin<-vector(length=n.ech)
for (i in 1:n.ech)
	{
	class.lin[i]<-(modlin$fitted.values[i]>0.5)		#TRUE si classe = 1 = red
	}
#comparaison aux véritables classes
error.lin<-(class.lin==data.class$ech[,3])
#% de bien classés
prop.lin<-length(error.lin[error.lin==T])/length(error.lin)

#représentation graphique
ech<-data.class$ech
plot(ech[ech[,3]==1,],col="red",xlim=c(min(ech[,1]),max(ech[,1])),ylim=c(min(ech[,2]),max(ech[,2])),xlab="x",ylab="y")
title(main="séparation linéaire des deux classes")
points(ech[ech[,3]==0,],col="blue")
points(t(data.class$m1),pch=23,col="red",bg="red")
points(t(data.class$m2),pch=23,col="blue",bg="blue")
n=200
x<-seq(-3,4,length.out=n)
y<-seq(-2,3,length.out=n)
plan<-expand.grid(x=x,y=y)
pred.lin<-predict(modlin,new=data.frame(x=plan[,1],y=plan[,2]))
pred.lin<-matrix(pred.lin,nrow=n,ncol=n)
indic1<-(pred.lin>0.5)
contour(x,y,indic1,add=T,lwd=2,nlevels=2)
contour(x,y,indic,add=T,lty=3,lwd=1,col="green",nlevels=2)
legend("topright",cex=0.8,legend=c("population 1","population 2","frontière estimée","frontière de Bayes"),lty=c(1,1,1,3),lwd=c(1,1,3,1),col=c("red","blue","black","green"))

#modèle linéaire de degré 2
###########################
modlin2<-lm(classe~I(x^2)+I(y^2)+(x+y)^2,data=data)

#classement obtenu avec un modèle linéaire
n.ech<-length(data.class$ech[,1])
class.lin2<-vector(length=n.ech)
for (i in 1:n.ech)
	{
	class.lin2[i]<-(modlin2$fitted.values[i]>0.5)     #TRUE si classe = 1 = red
	}
#comparaison aux véritables classes
error.lin2<-(class.lin2==data.class$ech[,3])
#% de bien classés
prop.lin2<-length(error.lin2[error.lin2==T])/length(error.lin2)

#représentation graphique
ech<-data.class$ech
plot(ech[ech[,3]==1,],col="red",xlim=c(min(ech[,1]),max(ech[,1])),ylim=c(min(ech[,2]),max(ech[,2])),xlab="x",ylab="y")
title(main=c("séparation des deux classes","modèle linéaire de degré 2"))
points(ech[ech[,3]==0,],col="blue")
points(t(data.class$m1),pch=23,col="red",bg="red")
points(t(data.class$m2),pch=23,col="blue",bg="blue")
n=200
x<-seq(-3,4,length.out=n)
y<-seq(-2,3,length.out=n)
plan<-expand.grid(x=x,y=y)
pred.lin2<-predict(modlin2,new=data.frame(x=plan[,1],y=plan[,2]))
pred.lin2<-matrix(pred.lin2,nrow=n,ncol=n)
indic2<-(pred.lin2>0.5)
contour(x,y,indic2,add=T,lwd=2,nlevels=2)
contour(x,y,indic,add=T,lty=3,lwd=1,col="green",nlevels=2)
legend("topright",cex=0.8,legend=c("population 1","population 2","frontière estimée","frontière de Bayes"),lty=c(1,1,1,3),lwd=c(1,1,3,1),col=c("red","blue","black","green"))

#modèle linéaire de degré 5
###########################
modlin3<-lm(classe~(x+y)^2+I(x^2)+I(x^3)+I(x^4)+I(x^5)+I(y^2)+I(y^3)+I(y^4)+I(y^5)+I(x^2*y)+I(x^3*y)+I(x^4*y)+I(x^5*y)+I(x*y^2)+I(x*y^3)+I(x*y^4)+I(x*y^5)+I(x^2*y^2)+I(x^3*y^2)+I(x^4*y^2)+I(x^5*y^2)+I(x^2*y^3)+I(x^2*y^4)+I(x^2*y^5)+I(x^3*y^3)+I(x^4*y^3)+I(x^5*y^3)+I(x^3*y^4)+I(x^3*y^5)+I(x^4*y^4)+I(x^5*y^4)+I(x^4*y^5)+I(x^5*y^5),data=data)

#classement obtenu avec un modèle linéaire
n.ech<-length(data.class$ech[,1])
class.lin3<-vector(length=n.ech)
for (i in 1:n.ech)
	{
	class.lin3[i]<-(modlin3$fitted.values[i]>0.5)     #TRUE si classe = 1 = red
	}
#comparaison aux véritables classes
error.lin3<-(class.lin3==data.class$ech[,3])
#% de bien classés
prop.lin3<-length(error.lin3[error.lin3==T])/length(error.lin3)

#représentation graphique
ech<-data.class$ech
plot(ech[ech[,3]==1,],col="red",xlim=c(min(ech[,1]),max(ech[,1])),ylim=c(min(ech[,2]),max(ech[,2])),xlab="x",ylab="y")
title(main=c("séparation des deux classes","modèle linéaire de degré 5"))
points(ech[ech[,3]==0,],col="blue")
points(t(data.class$m1),pch=23,col="red",bg="red")
points(t(data.class$m2),pch=23,col="blue",bg="blue")
n=200
x<-seq(-3,4,length.out=n)
y<-seq(-2,3,length.out=n)
plan<-expand.grid(x=x,y=y)
pred.lin3<-predict(modlin3,new=data.frame(x=plan[,1],y=plan[,2]))
pred.lin3<-matrix(pred.lin3,nrow=n,ncol=n)
indic3<-(pred.lin3>0.5)
contour(x,y,indic3,add=T,lwd=2,nlevels=2)
contour(x,y,indic,add=T,lty=3,lwd=1,col="green",nlevels=2)
legend("topright",cex=0.8,legend=c("population 1","population 2","frontière estimée","frontière de Bayes"),lty=c(1,1,1,3),lwd=c(1,1,3,1),col=c("red","blue","black","green"))


#############################################################################################
#Quatrième partie : utilisation de kNN, i.e. des plus proches voisins (vote)
#############################################################################################
library(class)
n<-200
x<-seq(-3,4,length.out=n)
y<-seq(-2,3,length.out=n)
plan<-expand.grid(x=x,y=y)

#k=1
###############
class.knn<-knn(train=data.class$ech[,1:2],test=plan,cl=data.class$ech[,3],k=1)

ech<-data.class$ech
plot(ech[ech[,3]==1,],col="red",xlim=c(min(ech[,1]),max(ech[,1])),ylim=c(min(ech[,2]),max(ech[,2])),xlab="x",ylab="y")
title(main=c("séparation des deux classes","kNN avec k=1"))
points(ech[ech[,3]==0,],col="blue")
points(t(data.class$m1),pch=23,col="red",bg="red")
points(t(data.class$m2),pch=23,col="blue",bg="blue")
indic.knn<-matrix(class.knn,nrow=n)
contour(x,y,indic.knn,add=T,lwd=2,nlevels=2)
contour(x,y,indic,add=T,lty=3,lwd=1,col="green",nlevels=2)
legend("topright",cex=0.8,legend=c("frontière de Bayes","frontière estimée"),lty=c(3,1),lwd=c(1,2),col=c("green","black"))


class.knn<-knn(train=data.class$ech[,1:2],test=data.class$ech[,1:2],cl=data.class$ech[,3],k=1)
error.knn<-(class.knn==data.class$ech[,3])
prop.knn1<-length(error.knn[error.knn==T])/length(error.knn)

#k=10
###############
class.knn<-knn(train=data.class$ech[,1:2],test=plan,cl=data.class$ech[,3],k=10)

ech<-data.class$ech
plot(ech[ech[,3]==1,],col="red",xlim=c(min(ech[,1]),max(ech[,1])),ylim=c(min(ech[,2]),max(ech[,2])),xlab="x",ylab="y")
title(main=c("séparation des deux classes","knn avec k=10"))
points(ech[ech[,3]==0,],col="blue")
points(t(data.class$m1),pch=23,col="red",bg="red")
points(t(data.class$m2),pch=23,col="blue",bg="blue")
indic.knn<-matrix(class.knn,nrow=n)
contour(x,y,indic.knn,add=T,lwd=2,nlevels=2)
contour(x,y,indic,add=T,lty=3,lwd=1,col="green",nlevels=2)
legend("topright",cex=0.8,legend=c("frontière de Bayes","frontière estimée"),lty=c(3,1),lwd=c(1,2),col=c("green","black"))


class.knn<-knn(train=data.class$ech[,1:2],test=data.class$ech[,1:2],cl=data.class$ech[,3],k=10)
error.knn<-(class.knn==data.class$ech[,3])
prop.knn10<-length(error.knn[error.knn==T])/length(error.knn)

#k=30
###############
class.knn<-knn(train=data.class$ech[,1:2],test=plan,cl=data.class$ech[,3],k=30)

ech<-data.class$ech
plot(ech[ech[,3]==1,],col="red",xlim=c(min(ech[,1]),max(ech[,1])),ylim=c(min(ech[,2]),max(ech[,2])),xlab="x",ylab="y")
title(main=c("séparation des deux classes","knn avec k=30"))
points(ech[ech[,3]==0,],col="blue")
points(t(data.class$m1),pch=23,col="red",bg="red")
points(t(data.class$m2),pch=23,col="blue",bg="blue")
indic.knn<-matrix(class.knn,nrow=n)
contour(x,y,indic.knn,add=T,lwd=2,nlevels=2)
contour(x,y,indic,add=T,lty=3,lwd=1,col="green",nlevels=2)
legend("topright",cex=0.8,legend=c("frontière de Bayes","frontière estimée"),lty=c(3,1),lwd=c(1,2),col=c("green","black"))

class.knn<-knn(train=data.class$ech[,1:2],test=data.class$ech[,1:2],cl=data.class$ech[,3],k=30)
error.knn<-(class.knn==data.class$ech[,3])
prop.knn30<-length(error.knn[error.knn==T])/length(error.knn)

#############################################################################################
#Cinquième partie : quantile du minimum des rayons de vecteurs uniformes dans R^d
#############################################################################################
#On simule n points uniformes dans [-1,1]^d et R est le minimum des normes infinies de ces points
#quantile q_alpha = (1-(1-alpha)^n)^(1/d)

quant.d<-function(alpha,nb.points,dim)
{
result=(1-(1-alpha)^(1/nb.points))^(1/dim)
}

nb.points<-200
dim=1
x<-(1:1000)/1000
y<-quant.d(x,nb.points,dim)
d<-rep(dim,1000)
data<-cbind(d,y)
plot(x,y,type="l",lty=1,xlab=expression(paste("probabilité ",alpha)),ylab="quantile")
title(main=c("quantile du rayon du plus proche voisin en fonction de la dimension","échantillon de 200 points"))

dim=5
y<-quant.d(x,nb.points,dim)
d<-rep(dim,1000)
data<-rbind(data,cbind(d,y))
lines(x,y,col="blue")

dim=10
y<-quant.d(x,nb.points,dim)
d<-rep(dim,1000)
data<-rbind(data,cbind(d,y))
lines(x,y,col="green")

dim=15
y<-quant.d(x,nb.points,dim)
d<-rep(dim,1000)
data<-rbind(data,cbind(d,y))
lines(x,y,col="orange")

dim=20
y<-quant.d(x,nb.points,dim)
d<-rep(dim,1000)
data<-rbind(data,cbind(d,y))
lines(x,y,col="red")

dim=40
y<-quant.d(x,nb.points,dim)
d<-rep(dim,1000)
data<-rbind(data,cbind(d,y))
lines(x,y,col="cyan")

legend("topleft",lty=1,col=c("black","blue","green","orange","red","cyan"),legend=c("dim=1","dim=5","dim=10","dim=15","dim=20","dim=40"),cex=0.8)

#second graphique
data<-data.frame(data)
boxplot(y~d,data=data,xlab="dimension",ylab="quantiles")
title(main=c("boxplot de la loi du rayon du point","le plus proche de 0"))

#############################################################################################
#Sixième partie : retour sur l'exemple simulé et choix de modèle
#############################################################################################

#Construction de l'ensemble de validation
data.class.test<-classimu2(1000,sigma=data.class$sigma,m1=data.class$m1,m2=data.class$m2)
test<-data.class.test$ech
indic<-matrix((test[,3]==1),nrow=2000)

#modèle linéaire de degré 1
pred.lin<-predict(modlin,new=data.frame(x=test[,1],y=test[,2]))
pred.lin<-matrix(pred.lin,nrow=2000)
indic1<-(pred.lin>0.5)
error<-(indic==indic1)
prop.ext.lin1<-length(error[error==T])/length(error)
plot(test[test[,3]==1,1:2],col="red",xlim=c(-3,4),ylim=c(-2,3),xlab="x",ylab="y")
title(main=c("validation sur un échantillon test","modèle linéaire de degré 1"))
points(test[test[,3]==0,1:2],col="blue")
x<-seq(-3,4,length.out=n)
y<-seq(-2,3,length.out=n)
plan<-expand.grid(x=x,y=y)
pred.lin1<-predict(modlin,new=data.frame(x=plan[,1],y=plan[,2]))
pred.lin1<-matrix(pred.lin1,nrow=n,ncol=n)
indic1<-(pred.lin1>0.5)
contour(x,y,indic1,add=T,lwd=2,nlevels=2)

#modèle linéaire de degré 2
pred.lin<-predict(modlin2,new=data.frame(x=test[,1],y=test[,2]))
pred.lin<-matrix(pred.lin,nrow=2000)
indic2<-(pred.lin>0.5)
error<-(indic==indic2)
prop.ext.lin2<-length(error[error==T])/length(error)
plot(test[test[,3]==1,1:2],col="red",xlim=c(-3,4),ylim=c(-2,3),xlab="x",ylab="y")
title(main=c("validation sur un échantillon test","modèle linéaire de degré 2"))
points(test[test[,3]==0,1:2],col="blue")
x<-seq(-3,4,length.out=n)
y<-seq(-2,3,length.out=n)
plan<-expand.grid(x=x,y=y)
pred.lin2<-predict(modlin2,new=data.frame(x=plan[,1],y=plan[,2]))
pred.lin2<-matrix(pred.lin2,nrow=n,ncol=n)
indic2<-(pred.lin2>0.5)
contour(x,y,indic2,add=T,lwd=2,nlevels=2)

#modèle linéaire de degré 5
pred.lin<-predict(modlin3,new=data.frame(x=test[,1],y=test[,2]))
pred.lin<-matrix(pred.lin,nrow=2000)
indic3<-(pred.lin>0.5)
error<-(indic==indic3)
prop.ext.lin3<-length(error[error==T])/length(error)
plot(test[test[,3]==1,1:2],col="red",xlim=c(-3,4),ylim=c(-2,3),xlab="x",ylab="y")
title(main=c("validation sur un échantillon test","modèle linéaire de degré 5"))
points(test[test[,3]==0,1:2],col="blue")
x<-seq(-3,4,length.out=n)
y<-seq(-2,3,length.out=n)
plan<-expand.grid(x=x,y=y)
pred.lin3<-predict(modlin3,new=data.frame(x=plan[,1],y=plan[,2]))
pred.lin3<-matrix(pred.lin3,nrow=n,ncol=n)
indic3<-(pred.lin3>0.5)
contour(x,y,indic3,add=T,lwd=2,nlevels=2)

#kNN avec k=30
class.knn<-knn(train=data.class$ech[,1:2],test=test[,1:2],cl=data.class$ech[,3],k=30)
error.knn<-(class.knn==test[,3])
prop.ext.knn30<-length(error.knn[error.knn==T])/length(error.knn)

class.knn<-knn(train=data.class$ech[,1:2],test=plan,cl=data.class$ech[,3],k=30)
plot(test[test[,3]==1,1:2],col="red",xlim=c(-3,4),ylim=c(-2,3),xlab="x",ylab="y")
title(main=c("validation sur un échantillon test","kNN avec k=30"))
points(test[test[,3]==0,1:2],col="blue")
indic.knn<-matrix(class.knn,nrow=n)
contour(x,y,indic.knn,add=T,lwd=2,nlevels=2)


#kNN avec k=10
class.knn<-knn(train=data.class$ech[,1:2],test=test[,1:2],cl=data.class$ech[,3],k=10)
error.knn<-(class.knn==test[,3])
prop.ext.knn10<-length(error.knn[error.knn==T])/length(error.knn)

class.knn<-knn(train=data.class$ech[,1:2],test=plan,cl=data.class$ech[,3],k=10)
plot(test[test[,3]==1,1:2],col="red",xlim=c(-3,4),ylim=c(-2,3),xlab="x",ylab="y")
title(main=c("validation sur un échantillon test","kNN avec k=10"))
points(test[test[,3]==0,1:2],col="blue")
indic.knn<-matrix(class.knn,nrow=n)
contour(x,y,indic.knn,add=T,lwd=2,nlevels=2)

#kNN avec k=1
class.knn<-knn(train=data.class$ech[,1:2],test=test[,1:2],cl=data.class$ech[,3],k=1)
error.knn<-(class.knn==test[,3])
prop.ext.knn1<-length(error.knn[error.knn==T])/length(error.knn)

class.knn<-knn(train=data.class$ech[,1:2],test=plan,cl=data.class$ech[,3],k=1)
plot(test[test[,3]==1,1:2],col="red",xlim=c(-3,4),ylim=c(-2,3),xlab="x",ylab="y")
title(main=c("validation sur un échantillon test","kNN avec k=1"))
points(test[test[,3]==0,1:2],col="blue")
indic.knn<-matrix(class.knn,nrow=n)
contour(x,y,indic.knn,add=T,lwd=2,nlevels=2)

#validation croisée
class.knn.cv<-knn.cv(train=data.class$ech[,1:2],cl=data.class$ech[,3],k=30)
error.knn.cv<-(class.knn.cv==data.class$ech[,3])
prop.cv.knn30<-length(error.knn.cv[error.knn.cv==T])/length(error.knn.cv)

class.knn.cv<-knn.cv(train=data.class$ech[,1:2],cl=data.class$ech[,3],k=10)
error.knn.cv<-(class.knn.cv==data.class$ech[,3])
prop.cv.kn10<-length(error.knn.cv[error.knn.cv==T])/length(error.knn.cv)

class.knn.cv<-knn.cv(train=data.class$ech[,1:2],cl=data.class$ech[,3],k=1)
error.knn.cv<-(class.knn.cv==data.class$ech[,3])
prop.cv.knn1<-length(error.knn.cv[error.knn.cv==T])/length(error.knn.cv)

taux<-matrix(nrow=3,ncol=11)
k<-c(1,2,3,5,10,15,20,25,30,40,50)
for (i in 1:11)
{
	class.knn<-knn(train=data.class$ech[,1:2],test=data.class$ech[,1:2],cl=data.class$ech[,3],k=k[i])
	error.knn<-(class.knn==data.class$ech[,3])
	taux[1,i]<-length(error.knn[error.knn==T])/length(error.knn)
	class.knn<-knn(train=data.class$ech[,1:2],test=test[,1:2],cl=data.class$ech[,3],k=k[i])
	error.knn<-(class.knn==test[,3])
	taux[2,i]<-length(error.knn[error.knn==T])/length(error.knn)
	class.knn<-knn.cv(train=data.class$ech[,1:2],cl=data.class$ech[,3],k=k[i])
	error.knn<-(class.knn==data.class$ech[,3])
	taux[3,i]<-length(error.knn[error.knn==T])/length(error.knn)
}
taux
plot(k,taux[1,],type="b",pch=3,xlab="nombre k",ylab="proportion de bien classés",ylim=c(0.75,1))
title(main="différentes mesures d'erreur")
lines(k,taux[2,],type="b",pch=3,col="green")
lines(k,taux[3,],type="b",pch=3,col="orange")
legend("topright",col=c("black","green","orange"),lty=1,lwd=2,legend=c("apprentissage","externe","validation croisée"))

#############################################################################################
#Septième partie : boostrap et bagging
#############################################################################################

#exemple de la loi normale N(m,1) pour le bootstrap
###################################################
n<-1OO
alpha<-0.05
x<-rnorm(n,mean=5,sd=1)

#IC théorique
a<-mean(x)+(1/sqrt(n))*qnorm(alpha/2)
b<-mean(x)+(1/sqrt(n))*qnorm(1-(alpha/2))

#IC par bootstrap
#NB : le package boot est beaucoup plus sophistiqué

B<-5000
ech.boot<-vector(length=B)
for (i in 1:B)
{
	temp<-sample(x,size=n,replace=T)
	ech.boot[i]<-mean(temp)
}
a.boot<-as.numeric(quantile(ech.boot,probs=alpha/2))
b.boot<-as.numeric(quantile(ech.boot,probs=1-alpha/2))

c(a,a.boot)
c(b,b.boot)


#bagging pour kNN
#################

#on montre que c'est inutile car le modèle est linéaire en les réponses
B<-10


#plan de test pour graphique
n<-200
x<-seq(-3,4,length.out=n)
y<-seq(-2,3,length.out=n)
plan<-expand.grid(x=x,y=y)

class.b1<-matrix(0,nrow=n,ncol=n)
for (i in 1:B)
{
	order<-sample(1:n,size=n,replace=T)
	temp<-data.class$ech[order,]
	class.knn<-as.numeric(knn(train=temp[,1:2],test=plan,cl=temp[,3],k=1))
	class.b1<-matrix(class.knn,nrow=n,ncol=n)+class.b1
}

indic.knn.b1<-(class.b1>3*B/2)
ech<-data.class$ech
plot(ech[ech[,3]==1,],col="red",xlim=c(min(ech[,1]),max(ech[,1])),ylim=c(min(ech[,2]),max(ech[,2])),xlab="x",ylab="y")
title(main=c("séparation des deux classes","kNN avec k=1"))
points(ech[ech[,3]==0,],col="blue")
points(t(data.class$m1),pch=23,col="red",bg="red")
points(t(data.class$m2),pch=23,col="blue",bg="blue")
contour(x,y,indic.knn.b1,add=T,lwd=2,nlevels=2)

class.knn1<-knn(train=data.class$ech[,1:2],test=plan,cl=data.class$ech[,3],k=1)
indic.knn1<-matrix(class.knn1,nrow=n)
contour(x,y,indic.knn,add=T,lty=3,lwd=1,col="orange",nlevels=2)
legend("topright",cex=0.8,legend=c("frontière kNN","frontière après bagging : B=50"),lty=c(1,3),lwd=c(2,1),col=c("black","orange"))