Simulation de Monte Carlo
Bonjour tout le monde,
Je viens écrire sur ce forum après des jours d'arrachage de cheveux, je suis un peu désespérée je dois dire ! Voilà mon problème : mon code me sert à faire une simulation de monte carlo. Mes données sont N points (par exemple 10000), et une certaine probabilité (l est un indice entier, qui fait varier p). Chaque "tirage" est indexé par k (k va de 1 à m), et pour chacun de ces tirages, pour chaque point, je tire un nombre aléatoire (Xi dans mon programme). Si Xi>p/m, alors le point continue ça route. Si Xi<=p/m le point est supprimé. Au tirage d'après, je recommence la même chose, avec N-(nombre de points supprimés) points... Ainsi de suite. Je tiens à préciser qu'il est nécessaire pour moi de faire des étapes indexées par k, que je ne peux pas juste faire un tirage pour chaque point en prenant en disant que la probabilité de survie d'un point est (1-p/m)^m (ce serait trop facile :().
Voici donc mon code :
Code:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
|
S=N !S et T servent à compter les points.
T=N
do k=1,m
do i=1,S
call RANDOM_SEED
call RANDOM_NUMBER(xi)
if (xi<=p/m) then
T=T-1 !(Position 1)
else
T=T
end if
end do
S=T !(Position 2)
end do
print*,S |
Je lance cette simulation pour plusieurs valeurs de p, et les résultats que j'obtiens ne sont pas du tout ceux que j'attends. Du coup, j'ai installé des tests pour comprendre ce qu'il se passe, du type
en particulier aux positions 1 et 2 que j'ai inscrit en commentaire dans mon code. Et là, il se passe des choses bien étranges ...
Quand j'observe en position 1, tout semble se dérouler comme j'en ai envie, le nombre T diminue de façon "régulière" (pour k->k+1 j'entends).
Quand j'observe en position 2, globalement k augmente, T reste constant, jusqu'à chuter brusquement pour k->k+1 donné, puis redevient constant pour un moment ... ect. Ce n'est pas du tout le comportement que j'attends!
Voilà pourquoi je demande votre aide, tout ça me laisse vraiment perplexe.
Merci d'avance pour votre temps et bonne journée,
Amélie
