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"|] *)