¿Qué trunca distribución decir?

Question

En un artículo de investigación sobre análisis de sensibilidad de un modelo de ecuación diferencial ordinaria de un sistema dinámico, el autor, siempre que la distribución de un modelo de parámetros como la distribución Normal (media=1e-4, std=3e-5) truncado en el intervalo [0.5 e-4 1.5 e-4]. Él, a continuación, utiliza muestras de esta truncada de distribución para las simulaciones del modelo. ¿Qué significa tener un truncado de distribución y muestra de esta truncada de distribución?

Me podía venir de dos maneras de hacer esto:

• Muestra de una distribución Normal, pero no tienen en cuenta todos los valores aleatorios que caen fuera del rango especificado antes de simulaciones.
• De alguna manera obtener un especial de "Truncado Normal" de distribución y obtener muestras de ello.

Son estos válida y equivalente enfoques?

Creo que en el primer caso, si uno se para de la parcela experimental de cdf/pdf de la muestra, no se vería como una distribución Normal, debido a que las curvas no se extienden a $\pm\infty$.

Answer 1

La simulación de la normal $\mathcal{N}(\mu,\sigma^2)$ distribución hasta que el resultado cae dentro de un intervalo de $(a,b)$ está bien cuando la probabilidad $$ \varrho = \int_a^b \varphi_{\mu\sigma^2}(x)\,\text{d} x $$ es lo suficientemente grande. Si es demasiado pequeño, este procedimiento es muy costoso ya que el número medio de sorteos para una aceptación es $1/\varrho$.

Como se describe en Monte Carlo Métodos Estadísticos (Capítulo 2, Ejemplo 2.2), así como en mi arXiv papel, de una manera más eficiente para simular este trunca normal es utilizar un aceptar-rechazar método basado en un exponencial $\mathcal{E}(\alpha)$ distribución.

Considerar, sin pérdida de la generalidad, en el caso de $\mu = 0$$\sigma = 1$. Al $b=+\infty$, un potencial instrumental de distribución es la traducción de la distribución exponencial, $\mathcal{E} (\alpha,{ a})$, con la densidad $$ g_{\alpha}(z) = \alpha e^{- \alpha(z - {a})} \; \mathbb{I}_{z \geq {a }} \;. $$ La relación de $$ p_{a,\infty}(z)/g_{\alpha}(z) \propto e^{- \alpha(z - a )}e^{-z^{2}/2} $$ es entonces limitada por $\exp(\alpha^{2}/2 - \alpha{a })$ si $\alpha > a$ y $\exp(- a^{2}/2)$ de lo contrario. La correspondiente (superior) es obligado $$ \begin{cases} 1/\alpha \; \exp (\alpha^{2}/2 - \alpha{a }) & \hbox{if } \alpha > a , \cr 1/\alpha \; \exp (- a^{2}/2) & \hbox{otherwise.} \cr \end{casos} $$ La primera expresión es minimizado por \begin{equation} \alpha^{*} = \frac{1}{2}a + \frac{1}{2} \sqrt{a^2 + 4}\;,\qquad (1) \end{equation} mientras que $\tilde\alpha = a $ minimiza el segundo enlazado. La elección óptima de $\alpha$ es por lo tanto (1).

Answer 2

Para truncar una distribución es restringir sus valores en un intervalo y re-normalización de la densidad, de modo que la integral sobre ese rango es 1.

Así, para truncar la $N(\mu, \sigma^{2})$ distribución para un intervalo de $(a,b)$ sería para generar una variable aleatoria que tiene una densidad de

$$ p_{a,b}(x) = \frac{ \phi_{\mu, \sigma^{2}}(x) }{ \int_{a}^{b} \phi_{\mu, \sigma^{2}}(y) dy } \cdot \mathcal{I} \{ x \in (a,b) \} $$

donde $\phi_{\mu, \sigma^{2}}(x)$ $N(\mu, \sigma^2)$ densidad. Puedes probar a partir de esta densidad en un número de maneras. Una forma (la forma más sencilla que puedo pensar) para hacer esto sería para generar $N(\mu, \sigma^2)$ valores y tirar las que caen fuera de la $(a,b)$ intervalo, como usted ha mencionado. Así que, sí, los dos balas que se enumeran podría cumplir el mismo objetivo. También, tienes razón en que la evidencia empírica de la densidad (o histograma) de las variables a partir de esta distribución no se haría extensivo a $\pm \infty$. Sería restringido a $(a,b)$, por supuesto.

Answer 3

Muestra de una distribución Normal, pero no tienen en cuenta todos los valores aleatorios caen fuera del rango especificado antes de simulaciones.

Este método es correcto, pero, como es mencionado por @Xi'an en su respuesta, tomaría mucho tiempo cuando el rango es pequeño (más precisamente, cuando su medida es pequeño debajo de la distribución normal).

Como cualquier otro tipo de distribución, se podría utilizar el método de inversión de $F^{-1}(U)$ (también se llama inversa de la transformación de muestreo), donde $F$ es la función acumulativa de la distribución de los intereses y $U\sim\text{Unif}(0,1)$. Al $F$ es la distribución obtenida mediante el truncamiento de algunos de distribución de $G$ en algunos intervalo de $(a,b)$, esto es equivalente a la de la muestra $G^{-1}(U)$$U\sim\text{Unif}\bigl(G(a),G(b)\bigr)$.

Sin embargo, y esto ya es mencionado por @Xi an en un comentario, para algunas situaciones de la inversión método requiere una precisa evaluación de la función cuantil $G^{-1}$, y yo añadiría también se requiere un rápido cálculo de $G^{-1}$. Al $G$ es una distribución normal, la evaluación de las $G^{-1}$ es bastante lento, y no es muy preciso, para que los valores de $a$ $b$ fuera de la "gama" de $G$.

Simular un truncado de distribución utilizando la importancia de muestreo

Una posibilidad es utilizar importancia de muestreo. Considere el caso en el estándar de la distribución Gaussiana ${\cal N}(0,1)$. Olvidar la anterior anotación, ahora vamos a $G$ ser la distribución de Cauchy. Las dos mencionadas cumplimiento de los requisitos para $G$ : simplemente se ha $\boxed{G(q)=\frac{\arctan(q)}{\pi}+\frac12}$$\boxed{G^{-1}(q)=\tan\bigl(\pi(q-\frac12)\bigr)}$. Por lo tanto, la truncada distribución de Cauchy es fácil de probar por la inversión de método y es una buena elección de la variable instrumental de la importancia de muestreo de la distribución normal truncada.

Después de un poco de simplificaciones, muestreo $U\sim\text{Unif}\bigl(G(a),G(b)\bigr)$ y teniendo en $G^{-1}(U)$ es equivalente a tomar $\tan(U')$$U'\sim\text{Unif}\bigl(\arctan(a),\arctan(b)\bigr)$:

a <- 1
b <- 5
nsims <- 10^5
sims <- tan(runif(nsims, atan(a), atan(b)))

Ahora uno tiene que calcular el peso de cada valor muestreado $x_i$, definido como el cociente $\phi(x)/g(x)$ de los dos densidades de hasta normalización, por lo tanto podemos tomar $$ w(x) = \exp(-x^2/2)(1+x^2), $$ pero podría ser más seguro para tener el registro de los pesos:

log_w <- -sims^2/2 + log1p(sims^2)
w <- exp(log_w) # unnormalized weights
w <- w/sum(w)

La muestra ponderada $(x_i,w(x_i))$ permite estimar la medida de cada intervalo de $[u,v]$ por debajo del objetivo de la distribución, sumando los pesos de cada valor muestreado de caer dentro del intervalo:

u <- 2; v<- 4
sum(w[sims>u & sims<v])
## [1] 0.1418

Esto proporciona una estimación de la meta acumulada de la función. Podemos obtener rápidamente y la trama con la spatsat paquete de:

F <- spatstat::ewcdf(sims,w)
# estimated F:
curve(F(x), from=a-0.1, to=b+0.1)
# true F:
curve((pnorm(x)-pnorm(a))/(pnorm(b)-pnorm(a)), add=TRUE, col="red")

# approximate probability of u<x<v:
F(v)-F(u)
## [1] 0.1418

Por supuesto, la muestra $(x_i)$ definitivamente no es una muestra de la distribución de destino, sino de la instrumental de Cauchy de distribución, y se obtiene una muestra de la distribución de destino mediante la realización ponderado de remuestreo, por ejemplo utilizando el muestreo multinomial:

msample <- rmultinom(1, nsims, w)[,1]
resims <- rep(sims, times=msample)
hist(resims) 

mean(resims>u & resims<v)
## [1] 0.1446

Otro método: rápido inversa de la transformación de muestreo

Olver y Townsend desarrollado un método de muestreo para una amplia clase de una distribución continua. Se implementa en el chebfun2 biblioteca para Matlab así como la ApproxFun biblioteca para Julia. Recientemente he descubierto esta biblioteca y suena muy prometedor (no sólo para el muestreo aleatorio). Básicamente esta es la inversión de método, pero el uso de potentes aproximaciones de la cdf y la inversa de la cdf. La entrada es el destino de la función de densidad hasta la normalización.

La muestra es simplemente generados por el siguiente código:

using ApproxFun
f = Fun(x -> exp(-x.^2./2), [1,5]);
nsims = 10^5;
x = sample(f,nsims);

Como se comprueba a continuación, se obtiene un estimado de la medida del intervalo de $[2,4]$ cerca el uno obtenidos anteriormente por la importancia de muestreo:

sum((x.>2) & (x.<4))/nsims
## 0.14191