---
title: "Régle et classifieur de Bayes"
author: "JMM"
date: "16 septembre 2020"
output: pdf_document
---

# On fixe les paramètres associés au modèle de mélange de 2 lois gausiennes que nous allons étudier

On s'intéresse à un problème de classification binéaire
$Y\in\{0,1\}$. On fixe $\pi_0=\mathbb{P}(Y=0)=3/10$ et 
$\pi_1=\mathbb{P}(Y=1)=7/10$. Par ailleurs, on suppose que la variable
$X$ conditionnellement à la classe $Y=k$ est distribuée suivant une loi
gausienne en dimension 2
$$
X|Y=k\sim\mathcal{N}_2\left(\mu_k,\Sigma_k\right)
$$
avec

$\mu_0=(3,3)$, $\mu_1=(5,5)$,
$\Sigma_0=\left[\begin{array}{cc}
 5   & 45/10 \\
 45/10 & 5
\end{array}\right]$ et 
$\Sigma_1=\left[\begin{array}{cc}
 5      & -45/10 \\
 -45/10 & 5
\end{array}\right]$.

```{r}
pi0 <- 0.3 ; pi1 <- 0.7
mu0 <- c(3,3) ; mu1 <- c(5,5)
Sigma0 <- matrix(c(5,45/10,45/10,5),2,2)
Sigma1 <- matrix(c(5,-45/10,-45/10,5),2,2)
```

# On trace la densité de la loi marginale de $X$

Pour un tel modèle la loi marginale de $X$ est un mélange des 2 lois
gausiennes décrites ci-dessus.

```{r}
library(mvtnorm)

len <- 20
x1 <- seq(-5,10,length=len) 
x2 <- seq(-5,10,length=len) 
z <- matrix(0,len,len)
for (i in 1:len) for (j in 1:len)
z[i,j] <- pi0*dmvnorm(c(x1[i],x2[j]),mu0,Sigma0)+
	pi1*dmvnorm(c(x1[i],x2[j]),mu1,Sigma1)
contour(x1,x2,z,nlevels=10)
```

# On trace la probabilité de Y=0 sachant X

```{r}
pi0x <- function(x)
{
c0 <- pi0*dmvnorm(x,mean=mu0,sigma=Sigma0)	
c1 <- pi1*dmvnorm(x,mean=mu1,sigma=Sigma1)	
c0/(c0+c1)
}

aff0 <- matrix(0,len,len)
for (i in 1:len) for (j in 1:len)
aff0[i,j] <- pi0x(c(x1[i],x2[j]))

print(pi0x(c(5,5)))

image(x1,x2,aff0,col=gray((200:0)/200))
# library(rgl)
# persp3d(x1,x2,aff0)
```

# On trace les frontières du classifieur de Bayes

```{r}
cbayes <- matrix(1,len,len)
cbayes[aff0>=1/2] <- 0
image(x1,x2,cbayes,col=gray((0:200)/200))
contour(x1,x2,aff0,levels=0.5,col="red",add=TRUE)
```

# Simulation de données

Un échantillon d'apprentissage 

```{r}
n <- 1000
y <- sample(c(0,1),n,prob=c(pi0,pi1),replace=TRUE)
x <- matrix(0,n,2)
for (i in 1:n)
{
if (y[i]==0) x[i,] <- rmvnorm(1,mu0,Sigma0)
else x[i,] <- rmvnorm(1,mu1,Sigma1)
}
traindat <- data.frame(y=y,x1=x[,1],x2=x[,2])
plot(traindat$x1,traindat$x2,col=traindat$y+1,pch="*",xlab=expression(x[1]),
     ylab=expression(x[2]))
```

Un échantillion de test

```{r}
ntest <- 500
y <- sample(c(0,1),ntest,prob=c(pi0,pi1),replace=TRUE)
x <- matrix(0,ntest,2)
for (i in 1:ntest)
{
if (y[i]==0) x[i,] <- rmvnorm(1,mu0,Sigma0)
else x[i,] <- rmvnorm(1,mu1,Sigma1)
}
testdat <- data.frame(y=y,x1=x[,1],x2=x[,2])
```


# Classification par regression

```{r}
model.reg <- lm(y~x1+x2,data=traindat)
a <- (0.5-model.reg$coefficients[1])/model.reg$coefficients[3]
b <- -model.reg$coefficients[2]/model.reg$coefficients[3]
contour(x1,x2,aff0,levels=0.5,col="red")
abline(a=a,b=b,col="blue",lwd=4)

res.reg <- predict(model.reg,testdat)
classif.reg <- rep(0,ntest)
classif.reg[res.reg > 0.5] <- 1
mean(classif.reg!=testdat$y)
```

# Méthodes des $k$-plus proches-voisins

```{r}
library(caret)
model.knn <- knn3(as.factor(y)~x1+x2,data=traindat,k=7)
predict(model.knn,data.frame(x1=-2,x2=4))

grille <- expand.grid(x1=x1,x2=x2)

res.knn <- predict(model.knn,grille,type="class")
contour(x1,x2,aff0,levels=0.5,col="red",lwd=4)
contour(x1,x2,matrix(as.numeric(res.knn)-1,len,len),add=TRUE,levels=0.5,col="pink",lwd=4)

res.knn <- predict(model.knn,testdat,type="class")
mean((as.numeric(res.knn)-1)!=testdat$y)

```

# Classifieur oracle

```{r}
res.oracle <- rep(0,ntest)
for (i in 1:ntest)
  res.oracle[i] <- pi0x(c(testdat$x1[i],testdat$x2[i]))
res.oracle <- (res.oracle <= 0.5)
mean(res.oracle!=testdat$y)
```

# knn choix de $k$ par validation croisée

```{r, cache=TRUE}
model.knn <- train(as.factor(y)~x1+x2,data=traindat,method="knn",metric="Accuracy",preProcess=c("center", "scale"),trControl=trainControl(method="repeatedcv",number=2,repeats=100),tuneGrid=data.frame(k=1:25))
```

```{r}
plot(model.knn)
print(model.knn)
res.knn <- predict(model.knn,testdat)
mean((as.numeric(res.knn)-1)==testdat$y)
```