Simulations de lois de probabilité avec CAML

Simulation de la loi binomiale avec Caml et comparaison avec les valeurs théoriques. Calcul récursif des valeurs de la loi de Poisson. Diverses expériences aléatoires. Lois géométriques et hypergéométriques.

OCaml
(*----
    Outils généraux  
----*)
 
 
 
(* une liste d'expériences aléatoires *)
let simul_exp (exp: unit -> 'issue) (occ: int): 'issue list =
  let rec aux_exp = fun n acc ->
    if n = occ then acc
    else aux_exp (n + 1) ((exp ()) :: acc)
  in aux_exp 0 [];;
 
 
(* pour clarifier : on définit un type polymorphe pour les
   dictionnaire avec un constructeur d'arité 1  *)
type ('cle,'valeur) dictionnaire = Dict of ('cle * 'valeur) list;;
 
(* extrait de la doc de la bibliothèque List:
 
val assoc : 'a -> ('a * 'b) list -> 'b
 
assoc a l returns the value associated with key a in the list of pairs l. That is, assoc a [ ...; (a,b); ...] = b if (a,b) is the leftmost binding of a in list l. Raise Not_found if there is no value associated with a in the list l.
 
val assq : 'a -> ('a * 'b) list -> 'b
 
Same as List.assoc, but uses physical equality instead of structural equality to compare keys.
 
val mem_assoc : 'a -> ('a * 'b) list -> bool
 
Same as List.assoc, but simply return true if a binding exists, and false if no bindings exist for the given key.
 
val mem_assq : 'a -> ('a * 'b) list -> bool
 
Same as List.mem_assoc, but uses physical equality instead of structural equality to compare keys.
 
val remove_assoc : 'a -> ('a * 'b) list -> ('a * 'b) list
 
remove_assoc a l returns the list of pairs l without the first pair with key a, if any. Not tail-recursive.
 
val remove_assq : 'a -> ('a * 'b) list -> ('a * 'b) list
 
Same as List.remove_assoc, but uses physical equality instead of structural equality to compare keys. Not tail-recursive.
 
*)
 
 
(* fonction qui renvoie la valeur associée à la clé d'un dictionnaire*)
let valeur (cle: 'key) (dic: ('key,'value) dictionnaire): 'value =
    List.assoc cle (list_of_dic dic);;
 
(* fonction qui renvoie la liste associée au dictionnaire *)
let list_of_dic (d: ('cle,'valeur) dictionnaire): ('cle * 'valeur) list =
  match d with Dict(li) -> li;;
 
(* renvoie la liste des couples (élément, fréquence en pourcentage) *)
let compte (liste: 'issue list) : ('issue,float) dictionnaire =
  let effectif_total = float_of_int (List.length liste) in
  let rec count (l: 'issue list) (accu: ('issue * float) list) : ('issue,float) dictionnaire =
    match l with
    | []    -> Dict(accu)
    |t :: q ->
      let c = List.partition (fun x -> x = t) l in
      count (snd c) ((t,100. *. (float_of_int (List.length (fst c))) /. effectif_total)::accu)
  in count liste [];;
 
(* on entre l'expérience et le nb de simulations et on obtient le dictionnaire des
   fréquences en pourcentage *)
let simul_dic (exp_alea: unit -> 'issue) (nb_exp: int): ('issue,float) dictionnaire =
  let li = list_of_dic (compte (simul_exp exp_alea nb_exp)) in 
  Dict(List.sort compare li);;
 
(* idem avec un joli affichage sous forme de chaîne de caractères *)
let pretty_dic (ft: ('i -> 's, unit, string) format) (dic: ('a,float) dictionnaire): ('a , string) dictionnaire =
  Dict(List.map (fun cpl -> (fst cpl, Printf.sprintf ft (snd cpl) ^ " %")) (list_of_dic dic));;
 
 
(*
 Duc de Toscane 
*)
 
 
(* vecteur aléatoire de n entiers entre min et max *)
let vec_alea (min: int) (max: int) (taille: int): int list =
  let rec aux_alea (size: int) (acc: int list): int list =
    if size = 0 then acc
    else aux_alea (size - 1) ((min + Random.int (max - min + 1)) :: acc)
  in aux_alea taille [];;
 
(* sommes des éléments d'une liste d'entiers *)
let somme (liste : int list): int =
  List.fold_left (+) 0 liste;;
 
 
(* simulations du lancer de trois dés : X = somme des faces  *)
let toscane ((): unit): int = 
  somme (vec_alea 1 6 3);; 
 
 
(* simule 100_000 expériences du Duc de Toscane *)
simul_dic toscane 100_000;;
 
pretty_dic "%.2f" (simul_dic toscane 100_000);;
 
(*
  lancers d'une pièce
*)
type piece = Pile | Face ;;
 
let pile_ou_face ((): unit): piece =
  if Random.bool () then Pile else Face ;; 
 
simul_dic pile_ou_face 100_000;;
 
(*
  lancers d'un dé
*)
type face = I | II | III | IV | V | VI ;;
 
let de: (int,face) dictionnaire = Dict([(1,I);(2,II);(3,III);(4,IV);(5,V);(6,VI)]);;
 
let lancer_de ((): unit): face =
  valeur (1 + Random.int 6) de;;
  
simul_dic lancer_de 100_000;;
 
 
(* SURCHARGΕ EΝ CAML cf <a href="http://vernoux.fr/2011/11/17/la-surcharge-des-operateurs-en-objective-caml/" rel="nofollow">http://vernoux.fr/2011/11/17/la-surcharge-des-operateurs-en-objective-caml/</a> *)
 
(* obtenir le typage *)
let typage x =
    let tag_x = Obj.tag (Obj.repr x) in
        if tag_x = Obj.int_tag then "int"
        else if tag_x = Obj.string_tag then "string"
        else if tag_x = Obj.double_tag then "float"
        else "unknown";;
 
(* conversion vers le type flottant *)
let to_float x =
  let x_untyped = Obj.magic x
  and type_x = typage x in
  match type_x with
  | "int"   -> float (Obj.magic x_untyped)
  | "float" -> Obj.magic x_untyped
  | _ -> failwith "Je ne sais pas caster ça en float !";;
 
(* espérance quand 'a est float ou int *)
 
let esperance (dic: ('a , float) dictionnaire): float =
  List.fold_left (fun acc cpl -> acc +. ((to_float (fst cpl)) *. snd cpl))  0.  (list_of_dic dic);;
 
(* variance avec K.H. *)
let variance (dic: ('a , float) dictionnaire): float =
  let li = list_of_dic dic in
  let e = esperance dic in
  (List.fold_left (fun acc cpl -> acc +. (((to_float (fst cpl)) ** 2.) *. snd cpl))  0.  li) -. (e ** 2.);;
 
 
 
 
(*---
 
Problème du lièvre et de la tortue (poly ex 2.37 )
 
  ---*)
 
type la_fontaine = Lievre | Tortue;;
 
let lievre_tortue (nb_tirages_max: int) ((): unit) : la_fontaine =
  let rec jeu (nb_tirages : int): la_fontaine =
    if nb_tirages = nb_tirages_max then
      Tortue
    else 
      let de = 1 + Random.int 6 in
      if de == 6 then
        Lievre
      else
        jeu (nb_tirages + 1)
  in jeu 0;;
 
(*---
  
  LOI BINOMIALE
  exo 2.36 : 5 jetons tirés successivement dans une urne avec 5 blancs et 9 noirs
  X = nb de blancs tirés.
  
--*)
 
(* renvoie 1 avec une probabilité p et 0 sinon *)
exception Pas_une_proba ;;
 
let bernoulli (p: float) ((): unit): int  =
  if ((p > 1.) || (p < 0.)) then raise Pas_une_proba
  else 
    if (Random.float 1. < p) then
      1
    else 0;;
 
(*  ou bien :
    let bernoulli (p: float) : unit -> int  =
       fun () -> if (Random.float 1. < p) then 1 else 0;;
*)
 
(* X ~ B(n,p) comme somme de var de Bernoulli B(1,p) *)
let binomial  (n: int) (p: float) (():  unit): int  =
  somme (simul_exp (bernoulli p) n);;
 
(*
ou bien :
  let binomial  (n: int) (p: float) :  unit -> int  =
     fun () -> som (simul_exp (bernoulli p) n);;
 
*)
 
 
simul_dic (binomial 5 (5. /. 14.)) 100_000;;
 
esperance (simul_dic (binomial 5 (5. /. 14.)) 100_000) ;;
 
(* Comparaison avec le calcul théorique *)
 
(* on utilise C(n,p) = (n-p+1)/p * C(n,p-1) -> acc_n = (num -1)*acc_{n-1}/(den + 1) et acc_0 = 1 *)
let binom (n: int) (p: int) : int =
   if p > n || n < 0 || p < 0 then
     0
   else
     let rec aux = fun acc num den ->
       if den <= p then
         aux ((acc * num) / den) (num - 1) (den + 1)
       else
         acc
     in aux 1 n 1;;
 
(* calcul de Px(k) quand X ~ B(n,p) *)
let proba_bin (n: int) (p: float) (k: int) : float =
  (float_of_int (binom n k)) *. (p ** (float_of_int k)) *. ((1. -. p)**(float_of_int (n - k)));;
 
(* renvoie le dictionnaire des (k, Px(k)) quand X ~ B(n,p) *)
let dic_bin (n: int) (p: float) : (int , float) dictionnaire =
  let rec aux = fun i acc ->
    if i > n then
      Dict(List.sort compare acc)
    else
      let pr = proba_bin n p i in
      aux (i + 1) ((i,pr) :: acc)
  in aux 0 [];;
 
dic_bin 5 (5. /. 14.);; 
 
esperance (dic_bin 5 (5. /. 14.));;
 
5. *. 5. /. 14.;;
 
 
 
 
 
 
(*----
 
  LOI DE POISSON
 
  ---*)
 
(* on utilise Px(k + 1) = (lambda / (k + 1)) * Px(k) *)
let poisson (lambda: float) (n: int) : float =
  let rec aux = fun k acc ->
    if k = n then acc
    else aux (k + 1) ((lambda /. (float_of_int (k + 1))) *. acc)
  in aux 0 (exp (-. lambda));;
 
(* fabrication de la table sous forme d'un tableau de chaînes de caractères *)
 
(* un tableau de longueur 10 rempli de 0 *)
let rang : int array =
  Array.make 10 0;;
 
(* les valeurs de lambda : tableau de 0.1 à 1.0 par pas de 0.1 *)
let rang_lambda : float array =
  Array.mapi (fun j x ->  (float_of_int (j + 1)) /. 10.) rang;;
 
(* une matrice constituée de 8 lignes identiques = aux valeurs de lambda *)
let tab : float array array =
  Array.make 8 rang_lambda;;
 
(* la table de Poisson brute : on applique la fonction poisson aux éléments de tab *)
let tab_poisson0 : string array array =
  Array.mapi (fun i r -> Array.map (fun la -> Printf.sprintf "%.3f" (poisson la i)) r) tab;;
 
(* La table de Poisson bordée des valeurs de k et de lambda *)
let tab_poisson : string array array =
  let r0 =
    Array.append [| "k / la" |] (Array.map (fun x -> " " ^ Printf.sprintf "%.1f" x^ " ") rang_lambda) in
  let t0 =
    Array.mapi (fun i r -> Array.append [| "  " ^ Printf.sprintf "%.1i" i ^ "   "|] r) tab_poisson0 in
 Array.append [|r0|] t0;;
 
 
(*
 
  [|"k / la"; " 0.1 "; " 0.2 "; " 0.3 "; " 0.4 "; " 0.5 "; " 0.6 "; " 0.7 "; " 0.8 "; " 0.9 "; " 1.0 "|];
  [|"  0   "; "0.905"; "0.819"; "0.741"; "0.670"; "0.607"; "0.549"; "0.497"; "0.449"; "0.407"; "0.368"|];
  [|"  1   "; "0.090"; "0.164"; "0.222"; "0.268"; "0.303"; "0.329"; "0.348"; "0.359"; "0.366"; "0.368"|];
  [|"  2   "; "0.005"; "0.016"; "0.033"; "0.054"; "0.076"; "0.099"; "0.122"; "0.144"; "0.165"; "0.184"|];
  [|"  3   "; "0.000"; "0.001"; "0.003"; "0.007"; "0.013"; "0.020"; "0.028"; "0.038"; "0.049"; "0.061"|];
  [|"  4   "; "0.000"; "0.000"; "0.000"; "0.001"; "0.002"; "0.003"; "0.005"; "0.008"; "0.011"; "0.015"|];
  [|"  5   "; "0.000"; "0.000"; "0.000"; "0.000"; "0.000"; "0.000"; "0.001"; "0.001"; "0.002"; "0.003"|];
  [|"  6   "; "0.000"; "0.000"; "0.000"; "0.000"; "0.000"; "0.000"; "0.000"; "0.000"; "0.000"; "0.001"|];
  [|"  7   "; "0.000"; "0.000"; "0.000"; "0.000"; "0.000"; "0.000"; "0.000"; "0.000"; "0.000"; "0.000"|]
 
 
*)

courtesy of webmatter.de