Simular la normal restringida en el límite inferior o superior en R

Question

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

Por ejemplo, podría querer simular una variable de una distribución normal con mean=3, sd= 2 y cualquier valor superior a 5 se vuelve a muestrear a partir de la misma distribución normal.

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

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

Entonces se me ocurrieron algunas ideas:

• Iterar un ifelse con un bucle que se repite hasta que todos los valores se encuentran dentro de los límites.
• Simula muchos más valores de los necesarios y toma el primero n que satisfagan la restricción.
• Evite los simuladores de variables normales vectorizadas y utilice en su lugar un bucle for con un do while dentro para simular cada observación de una en una y hacer un bucle cuando sea necesario.

Todo lo anterior parece un poco torpe.

Pregunta

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

Answer 1

Se denomina distribución normal truncada:

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

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

Robert, C.P. (1995) "Simulación de variables normales truncadas",
Statistics and Computing, Volume 5, Issue 2, June, pp 121-125

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

Aquí se discuten varias ideas diferentes para distintos puntos de truncamiento. No es la única manera de abordar estos de ninguna manera, pero por lo general tiene un rendimiento bastante bueno. Si quieres hacer un montón de diferentes normales truncadas con varios puntos de truncamiento, sería un enfoque razonable. Como usted ha señalado, msm::tnorm se basa en el planteamiento de Robert, mientras que truncnorm::truncnorm implementa el muestreador de aceptación-rechazo de Geweke (1991); está relacionado con el enfoque del artículo de Robert. Obsérvese que msm::tnorm incluye funciones de densidad, CDF y cuantiles (CDF inversa) en la forma habitual R moda.

Una referencia más antigua con un enfoque es El libro de Luc Devroye Desde que se agotó, ha recuperado los derechos de autor y lo ha puesto a disposición de los lectores como descarga.

Su ejemplo particular es el mismo que el muestreo de una normal estándar truncada en 1 (si $t$ es el punto de truncamiento, $(t-\mu)/\sigma = (5-3)/2 = 1$ ) y, a continuación, escalar el resultado (multiplicar por $\sigma$ y añada $\mu$ ).

En ese caso concreto, Robert sugiere que su idea (en la segunda o tercera encarnación) es bastante razonable. Obtienes un valor aceptable aproximadamente el 84% de las veces y así generas unos $1.19 n$ normales de media (se pueden calcular límites de forma que se generen suficientes valores utilizando un algoritmo vectorizado, digamos, el 99,5% de las veces, y luego, de vez en cuando, generar los últimos de forma menos eficiente, incluso de uno en uno).

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

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

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

http://arxiv.org/abs/1201.6140

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

Como era de esperar, su algoritmo resulta ser relativamente rápido.

A partir del gráfico del artículo, incluso el más lento de los algoritmos con los que se compara en los (para ellos) peores puntos de truncamiento está generando $10^8$ en unos 3 segundos, lo que sugiere que cualquiera de los algoritmos comentados puede ser aceptable si se aplica razonablemente bien.

Edición: Uno que no estoy seguro de que se menciona aquí (pero tal vez es en uno de los enlaces) es transformar (a través de cdf normal inversa) un truncado uniforme - pero el uniforme puede ser truncado por la simple generación de un uniforme dentro de los límites de truncamiento. Si la normal inversa cdf es rápida esto es a la vez rápido y fácil y funciona bien para una amplia gama de puntos de truncamiento.

Answer 2

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

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

Answer 3

Un ejemplo de uso de la CDF inversa (función cuantil) como sugiere @Glen_b

Puede utilizar runif para generar cuantiles aleatorios y luego pasar estos cuantiles a, por ejemplo qnorm (o cualquier otra distribución) para encontrar los valores a los que corresponden estos cuantiles para la distribución dada.

Si sólo genera cuantiles dentro de un intervalo específico, trunca la distribución. Podemos utilizar la FCD (por ejemplo pnorm ) para hallar los límites de los cuantiles, que corresponden a un truncamiento dado.

rtruncnorm <- function(n, mu, sigma, low, high) {
  # find quantiles that correspond the the given low and high levels.
  p_low <- pnorm(low, mu, sigma)
  p_high <- pnorm(high, mu, sigma)

  # draw quantiles uniformly between the limits and pass these
  # to the relevant quantile function.
  qnorm(runif(n, p_low, p_high), mu, sigma)
}

samples <- rtruncnorm(1000, 3, 2, low = -Inf, high = 5)

max(samples)
#> [1] 4.996336

hist(samples)