O que significa distribuição truncada?

14

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 ± .

Kavka
fonte

Respostas:

16

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 densidadeN(μ,σ2)(uma,b)

puma,b(x)=ϕμ,σ2(x)umabϕμ,σ2(y)dyEu{x(uma,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 )ϕμ,σ2(x)N(μ,σ2)N(μ,σ2)(uma,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.±(uma,b)

Macro
fonte
17

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 )N(μ,σ2)(uma,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 / ϱ .

ϱ=umabφμ,σ2(x)dx
1/ϱ

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 .E(α)

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 )μ=0 0σ=1b=+E(α,uma) 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α(z)=αe-α(z-uma)Euzuma.
pa,(z)/gα(z)eα(za)ez2/2
exp(α2/2αa)α>aexp(a2/2) A primeira expressão é minimizada por α=1
{1/αexp(α2/2αa)if α>a,1/αexp(a2/2)otherwise.
enquanto ˜ α = a minimiza o segundo limite. A escolha ótima de α é, portanto, (1).
α=12a+12a2+4,(1)
α~=aα
Xi'an
fonte
2
vocêUnif(Φ(uma),Φ(b))X=Φ-1(você) ? Isso não dá a distribuição desejada?
bnaul
2
a0
Xian
1
Xi'an está certo, @bnaul. Rodar qnormem um loop R não é uma boa ideia.
Stéphane Laurent
@ Xi'an: Isso é verdade, mas essas funções podem ser projetadas para ter precisão arbitrária.
Neil G
9

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).

F1(U)FUUnif(0,1)FG(a,b)G1(U)UUnif(G(a),G(b))

G1G1GG1abG

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

N(0,1)GGG(q)=arctan(q)π+12G1(q)=tan(π(q12)). 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 UUnif(G(a),G(b)) and taking G1(U) is equivalent to take tan(U) with UUnif(arctan(a),arctan(b)):

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 xi, defined as the ratio ϕ(x)/g(x) of the two densities up to normalization, hence we can take

W(x)=exp(-x2/2)(1+x2),
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 (xEu,W(xEu)) permite estimar a medida de cada intervalo [você,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")

ewcdf

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

Claro, a amostra (xEu)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) 

hist

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
Stéphane Laurent
fonte