¿Qué significa la distribución truncada?

En un artículo de investigación sobre el análisis de sensibilidad de un modelo de ecuación diferencial ordinaria de un sistema dinámico, el autor proporcionó la distribución de un parámetro del modelo como Distribución normal (media = 1e-4, std = 3e-5) truncada al rango [0.5e -4 1,5e-4]. Luego usa muestras de esta distribución truncada para simulaciones del modelo. ¿Qué significa tener una distribución truncada y una muestra de esta distribución truncada?

Podría idear dos formas de hacer esto:

• Muestra de una distribución Normal pero ignora todos los valores aleatorios que caen fuera del rango especificado antes de las simulaciones.
• De alguna manera, obtenga una distribución especial "Normal truncada" y obtenga muestras de ella.

¿Son estos enfoques válidos y equivalentes?

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

Truncar una distribución es restringir sus valores a un intervalo y volver a normalizar la densidad para que la integral sobre ese rango sea 1.

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

$$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 es la densidad N ( μ , σ 2 ) . Puede tomar muestras de esta densidad de varias maneras. Una forma (la forma más simple que se me ocurre) de hacer esto sería generar N ( μ , σ 2 ) valores y arrojar los que quedan fuera de ( a , b $\phi_{\mu, \sigma^{2}}(x)$$N(\mu, \sigma^2)$$N(\mu, \sigma^2)$$(a,b)$intervalo, como mencionaste. Entonces, sí, esas dos viñetas que enumeró cumplirían el mismo objetivo. Además, tiene razón en que la densidad empírica (o histograma) de las variables de esta distribución no se extendería a . Estaría restringido a ( a , b ) , por supuesto.$\pm \infty$$(a,b)$

Simular desde la distribución normal de hasta que el resultado cae dentro de un intervalo ( a , b ) está bien cuando la probabilidad ϱ = ∫ b a φ μ , σ 2 ( x $\mathcal{N}(\mu,\sigma^2)$$(a,b)$ es lo suficientemente grande. Si es demasiado pequeño, este procedimiento es demasiado costoso ya que el número promedio de sorteos para una aceptación es 1 / ϱ .

$$\varrho = \int_a^b \varphi_{\mu,\sigma^2}(x)\,\text{d} x$$ $1/\varrho$

Como se describe en los Métodos estadísticos de Monte Carlo (Capítulo 2, Ejemplo 2.2), así como en mi artículo arXiv , una forma más eficiente de simular esta normal truncada es utilizar un método de aceptación-rechazo basado en una distribución exponencial .$\mathcal{E}(\alpha)$

Considere, sin pérdida de generalidad, el caso y σ = 1 . Cuando b = + ∞ , una distribución instrumental potencial es la distribución exponencial traducida, E ( α , a ) , con densidad g α ( z ) = α e - α ( z - a$\mu = 0$$\sigma = 1$$b=+\infty$$\mathcal{E} (\alpha,{ a})$ La relación p un , ∞ ( z ) / g α ( z ) α e - α ( z - a ) e - z 2 / 2 está entonces limitada por exp ( α 2 / 2 - α un ) si α > una y por exp ( - un 2 / 2 ) de otro modo. El límite correspondiente (superior) es

$$g_{\alpha}(z) = \alpha e^{- \alpha(z - {a})} \; \mathbb{I}_{z \geq {a }} \;.$$ $$p_{a,\infty}(z)/g_{\alpha}(z) \propto e^{- \alpha(z - a )}e^{-z^{2}/2}$$ $\exp(\alpha^{2}/2 - \alpha{a })$$\alpha > a$$\exp(- a^{2}/2)$ La primera expresión se minimiza con α∗= $$\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{cases}$$ mientras que ˜ α = a minimiza el segundo límite. La elección óptima de α es por lo tanto (1). $$\begin{equation} \alpha^{*} = \frac{1}{2}a + \frac{1}{2} \sqrt{a^2 + 4}\;,\qquad (1) \end{equation}$$ $\tilde\alpha = a$$\alpha$

Muestra de una distribución Normal pero ignora todos los valores aleatorios que caen fuera del rango especificado antes de las simulaciones.

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

$F^{-1}(U)$$F$$U\sim\text{Unif}(0,1)$$F$$G$$(a,b)$$G^{-1}(U)$ con $U\sim\text{Unif}\bigl(G(a),G(b)\bigr)$.

Sin embargo, y esto ya lo menciona @ Xi'an en un comentario, para algunas situaciones el método de inversión requiere una evaluación muy precisa de la función cuantil$G^{-1}$, y agregaría que también requiere un cálculo rápido de$G^{-1}$. Cuando$G$ es una distribución normal, la evaluación de $G^{-1}$ es bastante lento y no es muy preciso para valores de $a$ y $b$ fuera del "rango" de $G$.

Simulate a truncated distribution using importance sampling

A possibility is to use importance sampling. Consider the case of the standard Gaussian distribution ${\cal N}(0,1)$. Forget the previous notations, now let $G$ be the Cauchy distribution. The two above mentionned requirements are fulfilled for $G$ : one simply has $\boxed{G(q)=\frac{\arctan(q)}{\pi}+\frac12}$ and $\boxed{G^{-1}(q)=\tan\bigl(\pi(q-\frac12)\bigr)}$. Therefore, the truncated Cauchy distribution is easy to sample by the inversion method and it is a good choice of the instrumental variable for importance sampling of the truncated normal distribution.

After a bit of simplifications, sampling $U\sim\text{Unif}\bigl(G(a),G(b)\bigr)$ and taking $G^{-1}(U)$ is equivalent to take $\tan(U')$ with $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)))

Now one has to calculate the weight for each sampled value $x_i$, defined as the ratio $\phi(x)/g(x)$ of the two densities up to normalization, hence we can take

$$w(x) = \exp(-x^2/2)(1+x^2),$$ but it could be safer to take the log-weights:

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

The weighted sample $(x_i,w(x_i))$ permite estimar la medida de cada intervalo $[u,v]$ debajo de la distribución objetivo, sumando los pesos de cada valor muestreado que cae dentro del intervalo:

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

Esto proporciona una estimación de la función acumulativa objetivo. Podemos obtenerlo y trazarlo rápidamente con el spatsatpaquete:

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 objetivo, sino de la distribución Cauchy instrumental, y uno obtiene una muestra de la distribución objetivo realizando un nuevo muestreo ponderado , 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: muestreo de transformación inversa rápida

Olver y Townsend desarrollaron un método de muestreo para una amplia clase de distribución continua. Se implementa en la biblioteca chebfun2 para Matlab , así como en la biblioteca ApproxFun para Julia . Recientemente descubrí esta biblioteca y suena muy prometedora (no solo para el muestreo aleatorio). Básicamente, este es el método de inversión, pero utiliza aproximaciones potentes del cdf y el cdf inverso. La entrada es la función de densidad objetivo hasta la normalización.

La muestra se genera simplemente mediante el siguiente código:

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

Como se verifica a continuación, produce una medida estimada del intervalo $[2,4]$ cercano al obtenido previamente por muestreo de importancia:

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