Simuler une normale contrainte sur la limite inférieure ou supérieure dans R

Je voudrais générer des données aléatoires à partir d'une distribution normale contrainte en utilisant R.

Par exemple, je pourrais vouloir simuler une variable à partir d'une distribution normale avec mean=3, sd= 2et toutes les valeurs supérieures à 5 sont rééchantillonnées à partir de la même distribution normale.

Ainsi, pour la fonction générale, je pourrais faire ce qui suit.

rnorm(n=100, mean=3, sd=2)

J'ai ensuite eu quelques réflexions:

• Itérer une ifelsefonction avec une boucle qui se répète jusqu'à ce que toutes les valeurs soient contraintes de se trouver dans les limites.
• Simulez beaucoup plus de valeurs que nécessaire et prenez la première nqui satisfait la contrainte.
• Évitez les simulateurs de variables normales vectorisés et utilisez plutôt une boucle for avec un do while inside pour simuler chaque observation une à la fois et boucle si nécessaire.

Tout ce qui précède semble un peu maladroit.

Question

• Quel est un moyen simple de simuler une variable normale aléatoire contrainte dans R à partir de la normale avec une moyenne = 3, sd = 2 et max = 5?
• Plus généralement, quelle est une bonne manière générale d'incorporer des contraintes dans des variables simulées dans R?

C'est ce qu'on appelle une distribution normale tronquée:

http://en.wikipedia.org/wiki/Truncated_normal_distribution

Christian Robert a écrit sur une approche pour le faire pour une variété de situations (en utilisant différentes selon l'endroit où se trouvaient les points de troncature) ici:

Robert, CP (1995) "Simulation de variables normales tronquées",
Statistics and Computing, Volume 5, Numéro 2, juin, pp 121-125

Document disponible sur http://arxiv.org/abs/0907.4010

Ceci discute un certain nombre d'idées différentes pour différents points de troncature. Ce n'est pas la seule façon de les aborder par tous les moyens, mais il a généralement de très bonnes performances. Si vous voulez faire beaucoup de normales tronquées différentes avec différents points de troncature, ce serait une approche raisonnable. Comme vous l'avez noté, msm::tnormest basé sur l'approche de Robert, alors qu'il truncnorm::truncnormmet en œuvre l'échantillonneur d'acceptation-rejet de Geweke (1991); cela est lié à l'approche de l'article de Robert. Notez que msm::tnormcomprend les fonctions densité, cdf et quantile (cdf inverse) de la manière habituelle R.

Une référence plus ancienne avec une approche est le livre de Luc Devroye ; depuis qu'il est épuisé, il a récupéré le copyright et l'a rendu disponible en téléchargement.

Votre exemple particulier est le même que l'échantillonnage d'une normale standard tronquée à 1 (si $t$ est le point de troncature, $(t-\mu)/\sigma = (5-3)/2 = 1$), puis mettre le résultat à l'échelle (multiplier par $\sigma$ et ajouter $\mu$).

Dans ce cas précis, Robert suggère que votre idée (dans la deuxième ou la troisième incarnation) est tout à fait raisonnable. Vous obtenez une valeur acceptable environ 84% du temps et générez donc environ$1.19 n$ les normales en moyenne (vous pouvez calculer les limites de manière à générer suffisamment de valeurs à l'aide d'un algorithme vectorisé, disons 99,5% du temps, puis générer de temps en temps les dernières moins efficacement - même une à la fois).

Il y a aussi une discussion sur une implémentation dans le code R ici (et dans Rccp dans une autre réponse à la même question, mais le code R est en fait plus rapide). Le code R simple génère 50000 normales tronquées en 6 millisecondes, bien que cette normale tronquée particulière ne coupe que les queues extrêmes, donc une troncature plus substantielle signifierait que les résultats seraient plus lents. Il met en œuvre l'idée de générer «trop» en calculant combien il devrait générer pour être presque certain d'en avoir assez.

Si j'avais besoin d'un seul type particulier de normale tronquée plusieurs fois, je chercherais probablement à adapter une version de la méthode ziggurat, ou quelque chose de similaire, au problème.

En fait, il semble que Nicolas Chopin ait déjà fait cela, donc je ne suis pas la seule personne à avoir pensé à:

http://arxiv.org/abs/1201.6140

Il discute de plusieurs autres algorithmes et compare le temps pour 3 versions de son algorithme avec d'autres algorithmes pour générer 10 ^ 8 normales aléatoires pour divers points de troncature.

Sans surprise, son algorithme s'avère relativement rapide.

D'après le graphique de l'article, même le plus lent des algorithmes avec lesquels il se compare aux (pour eux) les pires points de troncature génèrent $10^8$ valeurs en environ 3 secondes - ce qui suggère que n'importe lequel des algorithmes discutés peut être acceptable s'il est raisonnablement bien implémenté.

Edit: celui dont je ne suis pas certain est mentionné ici (mais peut-être qu'il se trouve dans l'un des liens) est de transformer (via un cdf normal inverse) un uniforme tronqué - mais l'uniforme peut être tronqué en générant simplement un uniforme dans les limites de la troncature . Si le cdf normal inverse est rapide, c'est à la fois rapide et facile et fonctionne bien pour une large gamme de points de troncature.

Dans la continuité des références de @ glen_b et en se concentrant exclusivement sur l'implémentation de R. Il existe quelques fonctions conçues pour échantillonner à partir d'une distribution normale tronquée:

• rtruncnorm(100, a=-Inf, b=5, mean=3, sd=2)dans le truncnormpaquet
• rtnorm(100, 3, 2, upper=5)dans le msmpaquet