Simulations...

Ca fait un moment que je n'ai pas posté de créations originales alors j'en propose une. Comme la rubrique "Statistique" n'est pas la plus populaire, je me suis dit que ca pourrait être bien de la remplir un peu...

Mon fichier consiste en la simulation de variables aléatoires :

• Loi uniforme sur [0, 1]
  
• Loi de Cauchy standard
  
• Loi normale centrée réduite

Pour la simulation de la loi uniforme, pas de problème : il suffit d'utiliser la fonction Rand().

Pour la simulation de la loi de Cauchy, j'ai utilisé la méthode de l'inverse de la fonction de répartition : c'est très simple a mettre en oeuvre!

En ce qui concerne la loi normale, j'ai utilisé la méthode de Box-Muller. Là encore, c'est très simple à mettre en place.

Je vous donne les résultats obtenus avec mon fichier :



Et voici le fichier source (à améliorer!!) :

Code:

% TeXgraph version 1.93
% Fenetre Xmin Xmax Ymin Ymax Xscale Yscale
100#-4#4#-0.05#0.5#1.5#20##
% Marges gauche droite haut bas cadre gestion_couleur comptgraph
101#0#0#0#0#0#1#7##
% Affectation des Variables theta et phi
18##[theta:=0.5236, phi:=1.0472, OriginalCoord(1)]##
% Déclaration des Variables Globales
15#A#SimNorm(1000)##
% Déclaration des Macros
16#Hist#[
$L:=%1, $h:=%2,
$m:=min(L), $M:=max(L),
Seq(sum(Map(Si($x<=z And z<x+h, 1), z, L)), x, m, M, h)
]
##
16#DrawHist#[
$L:=%1, $h:=%2,
$H:=Hist(L, h)/(Nops(L)*h),
$k:=1,
Ligne(Seq([rect(x, x+i*H[k], Si(x+h<max(L), x+h, max(L))), jump, Inc(k, 1)],
                  x, max(min(L), Xmin), min(Xmax, max(L)), h), 1)
]
##
16#SimCauchy#Seq(tan(pi*(Rand()-1/2)), $x, 1, %1)
##
16#dnorm#(1/sqrt(2*pi))*exp((-%1^2)/2)
##
16#dcauchy#1/(pi*(1+%1^2))
##
16#DrawHist2#[
$L:=%1, $h:=%2,
$m:=min(L), $M:=max(L),
A:=Seq(sum(Map(Si($x<=z And z<x+h, 1), z, L)), x, max(m, Xmax), min(M, Xmax), h)  ,
B:=Seq($x, x, max(m, Xmin), min(M, Xmax), h),
Ligne(Seq(rect(B[k], B[k]+i*A[k], B[k+1]), k, 1, Nops(A)-1))
]
##
16#SimNorm#Seq(sqrt(-2*ln(Rand()))*cos(2*pi*Rand()), $x, 1, %1)
##
% Déclaration des Eléments graphiques
% SimUnif (Utilisateur)
18##[IsVisible:=0]##
14#SimUnif#[
$N:=10000,
$h:=0.25,
$X:=Seq(Rand(), $x, 1, N),
FillStyle:=full, FillColor:=Light(gray, 0.3),
DrawHist(X, h),
FillStyle:=none, Width:=8, Color:=darkred,
Seg(i, 1+i)
]
#-1##
% SimCauchy (Utilisateur)
14#SimCauchy#[
$N:=10000,
$h:=0.25,
$X:=SimCauchy(N),
FillStyle:=full, FillColor:=Light(gray, 0.3),
DrawHist(X, h),
FillStyle:=none, Width:=8, Color:=darkred,
Courbe(t+i*\dcauchy(t))
]
#-1##
% SimNorm (Utilisateur)
18##[IsVisible:=1,ForMinToMax:=1]##
14#SimNorm#[
$N:=10000,
$h:=0.25,
$X:=SimNorm(N),
FillStyle:=full, FillColor:=Light(gray, 0.3),
DrawHist(X, h),
FillStyle:=none, Width:=8, Color:=darkred,
Courbe(t+i*\dnorm(t))
]
#-1##
% Axes (Axes)
18##[xyticks:=0.1,LabelSize:=1]##
1#Axes#[0+i*(0),1+i*(0.1),2+i*(2)]#0#0##


A noter que la simulation de la loi de Cauchy peut prendre pas mal de temps avec TeXgraph. En fait, cette loi prend trop souvent de grande valeurs (elle n'a même pas d'espérance!!) ce qui fait sévèrement ramer ma macro DrawHist.
Je pense que c'est d'ailleurs la macro à modifier dans mon fichier. Je l'ai fait en 2 min donc elle est loin d'être parfaite! (Ya même un soucis avec ; si je ne me trompe pas, elle trace pas tout a fait ce qu'elle devrait, mais bon : elle fait quand même bien son boulot!)

Salut Alphonse,

J'ai regardé ton fichier, et effectivement l'exécution est lente. Alors je te propose ces modifications qui améliorent la vitesse, l'astuce est de trier la liste de valeurs avant de calculer les fréquences, ce calcul se fait alors avec un seul parcourt de liste, ainsi que la construction des rectangles, ce qui fait que l'affichage est quasiment instantané maintenant:

Code:


% TeXgraph version 1.94 beta-7.3
% Fenetre Xmin Xmax Ymin Ymax Xscale Yscale
100#-4#4#-0.05#0.45#1.5#24##
% Marges gauche droite haut bas cadre gestion_couleur comptgraph
101#0#0#0#0#0#1#7##
% Affectation des Variables theta et phi et type de perspective
18##[theta:=0.5236, phi:=1.0472,OriginalCoord(1),IdMatrix(),IdMatrix3D(),
ModelView(ortho)]##
% Déclaration des Macros
16#Hist#[
$L:=%1, $h:=%2, $x:=L[1]+h, $val:=0,
for $z in L do
    if z<x then Inc(val,1)
          else val, val:=1, Inc(x,h)
    fi
od,
val
]##
16#DrawHist#[
$L:=%1, $h:=%2, $H:=Hist(L, h)/(Nops(L)*h),
$x:=L[1],
Ligne (for y in H do x,x+h,x+h+i*y,x+i*y,jump, Inc(x,h) od,1)
]##
16#SimCauchy#[$L:=for $k from 1 to %1 do
        repeat
              $x:=tan(pi*(Rand()-1/2))
        until (-5<=x) And (x<=5) od,
        x
    od,
 Sort(L),
 L]##
16#dnorm#(1/sqrt(2*pi))*exp((-%1^2)/2)##
16#dcauchy#1/(pi*(1+%1^2))##
16#SimNorm#[$L:=for $k from 1 to %1 do
        sqrt(-2*ln(Rand()))*cos(2*pi*Rand())
    od,
 Sort(L), L]##
% Déclaration des Eléments graphiques
% SimUnif (Utilisateur)
18##[IsVisible:=0]##
14#SimUnif#[  {sur [0,1]}
view(-0.05,1.1,-0.05,1.25), size(12,1),
$N:=10000,
$h:=0.01,
$X:=Seq(Rand(), $x, 1, N), Sort(X),
FillStyle:=full, FillColor:=Light(gray, 0.3),
DrawHist(X, h),
FillStyle:=none, Width:=8, Color:=darkred,
Seg(i, 1+i),
xyticks:=0.1, LabelSize:=scriptsize,
Color:=black, Axes(0,0.1+0.1*i)
]#
-1##
% SimCauchy (Utilisateur)
14#SimCauchy#[{sur [-4,4]}
view(-4,4,-0.05,0.4), size(12,0),
$N:=10000,
$h:=0.25,
$X:=SimCauchy(N),
FillStyle:=full, FillColor:=Light(gray, 0.3),
DrawHist(X, h),
FillStyle:=none, Width:=8, Color:=darkred,
Courbe(t+i*\dcauchy(t)),
xyticks:=0.1, LabelSize:=scriptsize,
Color:=black, Axes(0,1+0.1*i)
]#
-1##
% SimNorm (Utilisateur)
18##[IsVisible:=1,ForMinToMax:=1]##
14#SimNorm#[{sur [-4,4]}
view(-4,4,-0.05,0.45), size(12,0),
$N:=10000,
$h:=0.25,
$X:=SimNorm(N),
FillStyle:=full, FillColor:=Light(gray, 0.3),
DrawHist(X, h),
FillStyle:=none, Width:=8, Color:=darkred,
Courbe(t+i*\dnorm(t)),
xyticks:=0.1, LabelSize:=scriptsize,
Color:=black, Axes(0,1+0.1*i)
]#
-1##


Pour la simulation de Cauchy, j'ai forcé le tirage jusqu'à l'obtention d'une valeur entre -5 et 5, qu'en penses-tu? J'ai quand même l'impression que la simulation n'est pas bonne, j'ai testé avec N=100000 et ça délire!

Re,

Ca délire à partir de N=65536, mais c'est un bug de la commande Nops (elle travaille avec un type entier trop "court"!). Ce sera rectifié dans la prochaine version.

Salut Patrick!!!

J'ai regardé vite fait ce que tu as fait, mais j'ai pas encore testé. Je vais le faire dans les 10min qui viennent...

P.Fradin a écrit:Pour la simulation de Cauchy, j'ai forcé le tirage jusqu'à l'obtention d'une valeur entre -5 et 5, qu'en penses-tu? J'ai quand même l'impression que la simulation n'est pas bonne, j'ai testé avec N=100000 et ça délire!

Le fait d'imposer que la simulation soit comprise entr -5 et 5 doit fausser la simulation je pense. On a en fait affaire à une espèce de loi de Cauchy tronquée des deux côtés. Je pense que les fréquences sont alors modifiées et légèrement supérieures à ce qu'elles devraient être.

Je pense que cette sélection des valeurs de la simulation ne doit pas se faire dans la fonction SimCauchy, mais plutôt dans l'élément graphique permettant de dessiner l'histogramme (c'est ce que j'avais tenté de faire, mais c'était très lent) En gros, je suggère que SimCauchy renvoie toutes les valeurs simulées (même si elles sont trop grandes) et ensuite, on fait la sélection des valeurs qui nous intéressent dans l'élément graphique utilisateur, et pour les fréquences, on tient compte des valeurs simulées supérieures a 5 en valeur absolue.

Je sais pas si j'ai été très clair!!!

Alphonse Capriani a écrit:
Le fait d'imposer que la simulation soit comprise entr -5 et 5 doit fausser la simulation je pense. On a en fait affaire à une espèce de loi de Cauchy tronquée des deux côtés. Je pense que les fréquences sont alors modifiées et légèrement supérieures à ce qu'elles devraient être.

Oui très certainement, j'ai essayé en élargisant la fourchette à [-100,100] et c'est mieux.

Je pense que cette sélection des valeurs de la simulation ne doit pas se faire dans la fonction SimCauchy, mais plutôt dans l'élément graphique permettant de dessiner l'histogramme (c'est ce que j'avais tenté de faire, mais c'était très lent) En gros, je suggère que SimCauchy renvoie toutes les valeurs simulées (même si elles sont trop grandes) et ensuite, on fait la sélection des valeurs qui nous intéressent dans l'élément graphique utilisateur, et pour les fréquences, on tient compte des valeurs simulées supérieures a 5 en valeur absolue.

Je sais pas si j'ai été très clair!!!

Si, très clair. Je vais voir ce que l'on peut faire.

En effet, je viens de tester tes nouvelles macro et avec la simulation de la loi de Cauchy, on a clairement un soucis : fréquences trop élevées!

Cela dit, Tes macros sont incroyablement plus rapides que les miennes : félicitations!

Sinon, j'avais juste une petite question à propose de la fonction Rand().

Comment elle fait pour simuler suivant la loi uniforme? La première chose à laquelle j'ai pensé, c'est aux générateurs congruentiels, mais j'ai regardé et ca n'a pas l'air d'en être un. C'est un générateur de type registre? Un mix entre registre et congruentiel? Ou autre chose?

---

Concernant le problème de la loi de Cauchy avec 100 000 simulation, à mon avis, ca vient pas de texgraph mais de ce problème de fréquences qui sont réhaussées...

Poursuivons les simulations de variables aléatoires avec un joli mélange de gaussiennes :



Et voici le code du fichier Simulation.teg modifié :

Code:

% TeXgraph version 1.94 beta-7.3.5
% Fenetre Xmin Xmax Ymin Ymax Xscale Yscale
100#-3#9#-0.05#0.25#1#40##
% Marges gauche droite haut bas cadre gestion_couleur comptgraph
101#0#0#0#0#0#1#8##
% Affectation des Variables theta et phi et type de perspective
18##[theta:=0.5236, phi:=1.0472,OriginalCoord(1),IdMatrix(),IdMatrix3D(),
ModelView(ortho)]##
% Déclaration des Macros
16#Hist#[
$L:=%1, $h:=%2, $x:=L[1]+h, $val:=0,
for $z in L do
    if z<x then Inc(val,1)
          else val, val:=1, Inc(x,h)
    fi
od,
val
]##
16#DrawHist#[
$L:=%1, $h:=%2, $H:=Hist(L, h)/(Nops(L)*h),
$x:=L[1],
Ligne (for y in H do x,x+h,x+h+i*y,x+i*y,jump, Inc(x,h) od,1)
]##
16#SimCauchy#[$L:=for $k from 1 to %1 do
        repeat
              $x:=tan(pi*(Rand()-1/2))
        until (-5<=x) And (x<=5) od,
        x
    od,
 Sort(L),
 L]##
16#dnorm#[
$m:=Si(nil(%2), 0, %2), $s:=Si(nil(%3), 1, %3),
(1/sqrt(2*pi*s^2))*exp(-((%1-m)^2)/(2*s^2))
]##
16#dcauchy#1/(pi*(1+%1^2))##
16#SimNorm#[
$m:=Si(nil(%2), 0, %2), $s:=Si(nil(%3), 1, %3),
$L:=for $k from 1 to %1 do
        m+s*sqrt(-2*ln(Rand()))*cos(2*pi*Rand())
    od,
 Sort(L), L
 ]##
16#SimMelNorm#[
$N:=%1, $p:=%2,
$m1:=%3, $s1:=%4,
$m2:=%5, $s2:=%6,
$B:=sum(SimBernoulli(N, p)),
$L:=[SimNorm(B, m1, s1), SimNorm(N-B, m2, s2)],
Sort(L), L
]##
16#SimBernoulli#Seq(Si(Rand()<%2, 1, 0), $x, 1, %1)##
% Déclaration des Eléments graphiques
% SimUnif (Utilisateur)
18##[IsVisible:=0,tMin:=-4,tMax:=4]##
14#SimUnif#[  {sur [0,1]}
view(-0.05,1.1,-0.05,1.25), size(12,1),
$N:=10000,
$h:=0.01,
$X:=Seq(Rand(), $x, 1, N), Sort(X),
FillStyle:=full, FillColor:=Light(gray, 0.3),
DrawHist(X, h),
FillStyle:=none, Width:=8, Color:=darkred,
Seg(i, 1+i),
xyticks:=0.1, LabelSize:=scriptsize,
Color:=black, Axes(0,0.1+0.1*i)
]#
-1##
% SimCauchy (Utilisateur)
14#SimCauchy#[{sur [-4,4]}
view(-4,4,-0.05,0.4), size(12,0),
$N:=100000,
$h:=0.25,
$X:=SimCauchy(N),
FillStyle:=full, FillColor:=Light(gray, 0.3),
DrawHist(X, h),
FillStyle:=none, Width:=8, Color:=darkred,
Courbe(t+i*\dcauchy(t)),
xyticks:=0.1, LabelSize:=scriptsize,
Color:=black, Axes(0,1+0.1*i)
]#
-1##
% SimNorm (Utilisateur)
18##[ForMinToMax:=1,tMax:=8]##
14#SimNorm#[{sur [-4,4]}
{view(-4,4,-0.05,0.45), size(12,0),  }
$N:=10000,
$h:=0.25,
$m1:=0, s1:=1,
$m2:=4, s2:=1.5,
$X:=SimNorm(N, m1, s1),
$Y:=SimNorm(N, m2, s2),
FillStyle:=fdiag, FillColor:=Light(gray, 0.3),
DrawHist(X, h),
FillStyle:=bdiag, FillColor:=green,
DrawHist(Y, h),
FillStyle:=none, Width:=8, Color:=darkred,
Courbe(t+i*\dnorm(t, m1, s1)),
Courbe(t+i*\dnorm(t, m2, s2)),
xyticks:=0.1, LabelSize:=scriptsize,
Color:=black, Axes(0,1+0.1*i)
]#
-1##
% SimMelNorm (Utilisateur)
18##[IsVisible:=1,tMin:=-3,tMax:=9]##
14#SimMelNorm#[
view(-3,9,-0.05,0.25), size(12,0),
$N:=50000,
$h:=0.15,
$p:=0.25,
$m1:=0, $s1:=0.75, $m2:=4, $s2:=1.5,
$X:=SimMelNorm(N, p, m1, s1, m2, s2),
FillStyle:=full, FillColor:=Light(gray, 0.7),
DrawHist(X, h),
FillStyle:=none, Width:=8, Color:=darkred,
Courbe(t+i*(p*\dnorm(t, m1, s1)+(1-p)*\dnorm(t, m2, s2))),
xyticks:=0.1, LabelSize:=scriptsize,
Color:=black, Axes(0,1+0.1*i)
]#
-1##


Re,

Pour la fonction Rand(), je ne saurais te dire, elle est prédéfinie en Pascal. A priori, les tests montrent bien qu'elle est uniforme!

Pour Cauchy, si j'ai bien compris ta démarche, tu choisis un y dans [0,1] (loi uniforme) et tu renvoies x=F^(-1)(y). Plutôt que de filtrer x j'ai essaye en prenant y dans [0.01; 0.99] par exemple, ce qui donnerait pour la macro SimCauchy:

Code:


[$L:=for $k from 1 to %1 do  tan(pi*(0.98*Rand()+0.01-1/2))  od,
 Sort(L),  L]


J'obtiens ça avec N=500 et h=0.5:

avec N=50000 et h=0.25:

Oui! C'est une méthode qui donne toujours un résultat légèrement biaisé, mais c'est quand même mieux que ta première solution.

Mais ne serait il pas plus simple de sélectionner que les éléments compris entre -4 et 4 parmis les simulations? Ma question (que j'avais déjà formulé il me semble) amène plutot la suggestion suivante : il serait utile que TeXgraph soit muni d'une macro Select(<List>, <variable>, <Condition>) permettant de garder uniquement les complexes de la liste satisfaisant la condition mise en 3ème paramètre. (Finalement, c'est juste une variante de la macro Map) A mon avis, même avec ce genre de macro, l'affichage des simulations de la loi de Cauchy va prendre du temps car la sélection des éléments de la liste vérifiant la condition risque de prendre un peu de temps. Néanmoins, une telle fonction peut s'avérer utile.

Salut Alphonse,

Je ne sais pas si c'est l'âge, ou les bulles , mais cette variante de Map existe déjà, et devine qui me l'a suggérée?

Par contre il faut modifier la macro DrawHist pour ajouter l'effectif total N en paramètre (facultatif):

Code:


{DrawHist(L, h, N)}
[$L:=%1, $h:=%2, $N:=if %3=Nil then Nops(L) else %3 fi,
$H:=Hist(L, h)/(N*h),
$x:=L[1],
Ligne (for y in H do x,x+h,x+h+i*y,x+i*y,jump, Inc(x,h) od,1)
] 


Ne pas oublier de modifier la commande SimCauchy:

Code:


[{sur [-4,4]}
view(-4,4,-0.05,0.4), size(12,0),
$N:=100000,
$h:=0.25,
$X:=for z in SimCauchy(N) andif (-5<z) And (z<5) do z odfi,
FillStyle:=full, FillColor:=Light(gray, 0.3),
DrawHist(X, h, N),
FillStyle:=none, Width:=8, Color:=darkred,
Courbe(t+i*\dcauchy(t)),
xyticks:=0.1, LabelSize:=scriptsize,
Color:=black, Axes(0,1+0.1*i)
]


Tu vois où est la variante de Map ?

Ce qui me donne:

Salut tout le monde!

ET BONNE ANNEE 2009!!!

P.Fradin a écrit:Je ne sais pas si c'est l'âge, ou les bulles Laughing, mais cette variante de Map existe déjà, et devine qui me l'a suggérée?

Décidément, ca va pas fort!!! Vu que j'ai plus trop le temps de me pencher sur TeXgraph ces derniers mois, j'en oublie les suggestions et les tests que j'ai fait sur les version précédentes (et même sur la version 1.94). Il est temps que je trouve un peu de temps pour m'y reconsacrer...

Salut Alphonse et Bonne Année!

Pour te remonter le moral, je te fais remarquer que ta suggestion qui consiste à filtrer les valeurs après le tirage et non pendant, donne de bons résultats non?

Tout à fait!!!

Le résultat est nickel et bien plus rapide que ma première proposition de simulation!!!