Jai besoin de générer des nombres aléatoires suivant une distribution normale dans lintervalle $ (a, b) $. (Je travaille en R.)
Je sais que la fonction rnorm(n,mean,sd)
générera des nombres aléatoires suivant une distribution normale, mais comment définir les limites dintervalle à lintérieur de cela? Y a-t-il des fonctions R particulières disponibles pour cela?
Commentaires
Réponse
Il semble que vous souhaitiez simuler à partir dune distribution tronquée , et dans votre exemple spécifique , une normale tronquée .
Il existe une variété de méthodes pour le faire, certaines simples, dautres relativement efficace.
Je vais illustrer quelques approches sur votre exemple normal.
-
Voici une méthode très simple pour en générer une à la fois (dans une sorte de pseudocode ):
$ \ tt {repeat} $ génère $ x_i $ à partir de N (moyenne, sd) $ \ tt {jusquà} $ inférieur $ \ leq x_i \ leq $ supérieur
Si la plupart de la distribution est dans les limites, cest assez raisonnable mais cela peut devenir assez lent si vous générez presque toujours en dehors des limites.
Dans R, vous pouvez éviter la boucle un à la fois en calculant la zone dans les limites et en générant suffisamment de valeurs pour être presque certain quaprès avoir jeté les valeurs en dehors des limites, vous aviez toujours autant de valeurs que nécessaire.
-
Vous pouvez utiliser accepter-rejeter avec une fonction de majoration appropriée sur lintervalle (dans certains cas, être assez bon). Si les limites étaient raisonnablement étroites par rapport à la s.d. mais vous nétiez pas loin dans la queue, une majoration uniforme fonctionnerait bien avec la normale, par exemple.
-
Si vous avez un cdf raisonnablement efficace et un cdf inverse (comme
pnorm
etqnorm
pour le distribution normale dans R), vous pouvez utiliser la méthode inverse-cdf décrite dans le premier paragraphe de la section de simulation de la page Wikipédia sur la normale tronquée . [En vigueur cest la même chose que de prendre un uniforme tronqué (tronqué aux quantiles requis, ce qui ne nécessite en fait aucun rejet, puisque cest juste un autre uniforme) et dappliquer le cdf normal inverse à cela. Notez que cela peut échouer si vous « êtes loin dans la queue] -
Il existe dautres approches; la même page Wikipédia mentionne ladaptation de la méthode ziggurat , qui devrait fonctionner pour une variété de distributions.
Le même lien Wikipedia mentionne deux packages spécifiques (tous deux sur CRAN) avec des fonctions pour générer des normales tronquées:
Le package
MSM
dans R a une fonction,rtnorm
, qui calcule les tirages à partir dun tronqué normal. Le packagetruncnorm
dans R a également des fonctions pour dessiner à partir dune normale tronquée.
En regardant autour de vous, une grande partie de cela est couverte dans les réponses à dautres questions (mais pas exactement en double puisque cette question est plus générale que la normale tronquée) … voir la discussion supplémentaire dans
a. Cette réponse
b. La réponse de Xi « an » ici , qui contient un lien vers son article arXiv (avec dautres réponses intéressantes).
Réponse
Lapproche rapide et sale consiste à utiliser la règle 68-95-99.7 .
Dans une distribution normale, 99,7% des valeurs se situent à moins de 3 écarts types de la moyenne. Ainsi, si vous définissez votre moyenne au milieu de la valeur minimale et maximale souhaitée et que vous définissez votre écart-type sur 1/3 de votre moyenne, vous obtenez (principalement) des valeurs comprises dans lintervalle souhaité. Ensuite, vous pouvez simplement nettoyer le reste.
minVal <- 0 maxVal <- 100 mn <- (maxVal - minVal)/2 # Generate numbers (mostly) from min to max x <- rnorm(count, mean = mn, sd = mn/3) # Do something about the out-of-bounds generated values x <- pmax(minVal, x) x <- pmin(maxVal, x)
Jai récemment rencontré ce même problème, en essayant de générer notes aléatoires des élèves pour les données de test. Dans le code ci-dessus, jai « utilisé pmax
et pmin
pour remplacer les valeurs hors limites par les limites minimales ou maximales évaluer.Cela fonctionne dans mon but, car je génère des quantités de données assez petites, mais pour de plus grandes quantités, cela vous donnera des bosses notables aux valeurs min et max. Donc, en fonction de vos objectifs, il peut être préférable de supprimer ces valeurs, de les remplacer avec des NA
, ou « relancez-les » jusquà ce quils « soient dans les limites.
Commentaires
- Pourquoi se donner la peine de faire ça? Il est si simple de générer des nombres aléatoires normaux et de supprimer ceux qui nécessitent une troncature quil nest ‘ t nécessaire dêtre compliqué à ce sujet, sauf si la troncature souhaitée est proche de 100% de la surface de la densité.
- Jai peut-être ‘ mal interpréter la question initiale. Je suis tombé sur cette question en essayant de comprendre comment réaliser une tâche de programmation non directement liée aux statistiques dans R, et jai ‘ remarqué que cette page est une pile de statistiques , pas un échange de pile de programmation. 🙂 Dans mon cas, je voulais générer une quantité spécifique dentiers aléatoires, avec des valeurs allant de 0 à 100, et je voulais que les valeurs générées tombent sur une belle courbe en cloche dans cette plage. Depuis que jai écrit ceci, ‘ jai réalisé que
sample(x=min:max, prob=dnorm(...))
est peut-être un moyen plus simple de le faire. - @Glen_b Aaron Wells mentionne
sample(x=min:max, prob=dnorm(...))
qui semble un peu plus court que votre réponse. - Mais notez que lastuce
sample()
nest utile que si vous ‘ essayez de choisir des entiers aléatoires ou un autre ensemble de valeurs discrètes prédéfinies.
Réponse
Aucune des réponses ici ne donne une méthode efficace pour générer des variables normales tronquées qui nimpliquent pas le rejet de variables arbitrairement grandes nombre de valeurs générées. Si vous souhaitez générer des valeurs à partir dune distribution normale tronquée, avec des limites inférieure et supérieure spécifiées $ a < b $ , ceci peut être fait — sans rejet — en générant des quantiles uniformes sur la plage de quantiles autorisée par la troncature, et en utilisant échantillonnage par transformation inverse pour obtenir les valeurs normales correspondantes .
Soit $ \ Phi $ le CDF de la distribution normale standard. Nous voulons générer $ X_1, …, X_N $ à partir dune distribution normale tronquée (avec paramètre de moyenne $ \ mu $ et paramètre de variance $ \ sigma ^ 2 $ ) $ ^ \ dagger $ avec un et limites de troncature supérieures $ a < b $ . Cela peut être fait comme suit:
$$ X_i = \ mu + \ sigma \ cdot \ Phi ^ {- 1} (U_i) \ quad \ quad \ quad U_1, …, U_N \ sim \ text {IID U} \ Big [\ Phi \ Big (\ frac {a- \ mu} {\ sigma} \ Big), \ Phi \ Big (\ frac {b- \ mu} {\ sigma} \ Big) \ Big]. $$
Il ny a pas de fonction intégrée pour les valeurs générées à partir de la distribution tronquée, mais il est simple de programmer cette méthode en utilisant le fonctions ordinaires pour générer des variables aléatoires. Voici une simple R
fonction rtruncnorm
qui implémente cette méthode en quelques lignes de code.
rtruncnorm <- function(N, mean = 0, sd = 1, a = -Inf, b = Inf) { if (a > b) stop("Error: Truncation range is empty"); U <- runif(N, pnorm(a, mean, sd), pnorm(b, mean, sd)); qnorm(U, mean, sd); }
Il sagit dune fonction vectorisée qui générera des variables aléatoires N
IID à partir de la distribution normale tronquée. Il serait facile de programmer des fonctions pour dautres distributions tronquées via la même méthode. Il ne serait pas non plus trop difficile de programmer les fonctions de densité et de quantile associées pour la distribution tronquée.
$ ^ \ dagger $ Notez que la troncature modifie la moyenne et la variance de la distribution, donc $ \ mu $ et $ \ sigma ^ 2 $ sont et non la moyenne et la variance de la distribution tronquée.
Réponse
Trois méthodes ont fonctionné pour moi:
-
en utilisant sample () avec rnorm ():
sample(x=min:max, replace= TRUE, rnorm(n, mean))
-
en utilisant le package msm et la fonction rtnorm:
rtnorm(n, mean, lower=min, upper=max)
-
en utilisant rnorm () et en spécifiant les limites inférieure et supérieure, comme Hugh la indiqué ci-dessus:
sample <- rnorm(n, mean=mean); sample <- sample[x > min & x < max]
x <- rnorm(n, mean, sd); x <- x[x > lower.limit & x < upper.limit]