1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114 | ------------------------------------------------------------------------------
-- File: Sim_Alea.adb
-- Description: Simulation of random variables with ASCII-Art histograms
-- Date/version: 1-Feb-2005; 17-Apr-1997
-- Author: Gautier de Montmollin
------------------------------------------------------------------------------
with Ada.Text_IO; use Ada.Text_IO;
with Ada.Float_Text_IO; use Ada.Float_Text_IO;
with Ada.Numerics.elementary_functions; use Ada.Numerics.elementary_functions;
with F_Random; use F_Random;
procedure Sim_Alea is
use Ada.Numerics;
-- g: inverse de la fonction de repartition F de la v.a. X
function g_unif(y: Float) return Float is -- X v.a. uniforme: X ~ U(a;b)
a: constant:= -10.0;
b: constant:= 50.0;
begin
return a + y * (b-a);
end;
function g_expon(y: Float) return Float is -- X v.a. expon.: X ~ Exp(l)
l: constant Float:= 1.0/9.0;
begin
return -Log(1.0-y)/l;
end;
function g_weibull(y: Float) return Float is -- X v.a. de Weib.: X ~ W(l;p)
p: constant Float:= 2.0; -- X=T^p, T ~ Exp(l)
l: constant Float:= 1.0/9.0;
log_l: constant Float:= Log(l);
begin
return exp( p* (Log(-Log(1.0-y)) - log_l ) );
end;
function g_cauchy(y: Float) return Float is -- X v.a. de Cauchy
begin -- X=tan(T), T ~ U(-pi/2;pi/2)
return Tan( pi * (y-0.5) );
end;
type Fct is access function(x: Float) return Float;
procedure Tirages(g: Fct; l: String; n: Positive; b_inf,b_sup: Float) is
-- loi: choix de la loi;
-- n: nb de tirages
-- b_inf,b_sup: bornes de l'histogramme
t,inf,sup,moy,somme: Float;
mx: constant:= 79; mxm1: constant:= mx-1;
my: constant:= 20;
ti: array(0..mxm1) of Integer:= (others=> 0); -- tableau d'incidences
lbarre: constant Float:= (b_sup-b_inf) / Float(mx);-- largeur d'une barre
ilb: constant Float:= 1.0/lbarre;
k: Integer;
procedure Histogramme is
corr: constant Integer:= n/1000; -- corr. pour aff.
m: constant Integer:= Integer( (moy-b_inf)*ilb );
begin
for i in ti'Range loop -- correction
ti(i):= ti(i) / corr;
end loop;
Put_Line("Loi " & l);
for j in reverse 1..my loop
for i in ti'Range loop
if ti(i)>=j then if i=m then Put('|'); else Put('#'); end if;
else Put(' '); end if;
end loop;
New_Line;
end loop;
Put(b_inf,6,2,0); Put(b_sup,mx-24,2,0); New_Line;
Put_Line("Stat: inf, moy, sup:");
Put(inf,6,4,0); Put(moy,30,4,0); Put(sup,mx-59,4,0); New_Line;
Skip_Line;
end;
begin
Randomize; -- initialiser le generateur selon l'horloge
-- Tirages de la v.a.
-- avec calcul de la somme, de la moyenne et de l'etendue des valeurs
somme:= 0.0;
for i in 1..n loop
t:= g(Random);
somme:= somme + t;
if i=1 then
sup:= t;
inf:= t;
else
if t>sup then sup:= t; end if;
if t<inf then inf:= t; end if;
end if;
k:= Integer( (t-b_inf)*ilb );
if k in ti'Range then
ti(k):= ti(k)+1;
end if;
end loop;
moy:= somme / Float(n);
Histogramme;
end Tirages;
begin
Tirages( g_unif'Access, "uniforme", 50_000, -20.0, 60.0 );
Tirages( g_weibull'Access, "de Weibull", 300_000, 0.0, 324.0 );
Tirages( g_cauchy'Access, "de Cauchy", 500_000, -300.0, 300.0 );
end Sim_Alea;
|