Générer des nombres aléatoires suite à une distribution dans un intervalle

J'ai besoin de générer des nombres aléatoires suivant la distribution normale dans l'intervalle $(a,b)$ . (Je travaille dans R.)

Je sais que la fonction rnorm(n,mean,sd)générera des nombres aléatoires suivant la distribution normale, mais comment définir les limites d'intervalle à l'intérieur de cela? Existe-t-il des fonctions R particulières disponibles pour cela?

Il semble que vous souhaitiez simuler à partir d'une distribution tronquée , et dans votre exemple spécifique, une normale tronquée .

Il existe différentes méthodes pour le faire, certaines simples, d'autres relativement efficaces.

Je vais illustrer quelques approches sur votre exemple normal.

1. Voici une méthode très simple pour en générer un à la fois (dans une sorte de pseudocode):

   $\tt{repeat}$ générer$x_i$ partir de N (moyenne, sd)$\tt{until}$ inférieur$\leq x_i\leq$ supérieur

   

   Si la majeure partie de la distribution est dans les limites, c'est 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 une par une en calculant la zone dans les limites et générer suffisamment de valeurs pour être presque certain qu'après avoir jeté les valeurs en dehors des limites, vous avez toujours autant de valeurs que nécessaire.

2. Vous pouvez utiliser accepter-rejeter avec une fonction de majorisation appropriée sur l'intervalle (dans certains cas, l'uniforme sera suffisant). Si les limites étaient raisonnablement étroites par rapport au sd mais que vous n'étiez pas loin dans la queue, une majoration uniforme fonctionnerait bien avec la normale, par exemple.

   

3. Si vous avez un cdf raisonnablement efficace et un cdf inverse (comme pnormet qnormpour la 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 Wikipedia sur la normale tronquée . [En fait, cela revient à prendre un uniforme tronqué (tronqué aux quantiles requis, ce qui ne nécessite en fait aucun rejet, car ce n'est qu'un autre uniforme) et d'appliquer le cdf normal inverse à cela. Notez que cela peut échouer si vous êtes loin dans la queue]

   

4. Il existe d'autres approches; la même page Wikipedia mentionne l'adaptation de la méthode ziggurat , qui devrait fonctionner pour une variété de distributions.

Le même lien Wikipédia mentionne deux packages spécifiques (tous deux sur CRAN) avec des fonctions pour générer des normales tronquées:

Le MSMpackage dans R a une fonction rtnorm,, qui calcule les tirages à partir d'une normale tronquée. Le truncnormpackage dans R a également des fonctions pour puiser dans une normale tronquée.

En regardant autour, une grande partie de cela est couverte dans les réponses à d'autres questions (mais pas exactement en double car cette question est plus générale que la normale tronquée) ... voir la discussion supplémentaire dans

une. Cette réponse

b. La réponse de Xi'an ici , qui a un lien avec son article arXiv (avec d'autres réponses intéressantes).

L'approche rapide et sale consiste à utiliser la règle 68-95-99.7 .

Dans une distribution normale, 99,7% des valeurs se situent dans les 3 écarts-types de la moyenne. Donc, si vous définissez votre moyenne au milieu de votre valeur minimale et de votre valeur maximale souhaitées, et définissez votre écart-type à 1/3 de votre moyenne, vous obtenez (principalement) des valeurs qui se situent dans l'intervalle 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)

J'ai récemment fait face à ce même problème, essayant de générer des notes d'élèves aléatoires pour les données de test. Dans le code ci-dessus, j'ai utilisé pmaxet pminpour remplacer les valeurs hors limites par la valeur minimale ou maximale dans les limites. Cela fonctionne pour mon objectif, car je génère des quantités de données assez petites, mais pour des quantités plus importantes, cela vous donnera des bosses notables aux valeurs min et max. Ainsi, selon vos objectifs, il peut être préférable de supprimer ces valeurs, de les remplacer par des NAs ou de les "relancer" jusqu'à ce qu'elles soient dans les limites.

$a<b$

$\Phi$$X_1,...,X_N$$\mu$$\sigma^2$$^\dagger$$a<b$

Il n'y a pas de fonction intégrée pour les valeurs générées à partir de la distribution tronquée, mais il est trivial de programmer cette méthode en utilisant les fonctions ordinaires pour générer des variables aléatoires. Voici une Rfonction simple rtruncnormqui 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 s'agit d'une fonction vectorisée qui générera Ndes variables aléatoires IID à partir de la distribution normale tronquée. Il serait facile de programmer des fonctions pour d'autres 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$$\mu$$\sigma^2$