Co oznacza skrócona dystrybucja?

14

W artykule badawczym na temat analizy wrażliwości modelu równania różniczkowego zwyczajnego układu dynamicznego autor podał rozkład parametru modelu jako Rozkład normalny (średnia = 1e-4, std = 3e-5) obcięty do zakresu [0,5e -4 1,5e-4]. Następnie wykorzystuje próbki z tego obciętego rozkładu do symulacji modelu. Co to znaczy mieć obcięty rozkład i próbkę z tego obciętego rozkładu?

Mógłbym wymyślić dwa sposoby:

Próbkuj z rozkładu normalnego, ale przed symulacjami zignoruj wszystkie losowe wartości spoza określonego zakresu.
Jakoś uzyskać specjalny rozkład „Skrócona normalna” i pobrać z niego próbki.

Czy są to prawidłowe i równoważne podejścia?

Uważam, że w pierwszym przypadku, gdyby wykreślić eksperymentalny plik cdf / pdf próbki, nie wyglądałby on jak rozkład normalny, ponieważ krzywe nie rozciągają się do $\pm\infty$ .

distributions simulation truncation

— Kawka
źródło

16

Aby obciąć rozkład, należy ograniczyć jego wartości do przedziału i ponownie znormalizować gęstość, tak aby całka w tym zakresie wynosiła 1.

Tak więc obcięcie rozkładu $N(\mu, \sigma^{2})$ do przedziału $(a,b)$ oznaczałoby wygenerowanie zmiennej losowej o gęstości

p_{a, b} (x) = \frac{ϕ_{μ, σ^{2}} (x)}{\int_{a}^{b} ϕ_{μ, σ^{2}} (y) d y} \cdot I {x \in (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) \}$

gdzie $\phi_{\mu, \sigma^{2}}(x)$ to gęstość $N(\mu, \sigma^2)$ . Możesz próbkować z tej gęstości na wiele sposobów. Jednym ze sposobów (najprostszego, jaki mogę wymyślić), aby to zrobić, byłoby wygenerowanie wartości $N(\mu, \sigma^2)$ i wyrzucenie tych, które wypadają poza $(a,b)$ interwał, jak wspomniałeś. Tak, te dwie wymienione przez ciebie kule osiągnęłyby ten sam cel. Masz również rację, że gęstość empiryczna (lub histogram) zmiennych z tego rozkładu nie rozciągałaby się na . Oczywiście byłoby to ograniczone do . $\pm \infty$ $(a,b)$

— Makro
źródło

17

Symulowanie od normalnego rozkładu aż wynik mieści się w przedziale jest w porządku, gdy prawdopodobieństwo $\mathcal{N}(\mu,\sigma^2)$ $(a,b)$ jest wystarczająco duży. Jeśli jest zbyt mała, procedura ta jest zbyt kosztowna, ponieważ średnia liczba losowań dla jednej akceptacji wynosi .

ϱ = \int_{a}^{b} φ_{μ, σ^{2}} (x) d x

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

1 / ϱ

$1/\varrho$

Jak opisano w Metodach Statystycznych Monte Carlo (Rozdział 2, Przykład 2.2), a także w moim artykule arXiv , bardziej skutecznym sposobem symulacji tej okrojonej normy jest zastosowanie metody akceptacji-odrzucenia opartej na wykładniczym rozkładzie . $\mathcal{E}(\alpha)$

Rozważmy, bez utraty ogólności, przypadek i . Gdy , potencjalnym rozkładem instrumentalnym jest przetłumaczony rozkład wykładniczy o gęstości $\mu = 0$ $\sigma = 1$ $b=+\infty$ $\mathcal{E} (\alpha,{ a})$ Stosunek jest wówczas ograniczone jeśli i inaczej. Odpowiednia (górna) granica to

g_{α} (z) = α e^{- α (z - a)} I_{z \geq a} .

$g_{\alpha}(z) = \alpha e^{- \alpha(z - {a})} \; \mathbb{I}_{z \geq {a }} \;.$

p_{a, \infty} (z) / g_{α} (z) \propto e^{- α (z - a)} e^{- z^{2} / 2}

$p_{a,\infty}(z)/g_{\alpha}(z) \propto e^{- \alpha(z - a )}e^{-z^{2}/2}$

\exp (α^{2} / 2 - α a)

$\exp(\alpha^{2}/2 - \alpha{a })$

α > a

$\alpha > a$

\exp (- a^{2} / 2)

$\exp(- a^{2}/2)$

Pierwsze wyrażenie jest minimalizowane przez

{\begin{cases} 1 / α \exp (α^{2} / 2 - α a) & if α > a, \\ 1 / α \exp (- a^{2} / 2) & otherwise. \end{cases}

$\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}$

podczas gdy

minimalizuje drugą granicę. Optymalnym wyborem

jest zatem (1).

α^{*} = \frac{1}{2} a + \frac{1}{2} \sqrt{a^{2} + 4}, (1)

$\begin{equation} \alpha^{*} = \frac{1}{2}a + \frac{1}{2} \sqrt{a^2 + 4}\;,\qquad (1) \end{equation}$

\tilde{α} = a

$\tilde\alpha = a$

α

$\alpha$

— Xi'an
źródło

2

Być może czegoś mi brakuje, ale co jest nie tak z przyjmowaniem

i pozwalaniem

? Czy to nie daje pożądanej dystrybucji?

U \sim Unif (Φ (a), Φ (b))

$U \sim \text{Unif}(\Phi(a),\Phi(b))$

X = Φ^{- 1} (U)

$X = \Phi^{-1}(U)$

— bnaul

2

@bnaul: jest to całkowicie słuszne, ponieważ opiera się na odwrotnej transformacji cdf. Implikuje to jednak uzyskanie funkcji kwantylu rozkładu normalnego z bardzo dużą precyzją. Zwłaszcza, gdy

jest znacznie większe niż

.

a

$a$

0

$0$

— Xi'an,

1

Xi'an ma rację, @bnaul. Bieganie qnormw pętli R nie jest dobrym pomysłem.

— Stéphane Laurent

@ Xi'an: To prawda, ale takie funkcje mogą mieć dowolną precyzję.

— Neil G

9

Próbkuj z rozkładu normalnego, ale przed symulacjami zignoruj wszystkie losowe wartości spoza określonego zakresu.

Ta metoda jest poprawna, ale jak wspomniał @ Xi'an w swojej odpowiedzi, zajęłoby dużo czasu, gdy zasięg jest mały (a dokładniej, gdy jego miara jest mała w normalnym rozkładzie).

$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$ is a normal distribution, the evaluation of $G^{-1}$ is rather slow, and it is not highly precise for values of $a$ and $b$ outside the "range" of $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}),

$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))$ allows to estimate the measure of every interval $[u,v]$ under the target distribution, by summing the weights of each sampled value falling inside the interval:

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

This provides an estimate of the target cumulative function. We can quickly get and plot it with the spatsat package:

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

Of course, the sample $(x_i)$ is definitely not a sample of the target distribution, but of the instrumental Cauchy distribution, and one gets a sample of the target distribution by performing weighted resampling, for instance using the multinomial sampling:

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

hist

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

Another method: fast inverse transform sampling

Olver and Townsend developed a sampling method for a broad class of continuous distribution. It is implemented in the chebfun2 library for Matlab as well as the ApproxFun library for Julia. I have recently discovered this library and it sounds very promising (not only for random sampling). Basically this is the inversion method but using powerful approximations of the cdf and the inverse cdf. The input is the target density function up to normalization.

The sample is simply generated by the following code:

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

As checked below, it yields an estimated measure of the interval $[2,4]$ close to the one previously obtained by importance sampling:

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

— Stéphane Laurent
źródło