O que significa distribuição truncada?

Em um artigo de pesquisa sobre análise de sensibilidade de um modelo de equação diferencial ordinária de um sistema dinâmico, o autor forneceu a distribuição de um parâmetro de modelo como Distribuição normal (média = 1e-4, std = 3e-5) truncada no intervalo [0,5e -4 1,5e-4]. Ele então usa amostras dessa distribuição truncada para simulações do modelo. O que significa ter uma distribuição truncada e uma amostra dessa distribuição truncada?

Eu poderia criar duas maneiras de fazer isso:

• Amostra de uma distribuição Normal, mas ignore todos os valores aleatórios que estão fora do intervalo especificado antes das simulações.
• De alguma forma, obtenha uma distribuição especial "Truncated Normal" e obtenha amostras dela.

Essas abordagens são válidas e equivalentes?

Acredito que no primeiro caso, se alguém plotar o cdf / pdf experimental da amostra, não pareceria uma distribuição Normal porque as curvas não se estendem a $\pm\infty$ .

Truncar uma distribuição é restringir seus valores a um intervalo e re-normalizar a densidade para que a integral nesse intervalo seja 1.

Portanto, truncar a distribuição para um intervalo ( a , b ) seria gerar uma variável aleatória que tenha densidade$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) \}$$

onde é a densidade de N ( μ , σ 2 ) . Você pode amostrar dessa densidade de várias maneiras. Uma maneira (a maneira mais simples que eu consigo pensar) de fazer isso seria gerar valores de N ( μ , σ 2 ) e jogar fora os que ficam fora dos ( a , b $\phi_{\mu, \sigma^{2}}(x)$$N(\mu, \sigma^2)$$N(\mu, \sigma^2)$$(a,b)$intervalo, como você mencionou. Então, sim, essas duas balas que você listou atingiriam o mesmo objetivo. Além disso, você está certo de que a densidade empírica (ou histograma) das variáveis ​​dessa distribuição não se estenderia a . Seria restrito a ( a , b ) , é claro.$\pm \infty$$(a,b)$

Simular a distribuição normal de até que o resultado caia dentro de um intervalo ( a , b ) é bom quando a probabilidade ϱ = ∫ b a φ μ , σ 2 ( x $\mathcal{N}(\mu,\sigma^2)$$(a,b)$ é grande o suficiente. Se for muito pequeno, esse procedimento é muito caro, pois o número médio de empates para uma aceitação é 1 / ϱ .

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

Como descrito nos Métodos Estatísticos de Monte Carlo (Capítulo 2, Exemplo 2.2), bem como no meu artigo do arXiv , uma maneira mais eficiente de simular esse normal truncado é usar um método de aceitação / rejeição baseado em uma distribuição exponencial .$\mathcal{E}(\alpha)$

Considere, sem perda de generalidade, o caso e σ = 1 . Quando b = + ∞ , uma distribuição instrumental potencial é a distribuição exponencial traduzida, E ( α , a ) , com densidade g α ( z ) = α e - α ( z - a$\mu = 0$$\sigma = 1$$b=+\infty$$\mathcal{E} (\alpha,{ a})$ A proporção p um , ∞ ( z ) / g α ( z ) α e - α ( z - um ) e - z 2 / 2 é, em seguida, delimitada por exp ( α 2 / 2 - α um ) se α > um e pela exp ( - um 2 / 2 ) de outra forma. O limite (superior) correspondente é

$$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)$ A primeira expressão é minimizada por α∗= $$\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}$$ enquanto ˜ α = a minimiza o segundo limite. A escolha ótima de α é, portanto, (1). $$\begin{equation} \alpha^{*} = \frac{1}{2}a + \frac{1}{2} \sqrt{a^2 + 4}\;,\qquad (1) \end{equation}$$ $\tilde\alpha = a$$\alpha$

Amostra de uma distribuição Normal, mas ignore todos os valores aleatórios que estão fora do intervalo especificado antes das simulações.

Esse método está correto, mas, como mencionado por @ Xi'an em sua resposta, levaria muito tempo quando o intervalo é pequeno (mais precisamente, quando sua medida é pequena na distribuição normal).

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

$G^{-1}$$G^{-1}$$G$$G^{-1}$$a$$b$$G$

Simule uma distribuição truncada usando amostragem de importância

${\cal N}(0,1)$$G$$G$$\boxed{G(q)=\frac{\arctan(q)}{\pi}+\frac12}$$\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),$$ mas poderia ser mais seguro usar os pesos do log:

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

A amostra ponderada $(x_i,w(x_i))$ permite estimar a medida de cada intervalo $[u,v]$ na distribuição de destino, somando os pesos de cada valor amostrado dentro do intervalo:

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

Isso fornece uma estimativa da função cumulativa de destino. Podemos rapidamente obtê-lo e plotá-lo com o spatsatpacote:

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

Claro, a amostra $(x_i)$definitivamente não é uma amostra da distribuição de destino, mas da distribuição instrumental de Cauchy, e obtém-se uma amostra da distribuição de destino realizando uma reamostragem ponderada , por exemplo, usando a amostra multinomial:

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

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

Outro método: amostragem por transformação inversa rápida

Olver e Townsend desenvolveram um método de amostragem para uma ampla classe de distribuição contínua. Ele é implementado na biblioteca chebfun2 para Matlab , bem como na biblioteca ApproxFun para Julia . Eu descobri recentemente essa biblioteca e parece muito promissora (não apenas para amostragem aleatória). Basicamente, este é o método de inversão, mas usando aproximações poderosas do cdf e do cdf inverso. A entrada é a função de densidade de destino até a normalização.

A amostra é simplesmente gerada pelo seguinte código:

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

Como verificado abaixo, ele produz uma medida estimada do intervalo $[2,4]$ próximo ao obtido anteriormente por amostragem por importância:

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