problème avec une simulation [Résolu]

[Stéphane Laurent]

Bonjour,

Je suis face à ce problème : je simule les réalisations d'une distribution et en traçant la densité estimée avec la fonction density(), je constate que celle-ci diffère de la distribution théorique.

Voici le code :

Code : Tout sélectionner

a<-10 ; b<-10
aa<-3; bb<-15

nb.sim<-30000
sim.x <-rgamma(nb.sim, aa, bb); sim.y<-rgamma(nb.sim, a, b)
sim <- a/aa*bb/b*sim.x/sim.y

plot(density(sim, from=0, to=2))
xx<-seq(0, 2, by=0.01)
lines(xx, df(xx, 2*aa, 2*a), col=2)

Le simulations contenues dans le vecteur sim sont théoriquement distribuées selon une loi de Fisher de paramètre 2*aa et 2*a. Les deux courbes devraient donc se superposer.

Cela est de moins en moins le cas quand la valeur de a diminue (ce qui correspond à une augmentation de la variance de la Fisher).

Auriez-vous une idée d'où peut venir le problème ?

[Stéphane Laurent]

Le problème est bien plus clair ainsi, il n'y a pas besoin de passer par les Gamma pour simuler :

Code : Tout sélectionner

df1<-12 ; df2<-4

xx<-seq(0, 2, by=0.001)
plot(xx, df(xx, df1, df2), type="l")
nb.sim <- 40000
sim <- rf(nb.sim, df1, df2)
lines(density(sim), col=2)

PS : j'ai aussi essayé de simuler la Fisher en transformant des simulations de Beta, et il y a aussi des problèmes pour des petits df2

[Logez Maxime]

Bonjour,

Je ne pense pas que le problème viennent de df, mais de density :

Code : Tout sélectionner

hist(sim,prob=T,xlim=c(0,2),breaks=c(seq(0,2,le=30),max(sim)))
curve(df(x,df1,df2),add=TRUE,col=2)


La on voit bien que les simulations suivent une loi de Fisher. Je pense qu'il te faut modifier un peu les paramètres de density pour coller un peu plus aux simulations :

Code : Tout sélectionner

lines(density(sim,n=20000),col="green")

Maxime

[Stéphane Laurent]

Merci, je me suis rendu compte que le problème ne vient pas de df() et j'ai supprimé cette question avant de lire ta réponse.

Mais comment savoir quelle est la valeur appropriée de n dans la fonction density() ?

[Logez Maxime]

Re,

Je ne connais pas de moyen direct de calculer un n, mais déjà je pense que ça dépend du nombre de valeurs que tu as dans ta distribution ici tu as N=40000 et n=512 par défaut ce qui ne semble pas terrible et après je pense que le meilleur moyen est graphique. Si tu représentes n=512 tu te rends bien compte que ça ne va pas du tout, alors que si tu augmentes n alors tu vois que tu te rapproches de plus en plus de la distribution (par rapport à l'histogramme).

Maxime

[Stéphane Laurent]

Oui bonne idée de comparer avec l'histogramme. Merci.