Problème avec une simulation informatique

Guego

Bonsoir, je ne sais pas si je dois poster dans Probas ou dans Informatique, je tente dans Probas.

J'ai un couple de variables aléatoires compliqué $(X,Y)$ (pas indépendantes) à valeurs dans $]-\pi;\pi]^2$, que je sais simuler informatiquement (j'ai ma fonction $\texttt{SimulXY()}$ en Python qui a l'air de marcher). Je m'intéresse à l'espérance (conditionnelle) de $\cos(X)+\cos(Y)+\cos(X+Y)$ sachant $\sin(X)+\sin(Y)=\sin(X+Y)$.

Pour ça, j'ai fait le programme suivant :

T0=[]
n=0
while n<1000:
    (x,y)=SimulXY()
    if np.sin(x)+np.sin(y)==np.sin(x+y):
        n+=1
        T0.append(np.cos(x)+np.cos(y)+np.cos(x+y))

print(sum(T0)/1000)

J'obtiens environ 0.5.

Ensuite, je réfléchis, et je me dis que $\sin(X)+\sin(Y)=\sin(X+Y)$ équivaut à $X=0$ ou $Y=0$ ou $Y=-X$. De plus, étant donné la tête de mes variables aléatoires, ces 3 cas sont incompatibles, donc je me dis que je vais regarder ce qu'il se passe dans chacun de ces 3 sous-cas :

T1=[]
T2=[]
T3=[]
n=0
while n<1000:
    (x,y)=SimulXY()
    if np.sin(x)+np.sin(y)==np.sin(x+y):
        n+=1
        u=np.cos(x)+np.cos(y)+np.cos(x+y)
        if x==0:
            T1.append(u)
       if y==0:
            T2.append(u)
        if (x+y)==0:
            T3.append(u)

Je constate alors (en regardant la taille de mes 3 listes) que les $3$ sous-cas sont équiprobables et que j'obtiens environ 0.5 comme moyenne de chacune de mes 3 listes. Super. Jusqu'ici tout va bien.

Un peu plus tard, je ne retrouve plus mon fichier avec mes programmes. Je décide de re-simuler le cas $X=0$. Étant donné que $X=0$ implique $\sin(X)+\sin(Y)=\sin(X+Y)$, je n'ai pas besoin de mettre le $\texttt{if np.sin(x)+np.sin(y)==np.sin(x+y)}$. Je fais directement :

T4=[]
n=0
while n<1000:
    (x,y)=SimulXY()
    if x==0:
        n+=1
        T4.append(np.cos(x)+np.cos(y)+np.cos(x+y))

print(sum(T4)/1000)

Et là, boum ! J'obtiens 0.2, et donc pas du tout 0.5, comme précédemment ! Je me dis qu'il y a peut-être des fluctuations aléatoires dues aux simulations, que je n'en ai pas fait assez, etc. J'ai recommencé plusieurs fois, j'ai multiplié le nombre de simulations par 10. Systématiquement, avec T0, T1, T2, T3, j'obtiens une moyenne proche de 0.5, et systématiquement, avec T4, j'obtiens une moyenne proche de 0.2.

Je ne comprends vraiment pas pourquoi et ça me rend fou. Est-ce que quelqu'un aurait une explication ?

MrJ

Je n'arrive pas non plus à voir le problème en survolant les codes.

Pour comprendre ce qu'il se passe, je serai tenté de fabriquer toutes les listes dans une seule boucle while (et donc avec les mêmes simulations) pour les comparer à la fin.

T0 = []
T1 = []
T2 = []
T3 = []
T4 = []
n=0
while n<1000:
    (x,y)=SimulXY()
    if np.sin(x)+np.sin(y)==np.sin(x+y):
        n+=1        
        u=np.cos(x)+np.cos(y)+np.cos(x+y)
        T0.append(u)

        if x==0:
            T1.append(u)
        if y==0:
            T2.append(u)
        if (x+y)==0:
            T3.append(u)
    if x==0:
        T4.append(np.cos(x)+np.cos(y)+np.cos(x+y))

Édit : Dans une première réponse, j'avais mentionné la possibilité que l'erreur puisse venir de la présence de tests d'égalité avec des nombres flottants dont je ne suis pas sûr du comportement en Python.

Chalk

Sans parler du problème de l'égalité de nombres flottants soulevée par MrJ, le tout premier code me surprend. Sauf cas très particuliers (et dépendance vraiment très forte) sur la loi (X,Y), il est tout à fait improbable d'obtenir une égalité parfaite np.sin(x)+np.sin(y)==np.sin(x+y) même sur des millions d'itérations (d'autant plus que de manière standard le flottant Python est sur 64 bits).

Ce serait à mon avis utile d'en savoir plus sur (X,Y). Es-tu sûr que les 3 cas évoqués sont incompatibles ? As-tu essayé avec des else if et pas juste des if dans ton deuxième code pour voir la différence ?

Guego

MrJ a dit :

J'avais mentionné la possibilité que l'erreur puisse venir de la présence de tests d'égalité avec des nombres flottants dont je ne suis pas sûr du comportement en Python.

Je pense que c'est ça le problème. Les programmes que j'ai postés sont des simplifications. En réalité, je ne teste pas vraiment une égalité, mais je teste en mettant "si la valeur absolue de la différence est plus petite qu'un epsilon". Je ne l'ai pas signalé parce que j'ai testé avec plusieurs valeurs de epsilon (voire différentes valeurs à différents endroits du programme) et que ça ne changeait rien, mais ça vient peut-être quand même de là.

Guego

Version légèrement modifiée du programme de MrJ + Vérification qu'on n'a pas de recoupement entre les 3 cas :

T0 = []
T1 = []
T2 = []
T3 = []
T4 = []
n=0
while n<10000:
    (x,y)=SimulXY()
    if abs(np.sin(x)+np.sin(y)-np.sin(x+y))<0.001:
        n+=1        
        u=np.cos(x)+np.cos(y)+np.cos(x+y)
        T0.append(u)
        v=min([abs(x),abs(y),abs(x+y)])
        c=0
        if abs(x)==v:
            T1.append(u)
            c+=1
        if abs(y)==v:
            T2.append(u)
            c+=2
        if abs(x+y)==v:
            T3.append(u)
            c+=4
        if c not in [1,2,4]:
            print("Avertissement")
            print([c,(x,y,x+y,np.sin(x)+np.sin(y)-np.sin(x+y))])
    if abs(x)<0.001:
        T4.append(np.cos(x)+np.cos(y)+np.cos(x+y))

Résultat : aucun avertissement. Résultats numériques :

>>> print(np.mean(T0))
0.48943685220073985
>>> print(np.mean(T1))
0.497922408071928
>>> print(np.mean(T2))
0.4924253674913608
>>> print(np.mean(T3))
0.47794841317695635
>>> print(np.mean(T4))
0.18962312055618538

lourrran

Dans le dernier code que tu postes, tu testes : if abs(np.sin(x)+np.sin(y)-np.sin(x+y))<0.001 
ok.
Et ensuite, tu splittes en différents cas, mais tu testes avec x , y et x+y, au lieu de sin(x), sin(y) et sin(x+y) ; si $x$ est proche de $\pi$, tu ne récupères pas x, mais y avec l'instruction v = min ( abs(x), abs(y), abs(x+y) ) .

Je pense que quand tu affiches tes print(np.mean(T0)), etc, tu devrais aussi afficher le nombre d'éléments de T0 ; tu y gagnerais en visibilité.

MrJ

Le problème ne viendrait-il pas du fait que $X$ et $Y$ peuvent être proches de $\pm \pi$ simultanément ? Cela rendrait l'égalité $\sin(X)+\sin(Y) = \sin(X+Y)$ "presque" vraie.

Calli

Bonjour,
Est-ce que (X,Y) est une va continue ou discrète ? Et si tu peux donner SimulXY, ça pourrait nous permettre de faire des petits tests.

Guego

@Calli : C'est une v.a continue. Code de SimulXY, si vous voulez tester (en tout cas, merci de vous intéresser à mon problème)

import numpy as np
rng=np.random.default_rng()

def SimulXY():
    (X,Y,Z)=(0,0,1)
    constante=1/(24*(np.pi**2))
    while Z>constante*(6 - 4*np.cos(X - Y) - 4*np.cos(X + 2*Y) - 4*np.cos(2*X + Y) + 4*np.cos(3*X) + 4*np.cos(3*Y) - 2*np.cos(2*X - 2*Y) + 4*np.cos(3*X + 3*Y) - 2*np.cos(2*X + 4*Y) - 2*np.cos(4*X + 2*Y)):
        (X,Y,Z)=(rng.random()*2*np.pi-np.pi, rng.random()*2*np.pi-np.pi,rng.random()*0.12)
    return (X,Y)

lourrran

Affiche la taille de T4. Tu vas constater que T4 a environ 5000 lignes quand T0 en a 3300.
En fait, $abs(sin(x) ) <0.001$ n'entraine pas  $abs(sin(x) +sin(y) - sin(x+y) ) < 0.001$
Modifions un peu le code, pour voir ces lignes qui sont 'attrapées' dans T4, mais pas dans T0.

Edit : x ou y peuvent être très proches de $\pi$ ; J'ai donc corrigé la ligne v= ... ; l'ancienne version est en commentaire.  Les lignes if abs() qui suivent sont également modifiées.  

T0 = []
T1 = []
T2 = []
T3 = []
T4 = []
n=0
while n<10000:
    (x,y)=SimulXY()
    c=0
    if abs(np.sin(x)+np.sin(y)-np.sin(x+y))<0.001:
        n+=1        
        u=np.cos(x)+np.cos(y)+np.cos(x+y)
        T0.append(u)
        #v=min([abs(x), abs(y), abs(x+y) ])
        v=min([abs(np.sin(x)), abs(np.sin(y)), abs(np.sin(x+y))])
        #c=0
        if abs( np.sin(x) )==v:
            T1.append(u)
            c+=1
        if abs(np.sin(y) )==v:
            T2.append(u)
            c+=2
        if abs( np.sin(x+y) )==v:
            T3.append(u)
            c+=4
        if c not in [1,2,4]:
            print("Avertissement")
            print([c,(x,y,x+y,np.sin(x)+np.sin(y)-np.sin(x+y))])
    if abs( np.sin(x) )<0.001:
        T4.append(np.cos(x)+np.cos(y)+np.cos(x+y))
        if c==0 :
           print(x,y)

Pomme de terre

Même pour des valeurs de $x$ proches de $0$, il se peut que $|\sin(x)+\sin(y) - \sin(x+y)|$ soit deux fois plus grand que $|x|$. En conséquence $|x| < \epsilon$ n'implique pas $|\sin(x) + \sin(y) - \sin(x+y)| < \epsilon$, ce qui peut créer un biais.

Imaginons par exemple que $Y = D\pi - X$ avec $X$ de loi uniforme sur $[0,\pi]$ et $D$ de loi uniforme sur $\{0,1\}$, indépendante de $X$.

Lorsque $\epsilon \to 0$, l'espérance de $\cos(X)+\cos(Y) + \cos(X+Y)$ tendra alors vers $\frac53$ en conditionnant par $(|\sin(X)+\sin(Y)-\sin(X+Y)| < \epsilon) \cap (|X| < \epsilon)$, tandis qu'elle tendra vers $1$ en conditionnant  par $(|X| < \epsilon)$.

Pour éviter ce problème, il suffirait de considérer $|X| < \frac\epsilon2$.

Guego

Pomme de terre : Intéressant, c'est sûrement quelque chose comme ça.

J'ai fait des simulations, en ne retenant que les couples $(X,Y)$ pour lesquels $|\sin(X)+\sin(Y)-\sin(X+Y)| < 0.01$ d'une part, et j'ai refait des simulations en ne retenant que les couples $(X,Y)$ pour lesquels $|X| < 0.01$ d'autre part. Graphiquement, si on représente les couples $(X,Y)$ obtenus, ça ressemble à ça : 

(en bleu $|X|<0.01$, en rouge $|\sin(X)+\sin(Y)-\sin(X+Y)| < 0.01$)

Si on zoome sur la bande supérieure de $X=0$, ça donne :

On constate que la répartition des ordonnées n'est pas du tout la même : pour le bleu, le gros des points est entre 1.5 et 3, alors que pour le rouge, plus on monte, moins il y a de points. La loi de $Y$ sachant $|X|$ petit n'est pas du tout la même que la loi de $Y$ sachant $|\sin(X)+\sin(Y)-\sin(X+Y)|$ et $|X|$ petits.