```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, cache = TRUE)
```

## Exercice 4.1 : Factorielle

1. Construction de la fonction factorielle avec la fonction *prod*:

```{r,message=FALSE,warning=FALSE}
factorielle <- function(n) {
  if (n<0) stop("l'entier doit être strictement positif")
  if (n==0) return(1)
  if (floor(n)!=n) warning(paste("arrondi de",n,"en",floor(n)))
  resultat <- prod(1:n)
  return(resultat)
}
```

2. Construction de la fonction factorielle avec une boucle *for*:

```{r,message=FALSE,warning=FALSE}
factorielle <- function(n) {
  if (n<0) stop("l'entier doit être strictement positif")
  if (n==0) return(1)
  if (floor(n)!=n){
    warning(paste("arrondi de",n,"en",floor(n)))
    n <- floor(n)
  }
  resultat <- 1
  for (i in 1:n) resultat <- resultat*i
  return(resultat)
}
factorielle(5)
factorielle(4.23)
```
## Exercice 4.2 : Ventilation

1. Construisons la fonction de ventilation:
```{r,message=FALSE,warning=FALSE}
ventilation <- function(Xqual,p=0.05) {
  if (!is.factor(Xqual)) stop("Xqual doit etre un facteur \n")
  modalites <- levels(Xqual)
  if (length(modalites)<=1) stop("pas assez de niveaux \n")
  tabl <- table(Xqual)
  selecti <- (tabl/sum(tabl))<p
  if (!any(selecti)) return(Xqual) else {
   lesquels <- modalites[!selecti]
   prov <- factor(Xqual[(Xqual%in%lesquels)],levels=lesquels)
   prov <- table(prov)
   proba <- prov/sum(prov)
   for (j in modalites[selecti]) {
    ## tirages dans les modalites au hasard et remplacement
    if (length(lesquels)==1){
     warning("1 seule modalite\n")
     Xqual[Xqual==j]<-lesquels
    } else Xqual[Xqual==j]<-sample(lesquels,sum(Xqual==j),
             replace=TRUE,prob=proba)
   }
   Xqualvent <- factor(as.character(Xqual))
  }
  return(Xqualvent)
 }
```
2. Nous appliquons la fonction de la question précédente à chaque colonne du tableau qui est un facteur:
```{r,message=FALSE,warning=FALSE}
ventil.tab <- function (tab, seuil=0.05) {
 for (i in 1:ncol(tab)) {
  if (is.factor(tab[,i])) tab[,i]<-ventilation(tab[,i],p=seuil)
 }
 return(tab)
}
```

## Exercice 4.3 : Ventilation sur facteur ordonné

1. La fonction est identique ou presque à ce qui a été vu dans la correction de l'exercice 2.
Nous ajoutons la sortie (*return*) et quelques contrôles :
```{r,message=FALSE,warning=FALSE}
ventilation.ordonnee <- function(Xqual,p=0.05) {
  if (!is.ordered(Xqual)) stop("Xqual doit etre un ordonne \n")
  modalites <- levels(Xqual)
  if (length(modalites)<=1) stop("pas assez de niveaux \n")
  tabl <- table(Xqual)
  selecti <- (tabl/sum(tabl))<p
  if (!any(selecti)) return(Xqual) else {
    numero <- which(selecti)
    while(any((tabl/sum(tabl))<p)) {
      Xqual[,i] <- ventilation(Xqual[,i])
   }
 }
 return(Xqual)
}
```
2. Nous appliquons la fonction de la question précédente à chaque colonne du tableau qui est un facteur ordonné:
```{r,message=FALSE,warning=FALSE}
ventil.ordonnee.tab <- function (tab, seuil=0.05) {
 for (i in 1:ncol(tab)) {
  if (is.ordered(tab[,i])) tab[,i]<-ventilation.ordonnee(tab[,i],seuil)
 }
 return(tab)
}
```

## Exercice 4.4 : Parallélisation

```{r}
max_inc <- function(seq) {
  dseq <- seq[-1] > seq[-length(seq)]; cur_l <- 1+dseq[1]; max_l <- cur_l
  for (i in 1:(length(dseq)-1)) {
    if (dseq[i] && (dseq[i] == dseq[i+1])) {
    cur_l <- cur_l + 1;
    } else { cur_l <- 1+dseq[i+1] }
  max_l <- max(cur_l, max_l)
  }
  return(max_l)
}
```

1. Utilisation de sapply et table

```{r}
n <- 1000
k <- 100000
```

```{r}
system.time({
  set.seed(234)
  res <- sapply(1:k, function(aux){max_inc(runif(n))})
})
table(res)
```
2. Utilisation du package parallel

```{r,message=FALSE,warning=FALSE}
require(parallel)
system.time({
  cl <- makeCluster(detectCores()-1)    
  clusterExport(cl, varlist=c("n","k","max_inc"))
  res <- parSapply(cl, 1:k, function(aux){max_inc(runif(n))})
  stopCluster(cl)
})
table(res)
```

3. Pour accelérer encore les temps de calcul, il faut utiliser le package Rcpp.

```{Rcpp}
#include <Rcpp.h>
using namespace Rcpp;


// [[Rcpp::export]]
int max_inc_rcpp(NumericVector x) {
  int max_l = 0;
  int cur_l = 0;
  bool cur_g, prev_g;
  cur_g = x[1] > x[0];
  cur_l = 1 + cur_g;
  max_l = cur_l;
  for (int i = 1; i < (x.length()-1); i++) {
    prev_g = cur_g;
    cur_g = x[i+1] > x[i];
    if (cur_g && (prev_g == cur_g)) {
      cur_l = cur_l + 1;
    } else {
      cur_l = 1 + cur_g;
    }
    if (cur_l > max_l) { max_l = cur_l;};
  }
  return max_l;
}
```

```{r}
system.time(res <- sapply(1:k, function(i) { max_inc_rcpp(runif(n)) }))
```

## Exercice 4.5 : 

1. Simulons deux variables x et y:

```{r}
set.seed(1234)
x <- rnorm(50)
epsilon <- rnorm(50)
y <- 2+ 0.5*x + epsilon
```

2. On calule le coefficient entre x et y et le test associé, puis on compte le nombre de fois où le coefficient de  corrélation entre x et y permuté est plus grand que le coefficient de corrélation observé:

```{r}
cor.xy=cor(x,y)
cor.test(x,y)
nbsim=100000

system.time({
set.seed(234)
som=0
for (i in 1:nbsim) som=som+(cor(x,y[sample(length(y))])>cor.xy)
})
print(som/nbsim*2)
```
La probabilité critique obtenue par le test classique, 0.01484, est proche de celle calculée par permutation 0.01476.

3. On utilise le package parallel et la fonction parLapply

```{r,message=FALSE,warning=FALSE}
require(parallel)
cl <- makeCluster(3)
clusterExport(cl, varlist=c("x","y","cor.xy"))
system.time({
  res <- parSapply(cl, 1:nbsim, function(aux){1*(cor(x,y[sample(length(y))])>cor.xy)})
})
stopCluster(cl)
mean(res)*2
```
On multiplie par 2 le pourcentage obtenu car le test est bilatéral.
On note que le temps de calcul est en gros divisé par 2.6.

4. On relance les lignes de code avec $rchisq$ pour générer les aléas:

```{r}
set.seed(1234)
x <- rchisq(50, df=1)
epsilon <- rchisq(50, df=1)
y <- 2+ 0.5*x + epsilon

cor.xy=cor(x,y)
cor.test(x,y)
nbsim=100000

system.time({
set.seed(234)
som=0
for (i in 1:nbsim) som=som+(cor(x,y[sample(length(y))])>cor.xy)
})
print(som/nbsim*2)

require(parallel)
cl <- makeCluster(3)
clusterExport(cl, varlist=c("x","y","cor.xy"))
system.time({
  set.seed(234)
  res <- parSapply(cl, 1:nbsim, function(aux){1*(cor(x,y[sample(length(y))])>cor.xy)})
})
stopCluster(cl)
2*mean(res)
```
La probabilité critique obtenue avec le test classique, 0.007628, est cette fois assez différente de celle obtenue par le test de permutation. En effet, l'hypothèse de normalité n'est pas vérifiée ici et le test classique n'est pas adaptée. On aura donc plus confiance dans l'estimation de la probabilité critique obtenue par le test de permutation : 0.025