Simular normal restringido en límite inferior o superior en R

Me gustaría generar datos aleatorios a partir de una distribución normal restringida usando R.

Por ejemplo, es posible que desee simular una variable de una distribución normal con mean=3, sd= 2y cualquier valor mayor que 5 se muestrea de la misma distribución normal.

Por lo tanto, para la función general, podría hacer lo siguiente.

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

Entonces tuve algunos pensamientos:

• Itere una ifelsefunción con un bucle que se repite hasta que todos los valores estén restringidos a estar dentro de los límites.
• Simule muchos más valores de los requeridos y tome los primeros nque satisfagan la restricción.
• Evite los simuladores de variables normales vectorizadas y, en su lugar, use un bucle for con un do mientras está adentro para simular cada observación una a la vez y bucle donde sea necesario.

Todo lo anterior parece un poco torpe.

Pregunta

• ¿Cuál es una manera simple de simular una variable normal aleatoria restringida en R desde normal con media = 3, sd = 2 y máx = 5?
• En términos más generales, ¿cuál es una buena forma general de incorporar restricciones en variables simuladas en R?

Esto se llama distribución normal truncada:

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

Christian Robert escribió sobre un enfoque para hacerlo en una variedad de situaciones (usando diferentes dependiendo de dónde estaban los puntos de truncamiento) aquí:

Robert, CP (1995) "Simulación de variables normales truncadas",
Estadística e informática, Volumen 5, Número 2, junio, pp 121-125

Documento disponible en http://arxiv.org/abs/0907.4010

Esto discute una serie de ideas diferentes para diferentes puntos de truncamiento. No es la única forma de abordarlos de ninguna manera, pero generalmente tiene un rendimiento bastante bueno. Si desea hacer muchas normales truncadas diferentes con varios puntos de truncamiento, sería un enfoque razonable. Como notó, msm::tnormse basa en el enfoque de Robert, mientras truncnorm::truncnormimplementa el muestreador de aceptación / rechazo de Geweke (1991); Esto está relacionado con el enfoque en el artículo de Robert. Tenga en cuenta que msm::tnormincluye las funciones de densidad, cdf y cuantil (cdf inverso) de la Rmanera habitual .

Una referencia más antigua con un enfoque es el libro de Luc Devroye ; desde que se agotó, recuperó los derechos de autor y lo puso a disposición como descarga.

Su ejemplo particular es lo mismo que muestrear un estándar normal truncado en 1 (si es el punto de truncamiento, ), y luego escalar el resultado (multiplicar por y agregar ).$t$$(t-\mu)/\sigma = (5-3)/2 = 1$$\sigma$$\mu$

En ese caso específico, Robert sugiere que su idea (en la segunda o tercera encarnación) es bastante razonable. Obtiene un valor aceptable aproximadamente el 84% del tiempo y, por lo tanto, genera aproximadamente normales en promedio (puede calcular límites para generar suficientes valores usando un algoritmo vectorizado, digamos el 99.5% del tiempo, y luego de vez en cuando generar los últimos menos eficientemente, incluso uno a la vez).$1.19 n$

También se discute una implementación en el código R aquí (y en Rccp en otra respuesta a la misma pregunta, pero el código R allí es en realidad más rápido). El código R simple allí genera 50000 normales truncadas en 6 milisegundos, aunque esa normal truncada en particular solo corta las colas extremas, por lo que un truncamiento más sustancial significaría que los resultados fueron más lentos. Implementa la idea de generar "demasiados" calculando cuántos debería generar para estar seguro de obtener suficiente.

Si necesitara un tipo particular de normal truncado muchas veces, probablemente buscaría adaptar una versión del método zigurat, o algo similar, al problema.

De hecho, parece que Nicolas Chopin ya lo hizo, así que no soy la única persona que se le ocurrió:

http://arxiv.org/abs/1201.6140

Analiza varios otros algoritmos y compara el tiempo para 3 versiones de su algoritmo con otros algoritmos para generar 10 ^ 8 normales al azar para varios puntos de truncamiento.

Quizás, como era de esperar, su algoritmo resulta ser relativamente rápido.

Del gráfico en el documento, incluso el algoritmo más lento con el que se compara en (para ellos) los peores puntos de truncamiento están generando valores en aproximadamente 3 segundos, lo que sugiere que cualquiera de los algoritmos discutidos allí puede ser aceptable si es razonable Bien implementado.$10^8$

Editar: Uno que no estoy seguro se menciona aquí (pero tal vez está en uno de los enlaces) es transformar (a través de cdf normal inverso) un uniforme truncado, pero el uniforme se puede truncar simplemente generando un uniforme dentro de los límites de truncamiento . Si el cdf normal inverso es rápido, esto es rápido y fácil y funciona bien para una amplia gama de puntos de truncamiento.

Siguiendo las referencias de @ glen_b y centrándose exclusivamente en la implementación de R. Hay un par de funciones diseñadas para muestrear desde una distribución normal truncada:

• rtruncnorm(100, a=-Inf, b=5, mean=3, sd=2)en el truncnormpaquete
• rtnorm(100, 3, 2, upper=5)en el msmpaquete