#############"Isomorphismes" D'EXPÉRIENCES ALÉATOIRES SIMULÉES ##########
#********LIGNES DE COMMANDE pile ou face pas forcément équiprobable **********
piece1 <- sample(c("Pile", "Face"), size = 1000, prob = c(.4, .6), replace = TRUE)
(distpiece1 <- table(piece1))
sum(piece1 == "Pile") / 1000
barplot(distpiece1 / 1000)

#--------------------------------------
#******LIGNES DE COMMANDE pile ou face avec 2 pièces différentes *****
piece1 <- sample(c("Pile", "Face"), size = 1000, prob = c(.4, .6), replace = T)
piece2 <- sample(c("Pile", "Face"), size = 1000, prob = c(.5, .5), replace = T)
deuxpieces <- paste(piece1,piece2, sep = "")
table(deuxpieces)
barplot(table(deuxpieces) / 1000)
#-----------------------------------------------
#*****FONCTION jet simultané de 2 pièces identiques équilibrées****
pileface <- function(nbsim = 2000){
 resultats <- vector(length = 3)
 names(resultats) <- c("deuxpils", "deuxfaces", "autre")
 for(i in 1:nbsim){
 pieceA <- sample(c("Pile", "Face"), 1)
 pieceB <- sample(c("Pile", "Face"), 1)
 if(pieceA == "Pile" & pieceB == "Pile") {resultats[1] <- resultats[1] + 1} else {
   if(pieceA == "Face" & pieceB == "Face") {resultats[2] <- resultats[2] + 1} else {
    resultats[3] <- resultats[3] + 1                                               }
                                                                                 }
                  }
 print(resultats)
 print(resultats / nbsim)
 barplot(resultats / nbsim)
}

#-----------------------------------------------
#***Le problème historique du grand duc de Toscane ("Somme" de 3 dés)**
#*****LIGNES DE COMMANDE Simulation GrandDuc*****
de1 <- sample(c(1:6), 1000, replace = TRUE)
  (distde1 <- table(de1))
  barplot(distde1 / 1000)
de2 <- sample(c(1:6), 1000,replace = T)
  (distde2 <- table(de2))
  dev.new()
  barplot(distde2 / 1000)
de3 <- sample(c(1:6), 1000,replace = T)
  (distde3 <- table(de3))
  dev.new()
  barplot(distde3 / 1000)
de <- de1 + de2 + de3
(distde <- table(de))
  dev.new()
  barplot(distde / 1000)
  nbneuf <- sum(de == 9)
  nbdix <- sum(de == 10)
  cat("Fréquence des neuf =", nbneuf / 1000, "\n")
  cat("Fréquence des dix =", nbdix / 1000, "\n")
barplot(distde, xlab = "Somme des numéros des 3 faces",
 ylab = "Effectifs simulés",
 main = paste("Diagramme en barre de 1000 simulations\n du jet de 3 dés équilibrés"))
#--------------------------------------------------
#*****FONCTION simulation Grand Duc*************
simgrandduc <- function(nbsim=1000){
 de1 <- sample(c(1:6), nbsim, replace = TRUE)
 de2 <- sample(c(1:6), nbsim, replace = TRUE)
 de3 <- sample(c(1:6), nbsim, replace = TRUE)
 de <- de1 + de2 + de3
 distde <- table(de)
 print(distde)
 barplot(distde / nbsim)
}

#---------------------------------------------------
#**LIGNES DE COMMANDE probabilité Grand Duc********
#**Somme des valeurs des faces obtenues en jetant 3 dés**
#***Calculer avec le modèle mathématique "exact"****
#****Construire l'univers correspondant à cette expérience***
serieS3de <- array(data = NA, dim = c(6, 6, 6))
 for(i in 1:6){
   for(j in 1:6){
     for(k in 1:6){
        serieS3de[i, j, k] <- i + j + k
                  }
                }
            }
 serieS3de
 (distS3de <- table(serieS3de))
 nbneuf <- sum(serieS3de == 9)
 nbdix <- sum(serieS3de == 10)
 cat("Probabilité de neuf =",nbneuf / 216,"\n")
 cat("Probabilité de dix =",nbdix / 216,"\n")
 dev.new()
 barplot(distS3de)
 graphics.off()
#-------------------------------------------------
#**Fonction probabilité Grand Duc********
probgranduc <- function(){
serieS3de <- array(data = NA, dim = c(6, 6, 6))
 for(i in 1:6){
   for(j in 1:6){
     for(k in 1:6){
        serieS3de[i, j, k] <- i + j + k
                  }
                }
            }
 serieS3de
 distS3de <- table(serieS3de)
 nbneuf <- sum(serieS3de == 9)
 nbdix <- sum(serieS3de == 10)
 cat("Probabilité de neuf =",nbneuf / 216,"\n")
 cat("Probabilité de dix =",nbdix / 216,"\n")
 print(distS3de)
 barplot(distS3de / 216)
}

#---------------------------------------
# **LIGNES DE COMMANDES ****** SIMULATIONS NUMÉRIQUES *****
# Illustration graphique de la loi des grands nombres :
# Lorsque n augmente, on observe les suites de distributions
# Les fréquences tendent vers une valeur limite : la probabilité
# Les écarts à cette valeur limite sont de plus en plus faibles
#
nechant <- rep(c(10, 20, 50, 100, 500, 1000),
               c(20, 20, 20, 20, 20, 20))
simfreqar <- c(rbinom(20, 10, .3)/10,
             rbinom(20, 20, .3)/20, rbinom(20, 50, .3)/50,
             rbinom(20, 100, .3)/100, rbinom(20, 500, .3)/500,
             rbinom(20, 1000, .3)/1000)
dev.off()
plot(nechant, simfreqar, xaxp = c(0, 1000, 10),
  main = "Distributions des fréquences des succès")
dev.new()
plot(as.factor(nechant), simfreqar,
  xlab = "Échantillons par tailles",
  ylab = "Fréquence des succès",
  main = "Résumé des distributions")

################ CALCULS DE PROBABILITÉS ################
#-------------- Loi binomiale ---------------
# Calcul de P(A = X = B), X étant une v.a. de distribution binomiale
# de paramètres n=100 et p=0,52.
# Les exemples choisis peuvent servir de base à une réflexion sur
# les différentes façons de déterminer un intervalle de fluctuation, à partir
# de l'exemple 1 (Monsieur Z du document d'inspection).
# P(42<=X<=62):
sum(dbinom(42:62, 100, .52))
# P(43<=X<=62):
sum(dbinom(43:62, 100, .52))
# P(42<=X<=61):
sum(dbinom(42:61, 100, .52))
# P(X<=41) ; P(X<=42) ; P(X<=43):
pbinom(41:43, 100, .52)
#-------------- Combinaisons et loi hypergéométrique ---------------
# Calcul de P(X=3) ; x=3 ; X de loi hypergéométrique de paramètres m=3, n=5, k=4,
(proba <- choose(3, 3)*choose(5, 4-3)/choose(3+5, 4))
# Pour vérification :
(proba <- dhyper(x = 3,m = 3, n = 5, k = 4))
#------------------