UP | HOME

Séance 7 : Algo & Prog avec R

Table des matières

1 Recherche d’éléments

  1. Déterminez le nombre de 0 présent dans un vecteur.
  2. Déterminer l’indice du premier 0 d’un vecteur.
  3. Déterminez le plus petit élément non nul d’un vecteur.
  4. Déterminer la valeur médiane d’un vecteur.
vec <- sample(0:10,12, replace=TRUE)
print(vec)
sum(vec == 0)
head(which(vec == 0),1)
min(vec[ vec != 0])
median(vec)

2 Calcul de l’hypoténuse

  1. Calculez l’hypoténuse d’un triangle connaissant les deux côtés de l’angle droit (a et b).
  2. Demander à l’utilisateur de saisir les valeurs a et b à l’aide de la fonction scan()
  3. Afficher et formatter le résultat avec les fonctions cat et paste.

    cat("Veuillez entrer les longueurs des deux côtés de l’angle droit :\n")
    ab <- scan( nmax=2)
    if(length(ab) == 2) {
      hypo <- sqrt(ab[1]**2 + ab[2]**2)
      print(paste("Longueur de l'hypoténuse :", hypo))
    } else {
      print("Erreur de saisie")
    }
    

3 Jeu de hasard

Virginie lance trois dés numérotés de 1 à 6.

  • Si elle obtient une somme de 18, elle gagne 50 euros,
  • entre 10 et 17, elle gagne 5 euros,
  • sinon elle ne gagne rien.
  1. Écrivez une fonction JeuHasard utilisant la fonction sample pour simuler un lancer de dés, puis renvoyant le gain.
    JeuHasard <- function() {
      somme <- sum(sample(1:6, 3, replace=TRUE))
      if(somme <10) {return(0)}
      else if(somme <18) {return(5)}
      else {return(50)}
    }
    
  2. Écrire une simulation où Virginie joue jusqu’à ce que son gain dépasse 50.
    gain <- 0
    partie <- 0
    while(gain < 50) {
      gain <- gain + JeuHasard()
      partie <- partie + 1
      cat(paste("Partie", partie, ":", gain,"\n"))
    }
    
    Partie 1 : 0 
    Partie 2 : 0 
    Partie 3 : 5 
    Partie 4 : 5 
    Partie 5 : 10 
    Partie 6 : 10 
    Partie 7 : 15 
    Partie 8 : 20 
    Partie 9 : 25 
    Partie 10 : 30 
    Partie 11 : 30 
    Partie 12 : 35 
    Partie 13 : 35 
    Partie 14 : 40 
    Partie 15 : 40 
    Partie 16 : 40 
    Partie 17 : 45 
    Partie 18 : 50
    
    tirages <- expand.grid(1:6, 1:6, 1:6)
    sommes <- rowSums(tirages)
    n5 <- sum( sommes >= 10) - 1
    cat(paste("Esperance :", 50/6**3 + 5*n5/nrow(tirages)))
    
    Esperance : 3.33333333333333
    
  3. Quelle est la probabilité de gagner 50 euros ? Quelle est l’espérance de gain ? Proposer un tarif pour jouer à ce jeu ? Justifier.
    prob50 <- 1 / 6**3
    ## Estimation de l'espérance par simulation
    n <- 10000
    gains <- replicate(n, JeuHasard())
    cat(paste("Esperance simulée :", sum(gains)/n, "\n")) 
    ## Calcul théorique de l'espérance
    tirages <- expand.grid(1:6, 1:6, 1:6)
    sommes <- rowSums(tirages)
    prob5 <- (sum(sommes >= 10) - 1)/ nrow(tirages)
    cat(paste("Esperance théorique", 50*prob50 + 5*prob5), '\n')
    
    Esperance simulée : 3.3625
    Esperance théorique 3.33333333333333
    

4 La baguette chanceuse

J’ai dans ma poche 2 pièces de 10 centimes, 1 pièce de 20 centimes, 3 pièces de 50 centimes, 1 pièce de 1 euro et 2 pièces de 2 euro. Quelle est la probabilité pour qu’en sortant deux pièces au hasard de ma poche, je puisse payer une baguette à 1 euro. ?

  1. Estimez la probabilité par simulation. Indice: replicate
  2. Calculez la probabilité théorique.
money <- c(10, 10, 20, 50, 50, 50, 100, 200, 200)
can_pay <- function(take, money, amount) {
  sum(sample(money, take, replace=FALSE)) >= amount
} 
n <- 10**5
proba <- mean(replicate(n, can_pay(2, money, 100)))
cat(paste("Simulation probability :", round(proba,4)),"\n")

## Theory : 
## Winning combinations : 
## 100 or 200 x 1 of the 8 remaining coins
wins <- choose(3,1)*choose(8,1) 
## Combinations with only 100 or 200 are counted twice
wins <- wins - choose(3,2) 
## 2 x 50.
wins <- wins + choose(3,1) 
total <- choose(length(money), 2)
cat(paste("Theoritical probability:", round(wins/total,4)),"\n")
Simulation probability : 0.6673 
Theoritical probability: 0.6667 

5 Nombres aléatoires absents

  1. Utilisez la fonction sample pour créer une liste de 100 entiers tirés au hasard entre 0 et 99.
  2. Calculer le nombre d’éléments compris entre 0 et 99 qui n’appartiennent pas à cette liste.
  3. Recommencer cette expérience un grand nombre de fois et calculer la moyenne du nombre d’absents.
  4. Comparer à la valeur théorique.
n <- 100 
elems <- 0:99

xp <- 1000
absents <- 0
for(i in 1:xp) {
  li <- sample(elems, n, replace=TRUE)
  absents <- absents + sum( !(elems %in% li) )
}
print(absents/xp)
[1] 36.697

6 Coefficient de Gini

  1. Entrer les données ci-dessous dans R grâce à la fonction scan.

    Nous considérons trois échantillons aléatoires d’une population :

    1. le nombre d’heures de charge pour différents modèles de téléphone mobile ;
      45.8   41.1   55.9   46.6   57.0   45.0   58.5   46.7   49.3   52.7   
      54.9   48.5   40.4   44.4   51.0   44.2   59.1   46.9   50.7   43.7   
      41.7   52.8   60.5   38.5   60.4   53.8   47.3   50.2   58.8   50.7
      
    2. le nombre d’heures passées devant la télévision pour 34 foyers ;
      23.1   15.9   21.0   26.0   25.1   14.7   24.2   16.6   18.2   16.5   
      20.7   15.3   17.7   19.1   22.7   21.9   14.6   26.3   25.8    9.4   
      17.0   21.2   17.9   24.7   21.1   17.2   19.1   22.7   24.0   24.7   
      22.5    8.3   2.5    30.4
      
    3. le nombre de km pour aller et revenir du travail de 12 employés.
      3.7    14.3   11.0   26.5    5.2    4.8   24.2   16.9   8.2    26.5   
      40.7   5.3
      
  2. Le coefficient de Gini est donné par la formule (\(x_i \leq x_{i+1}\)):

    \[ \frac{2}{n} \frac{ \sum_1^n i \times x_i}{\sum_1^n x_i} - \frac{n+1}{n}. \] Écrire une fonction GiniIndex qui calcule le coefficient de Gini. Calculer les coefficients de Gini pour les trois échantillons ci-dessus. Qu’en concluez-vous ?

    GiniIndex(phones)
    GiniIndex(tv)
    GiniIndex(dist)
    
    [1] 0.07121991
    [1] 0.1539792
    [1] 0.3894376
    
  3. Implémenter les fonctions LorenzCurve(data) renvoyant une matrice (ou dataframe) contenant les coordonnées des points de la courbe de Lorenz pour l’échantillon data.
    LorenzCurve(dist)
    
                x          y
    1  0.00000000 0.00000000
    2  0.08333333 0.01975440
    3  0.16666667 0.04538174
    4  0.25000000 0.07314469
    5  0.33333333 0.10144154
    6  0.41666667 0.14522157
    7  0.50000000 0.20395088
    8  0.58333333 0.28029899
    9  0.66666667 0.37052856
    10 0.75000000 0.49973305
    11 0.83333333 0.64121730
    12 0.91666667 0.78270155
    13 1.00000000 1.00000000
    
  4. Utiliser l’extension ggplot2 pour produire des graphiques de la courbe de Lorenz.

    Tout d’abord, il faut installer quelques extensions.

    install.packages(c("ggplot2", "ggthemes", "Rmisc"))
    

    Ensuite, définir les fonctions construisant les graphiques à partir d’un jeu de données.

    library(ggplot2)
    library(ggthemes)
    PlotLorenz <- function(lorenzCurve, giniIndex, dataName) {
      ggplot() + 
        geom_line(aes(y = y, x = x),
                  data = as.data.frame(lorenzCurve), stat="identity", 
                  color =  ptol_pal()(1), size = 2) +
        geom_segment(aes(x = 0, y = 0, xend = 1, yend = 1),linetype="dotted", size = 1.5) +
        theme_gdocs() + scale_colour_ptol() +
        ggtitle(paste("Courbe de Lorenz - Gini Index =", format(giniIndex, digits=2)), 
                sub = dataName) +
        labs(x="Part cumulée de la population", y=paste("Part cumulée :", dataName)) +  
        theme(text = element_text(size=14), plot.title = element_text(size=14))
    }
    
    PlotLorenzFromSample <- function(data, dataName) {
      PlotLorenz(LorenzCurve(data), GiniIndex(data), dataName)
    }
    
    phonesName <- "temps de charge"
    tvName <- "temps passé devant la TV"
    distName <- "distance du lieu de travail"
    

    Ensuite, exécuter le code ci-dessous pour afficher les courbes de Lorenz des trois échantillons.

    library(Rmisc)
    multiplot(
      PlotLorenzFromSample(phones, phonesName),
      PlotLorenzFromSample(tv, tvName),
      PlotLorenzFromSample(dist, distName), cols =3
    )
    

    lorenz.jpg

  5. Utiliser l’extension shiny pour produire une application web interactive.

    Installer d’abord l’extension.

    install.packages("shiny")
    

    Ensuite, rassembler le code définissant les nécessaire dans un fichier et exécuter le. N’hésitez pas à lire l’introduction à shiny.

    library(shiny)
    if (interactive()) {
      ## Only run examples in interactive R sessions
      choices <- c("phones", "tv", "dist")
      names(choices) <-c(phonesName, tvName, distName) 
    
      ## user interface object
      ui <- fluidPage(
        titlePanel("Ma première application web avec shiny!"),
        sidebarLayout(
          ## Sidebar with radio buttons
          sidebarPanel(
            radioButtons("data", "Données : ", choices)
          ),
          ## Main panel with graphic
          mainPanel(
            plotOutput("dataPlot")
          )
        )
      )
    
      ## a server function
      server <- function(input, output) {
        ## Build plot objects in advance
        phonesPlot <- PlotLorenzFromSample(phones, phonesName)
        tvPlot <- PlotLorenzFromSample(tv, tvName)
        distPlot <- PlotLorenzFromSample(dist, distName)
    
        ## Set the right plot object according to the user choice
        output$dataPlot <- renderPlot({
          switch(input$data,
                 phones = phonesPlot,
                 tv = tvPlot,
                 dist = distPlot,
                 phonesPlot)
        })
      }
      ## a call to the shinyApp function
      shinyApp(ui, server)
    }
    

7 Racines d’un trinôme   HARD

Il est fortement conseillé de lire la page wikipedia sur les équations du second degré ou encore mieux Scilab is not naive.

7.1 Méthode naïve

Écrire une fonction triroot(a,b,c) prenant en paramètre les coefficients a, b et c d’un trinôme et renvoyant les racines réelles de l’équation \(ax^2 + bx + c =0\). Plus précisément, la fonction renvoie :

  • NULL si l’équation n’admet pas de racines réelles ;
  • un scalaire si l’équation admet une racine double ;
  • un vecteur à deux éléments si l’équation admet deux racines distinctes.

    triroot <- function(a, b, c) {
      delta <- b**2 - 4*a*c
      if(delta < 0) {
        return(NULL)
      } else if(delta > 0) {
        s <- sqrt(delta)
        return(c( -b + s, -b - s)/(2*a))
      } else {
        return -b/(2*a)
      }
    }
    

Vérifiez votre programme en utilisant la fonction prédéfinie polyroot.

triroot(1,-3,2)
polyroot(c(2,-3,1))
[1] 2 1
[1] 1+0i 2-0i

7.2 Autour de la validité de la comparaison à 0

Tester votre programme avec le code ci-dessous. Quelle conclusion en tirez-vous?

triroot(0.01,0.2,1)
polyroot(c(1,0.2,0.01))
triroot(0.011025,0.21,1)
polyroot(c(1,0.21,0.01025))
[1] -10 -10
[1] -10-0i -10+0i
NULL
[1]  -7.52792+0i -12.95989-0i

Proposer une amélioration permettant d’éviter le problème ci-dessus grâce à la fonction all.equal.

7.3 Autour de l’annulation massive

On va analyser l’erreur d’arrondi pendant le calcul du discriminant quand \(b^2 >> 4ac\) en étudiant le trinôme \(\epsilon x + \frac{x}{\epsilon} - \epsilon\).

li <- head(seq(0.0001,0,length.out=4),-1)
for (epsilon in li) {
  e1 <- epsilon**2
  #print(paste("Epsilon:", epsilon))
  print(paste("Expected root:",e1))
  r1 <- triroot(epsilon, 1/epsilon, -epsilon)[1]
  print(paste("Naive method:",r1, "error=", format(abs(1-r1/e1), digits=3)))
  r2 <-  Re(polyroot(c(-epsilon,1/epsilon,epsilon))[1])
  print(paste("R method:",r2, "error=", format(abs(1-r2/e1),digits=3)))
}
[1] "Expected root: 1e-08"
[1] "Naive method: 9.09494701772928e-09 error= 0.0905"
[1] "R method: 1e-08 error= 0"
[1] "Expected root: 4.44444444444445e-09"
[1] "Naive method: 1.36424205265939e-08 error= 2.07"
[1] "R method: 4.44444444444445e-09 error= 0"
[1] "Expected root: 1.11111111111111e-09"
[1] "Naive method: 0 error= 1"
[1] "R method: 1.11111111111111e-09 error= 0"
triroot <- function(a, b, c) {
  b <- b/2
  d1 <- b**2 
  d2 <- a*c
  if(isTRUE(all.equal(d1,d2))) {
    return(-b/a)
  } else if(d1 < d2) {
    return(NULL)
  } else {
    h <- -(b + sign(b)*sqrt(d1-d2))
    return(c(c/h, h/a))
  } 
}

7.4 TODO Autour du dépassement de capacité

7.5 STARTED Extension aux racines complexes

Vous avez de la chance, il existe une classe complex d’emblée intégrée à R, donc sans import.

  1. Cherchez dans la documentation ou sur le web comment définir les nombres complexes \(z_1=3-2i\), \(z_2=5+i\) et \(z_3=i\).
  2. Calculez l’addition \(z_1 + z_2\), le produit \(z_1 \times z_2\) et l’inverse $\frac{1}{z_3}. Vérifiez que \(z_3\) est bien la racine carrée de -1 …
  3. Reprogrammez la méthode triroot retournera aussi les racines complexes.
  4. Calculez les racines des polynômes \(2x^2 - 3x -2\), \(x^2+x+1\), et \(x^2+1\).
z1 <- 3-2i
z2 <- 5+1i
z3 <- 1i
z1 + z2
z1 * z2
1/z3
[1] 8-1i
[1] 17-7i
[1] 0-1i

8 Jeu de Nim (Fibonacci)   HOME HARD

La page de Jean-Paul Dalavan peut vous aider.

8.1 Décomposition de Zeckendorf   KEY

Étant donné un entier naturel \(n\), déterminer la liste des facteurs de la décomposition de Zeckendorf de \(n\), c’est-à-dire les uniques nombres de Fibonacci deux à deux distincts et non consécutifs de somme égale à \(n\)

8.2 Jeu de Nim

Deux joueurs tirent à tour de rôle des allumettes d’une boîte, avec les règles suivantes :

  • Chaque joueur tire à chaque fois au moins une allumette.
  • Le premier joueur ne retire pas la totalité des allumettes au premier tour.
  • Un joueur tire au plus deux fois le nombre d’allumettes tirées par le joueur précédent.
  • Le joueur qui retire la dernière allumette a gagné.

On peut montrer que si le nombre initial d’allumettes n’est pas un nombre de Fibonacci, une stratégie gagnante pour le joueur 1 consiste à tirer autant d’allumettes que le plus petit terme de la décomposition de Zeckendorf du nombre d’allumettes.

Écrire une fonction prenant en argument un entier n (nombre d’allumettes), et mettant en place la stratégie gagnante si c’est possible (on alternera le jeu de l’ordinateur, qui commence, et le jeu de l’utilisateur, à qui on demandera une valeur, en précisant dans quel intervalle il a le droit de tirer).

zeckendorf <- function() {
  fibo <- c(1,1)
  z <- function(x){
    ## Compute mandatory fibonnacci terms
    while(x > tail(fibo,1) ) {
      ## Here, <<- will assign to the parent environment
      fibo <<- append(fibo, sum(tail(fibo,2)))
    }
    ## Greedy computation of zeckendorf decomposition
    i <- length(fibo)
    zeck <- c()
    while( x > 0) {
      ## find largest Fibonacci number smaller than x
      while(fibo[i] > x) {
        i <- i - 1
      }
      zeck <- append(zeck, fibo[i])
      x <- x - fibo[i]
      i <- i - 2
    }
    return(zeck)
  }
  return(z)
}

z <- zeckendorf()
z(20)
z(30)
z(8)
z(100)
[1] 13  5  2
[1] 21  8  1
[1] 8
[1] 89  8  3
play_game <- function(n, player1, player2, verbose=TRUE) {
  players <- c(player1, player2)
  pidx = 1
  l <- 0; ##last stroke
  repeat {
    ## Player pidx plays 
    lim <- ifelse(l == 0, n-1,min(2*l, n))
    take <- players[[pidx]](n,lim)
    if(1 <= take && take <= lim) {
      ## valid stroke
      n <- n - take
      l <- take
      if(verbose) {
        cat(paste("Player", pidx, "takes", take,"matche(s) :", n, "left.\n"))
      }
      ## Termination condition 
      if(n <= 0) {
        break
      }
    } else {
      cat(paste("Player", pidx, "takes", take,"matche(s) : ERROR!\n"))
      ## Swtich players
      pidx <- ifelse(pidx == 1, 2, 1) 
      break
    }

    ## Swtich players
    pidx <- ifelse(pidx == 1, 2, 1) 
  }
  if(verbose) {
    cat(paste("Player", pidx, "wins!\n"))
  }
  return(pidx)
}

player.one <- function(n, max) 1

player.random <- function(n, max) {
  ifelse(max > n, n, sample(1:max,1))
}

player.human <- function(n, max) {
  take <- ""
  while( ! grepl("^[0-9]+$",take)) {
    take <- readline(prompt="Enter an integer: ")
  }
  return(as.integer(take))
}

z <- zeckendorf()
player.zeckendorf <- function(n, max) {
  if(max >= n) {return(n)}
  else {
    take <- tail(z(n),1)
    if(take <= max) {return(take)}
    else {return(1)}
  }
}

Created: 2017-12-08 ven. 11:01