Czy to jest poprawne ? (generowanie skróconej normy-wielowymiarowej-gaussowskiej)

Jeśli $X\in\mathbb{R}^n,~X\sim \mathcal{N}(\underline{0},\sigma^2\mathbf{I})$ tj.

$$f_X(x) = \frac{1}{{(2\pi\sigma^2)}^{n/2}} \exp\left(-\frac{||x||^2}{2\sigma^2}\right)$$

Chcę analogicznej wersji skróconego rozkładu normalnego w przypadku wielu odmian.

Dokładniej, chcę wygenerować ograniczony normą (do wartości $\geq a$ ) wielowymiarowy Gaussowski $Y$ st

$$f_Y(y) = \begin{cases} c.f_X(y), \text{ if } ||y||\geq a \\[2mm] 0, \text{ otherwise }. \end{cases}$$ gdzie $c=\frac{1}{Prob\big\{||X||\geq a\big\}}$

---

Teraz obserwuję:

Jeśli $x=(x_1,x_2,\ldots,x_n)$ ,$||x||\geq a$

$\implies |x_n|\geq T\triangleq \sqrt{\max\left(0,\left(a^2-\sum_1^{n-1}x_i^2\right)\right)}$

Dlatego wybierając jako próbki Gaussa, można ograniczyć jako próbkę poza rozkładem obciętym -normalnym (zgodnie z rozkładem ogona Gaussa ) , z wyjątkiem jego losowo wybranego znaku z prawdopodobieństwem . x n ≥ T N T ( 0 , σ 2 ) 1 /$x_1,\ldots,x_{n-1}$$x_n$$\geq T$$\mathcal{N}_T(0,\sigma^2)$$1/2$

Teraz moje pytanie brzmi:

Jeśli wygeneruję każdą próbkę wektorową z( X 1 , … , X n$(x_1,\ldots,x_n)$$(X_1,\ldots,X_n)$ jako,

$x_1,\ldots,x_{n-1}\sim \mathcal{N}(0,\sigma^2)$

i

$x_n = Z_1 *Z_2~$ gdzie, , , (tj. ścięty-skalar-normalny RV zZ 2 ~ N T ( 0 , σ 2 ) t ( x 1 , ... , x n - 1 ) ≜ $~Z_1\sim\{\pm1 ~\text{w.p.}~ 1/2\}$$Z_2\sim\mathcal{N}_T(0,\sigma^2)$$T(x_1,\ldots,x_{n-1})\triangleq \sqrt{\max\left(0,\left(a^2-\sum_1^{n-1}x_i^2\right)\right)}$

Czy będzie normatywnym ( ) wielowymiarowym gaussowskim? (tj. taki sam jak zdefiniowany powyżej ). Jak powinienem zweryfikować? Wszelkie inne sugestie, jeśli tak nie jest?≥ a $(X_1,X_2,\ldots,X_n)$$\geq a$$Y$

EDYTOWAĆ:

Oto wykres punktowy punktów w przypadku 2D z normą obciętą do wartości powyżej „1”

Uwaga: Poniżej znajduje się kilka świetnych odpowiedzi, ale brakuje uzasadnienia, dlaczego ta propozycja jest błędna. W rzeczywistości to główny punkt tego pytania.

Wielowymiarowy rozkład normalny jest sferycznie symetryczny. Poszukiwany rozkład obcina promień ρ = | | X | | 2 poniżej o . Ponieważ to kryterium zależy tylko od długości X , skrócony rozkład pozostaje sferycznie symetryczny. Ponieważ ρ jest niezależne od sferycznego kąta X / | | X | | i $X$$\rho=||X||^2$$a$$X$$\rho$$X/||X||$ marozkład χ ( n ) , dlatego możesz wygenerować wartości z rozkładu obciętego w kilku prostych krokach:$\rho\,\sigma$$\chi(n)$

1. Wygeneruj .$X \sim \mathcal{N}(0,\mathbb{I}_n)$

2. Wygeneruj jako pierwiastek kwadratowy rozkładu χ 2 ( d ) obciętego w ( a / σ ) 2 .$P$$\chi^2(d)$$(a/\sigma)^2$

3. Niech .$Y = \sigma P\, X/||X||$

W kroku 1 jest uzyskiwany jako sekwencja d niezależnych realizacji standardowej zmiennej normalnej.$X$$d$

W etapie 2, jest łatwo wytwarzany przez odwracanie odwrotna dystrybuanta F - 1 o × 2 ( d ) rozkład: wygenerować jednorodną zmienna U obsługiwane w zakresie (kwantyli) pomiędzy F ( ( / Ď ) 2 ) i 1 i ustaw P = $P$$F^{-1}$$\chi^2(d)$$U$$F((a/\sigma)^2)$$1$ .$P = \sqrt{F(U)}$

Oto histogram takich niezależnych realizacji σ P dla σ = 3 w n = 11 wymiarów, skrócony poniżej przy a = 7 . Wygenerowanie zajęło około sekundy, co świadczy o wydajności algorytmu.$10^5$$\sigma P$$\sigma=3$$n=11$$a=7$

Czerwona krzywa jest gęstością skróconego rozkładu skalowanego przez σ = 3 . Jego ścisłe dopasowanie do histogramu świadczy o ważności tej techniki.$\chi(11)$$\sigma=3$

Aby uzyskać intuicję dotyczącą obcięcia, rozważ przypadek , σ = 1 in n = 2 wymiary. Oto wykres rozrzutu Y 2 w stosunku do Y 1 (dla 10 4 niezależnych realizacji). Wyraźnie pokazuje otwór w promieniu a :$a=3$$\sigma=1$$n=2$$Y_2$$Y_1$$10^4$$a$

Na koniec zauważ, że (1) komponenty muszą mieć identyczne rozkłady (ze względu na symetrię sferyczną) i (2), z wyjątkiem sytuacji, gdy a = 0 , ten wspólny rozkład nie jest normalny. W rzeczywistości, jak rośnie duża, szybki spadek (jednoczynnikowa) rozkładu normalnego powoduje, że większość prawdopodobieństwa sferycznie obcięty wielowymiarowa normalne klaster w pobliżu powierzchni z n - 1 -sphere (o promieniu A ). Dystrybucja krańcowa musi zatem być zbliżona do skalowanej symetrycznej wersji Beta ( ( n - 1 ) / 2 , ( n $X_i$$a=0$$a$$n-1$$a$. rozkład skoncentrowany w przedziale ( - a , a ) . Jest to widoczne w poprzednim wykresie rozrzutu, w którym a = 3 σ jest już duży w dwóch wymiarach: punkty otaczają pierścień (sfera 2 - 1 ) o promieniu 3 $((n-1)/2,(n-1)/2)$$(-a,a)$$a=3\sigma$$2-1$$3\sigma$

Oto histogramy krańcowych z symulacji rozkładu wielkości w 3 wymiarach z a = 10 , σ = 1 (dla których aproksymującego beta ( 1 , 1 ) Rozkład jest jednolita)$10^5$$3$$a=10$$\sigma=1$$(1,1)$

Ponieważ pierwsze marginesów procedury opisanej w pytaniu są normalne (z założenia), procedura ta nie może być poprawna.$n-1$

---

Poniższy Rkod wygenerował pierwszą cyfrę. Jest on zbudowany z równoległych etapów 1-3 do wytwarzania . Został zmodyfikowany w celu wytworzenia drugiego postać zmieniając zmiennych , , i i wydający polecenia działki po$Y$adnsigmaplot(y[1,], y[2,], pch=16, cex=1/2, col="#00000010")y został wygenerowany.

Generowanie jest modyfikowane w kodzie w celu uzyskania wyższej rozdzielczości numerycznej: kod faktycznie generuje 1 - U i używa go do obliczenia $U$$1-U$$P$ .

Tę samą technikę symulowania danych zgodnie z domniemanym algorytmem, podsumowania ich za pomocą histogramu i nałożenia histogramu można zastosować do przetestowania metody opisanej w pytaniu. Potwierdzi to, że metoda nie działa zgodnie z oczekiwaniami.

a <- 7      # Lower threshold
d <- 11     # Dimensions
n <- 1e5    # Sample size
sigma <- 3  # Original SD
#
# The algorithm.
#
set.seed(17)
u.max <- pchisq((a/sigma)^2, d, lower.tail=FALSE)
if (u.max == 0) stop("The threshold is too large.")
u <- runif(n, 0, u.max)
rho <- sigma * sqrt(qchisq(u, d, lower.tail=FALSE)) 
x <- matrix(rnorm(n*d, 0, 1), ncol=d)
y <- t(x * rho / apply(x, 1, function(y) sqrt(sum(y*y))))
#
# Draw histograms of the marginal distributions.
#
h <- function(z) {
  s <- sd(z)
  hist(z, freq=FALSE, ylim=c(0, 1/sqrt(2*pi*s^2)),
       main="Marginal Histogram",
       sub="Best Normal Fit Superimposed")
  curve(dnorm(x, mean(z), s), add=TRUE, lwd=2, col="Red")
}
par(mfrow=c(1, min(d, 4)))
invisible(apply(y, 1, h))
#
# Draw a nice histogram of the distances.
#
#plot(y[1,], y[2,], pch=16, cex=1/2, col="#00000010") # For figure 2
rho.max <- min(qchisq(1 - 0.001*pchisq(a/sigma, d, lower.tail=FALSE), d)*sigma, 
               max(rho), na.rm=TRUE)
k <- ceiling(rho.max/a)
hist(rho, freq=FALSE, xlim=c(0, rho.max),  
     breaks=seq(0, max(rho)+a, by=a/ceiling(50/k)))
#
# Superimpose the theoretical distribution.
#
dchi <- function(x, d) {
  exp((d-1)*log(x) + (1-d/2)*log(2) - x^2/2 - lgamma(d/2))
}
curve((x >= a)*dchi(x/sigma, d) / (1-pchisq((a/sigma)^2, d))/sigma, add=TRUE, 
      lwd=2, col="Red", n=257)

Napisałem to, zakładając, że nie chcesz, aby jakiekolwiek punkty posiadały || y || > a, który jest odpowiednikiem zwykłego obcięcia jednowymiarowego. Jednak napisałeś, że chcesz zachować punkty posiadające | y || > = a i wyrzuć pozostałe. Niemniej jednak można dokonać oczywistej korekty mojego rozwiązania, jeśli naprawdę chcesz zachować punkty mające | y || > = a.

Najprostszym sposobem, który okazuje się być bardzo ogólną techniką, jest użycie odrzucenia akceptacji https://en.wikipedia.org/wiki/Rejection_sampling . Będzie dość szybki, dopóki Prob (|| X ||> a) będzie dość niski, ponieważ wtedy nie będzie wielu odrzuceń.

Wygeneruj wartość próbki x z nieograniczonej normalnej zmiennej wielowymiarowej (nawet jeśli twój problem mówi, że normalna wielowymiarowa jest sferyczna, można zastosować tę technikę, nawet jeśli nie jest). Jeśli || x || <= a, zaakceptuj, tzn. użyj x, w przeciwnym razie odrzuć go i wygeneruj nową próbkę. Powtarzaj ten proces, aż uzyskasz tyle zaakceptowanych próbek, ile potrzebujesz. Skutkiem zastosowania tej procedury jest wygenerowanie y w taki sposób, że jego gęstość wynosi c * f_X (y), jeśli || y || <= a, a 0 jeśli || y || > a, według mojej poprawki do części otwierającej pytania. Nigdy nie musisz obliczać c; jest to automatycznie określane przez algorytm na podstawie częstotliwości odrzucania próbek.

To niezła próba, ale nie działa z powodu „stałej normalizacji”: jeśli weźmie się pod uwagę gęstość połączenia rozkład fX(x)∝

$$f_X(x) \propto \frac{1}{{(2\pi\sigma^2)}^{n/2}} \exp\left(-\frac{||x||^2}{2\sigma^2}\right)\mathbb{I}_{||x||>a}=\frac{1}{{(2\pi\sigma^2)}^{n/2}} \exp\left(-\frac{x_1^2+\ldots+x_n^2}{2\sigma^2}\right)\mathbb{I}_{||x||>a}$$ $$f_X(x) \propto \frac{1}{{(2\pi\sigma^2)}^{(n-1)/2}} \exp\left(-\frac{||x_{-n}||^2}{2\sigma^2}\right)\frac{1}{{(2\pi\sigma^2)}^{1/2}} \exp\left(-\frac{x_n^2}{2\sigma^2}\right)\mathbb{I}_{||x||>a}$$ $$=\frac{1}{{(2\pi\sigma^2)}^{(n-1)/2}} \exp\left(-\frac{||x_{-n}||^2}{2\sigma^2}\right)\frac{1}{{(2\pi\sigma^2)}^{1/2}} \exp\left(-\frac{x_n^2}{2\sigma^2}\right)\mathbb{I}_{||x_{-n}||^2+x_n^2>a^2}$$ $$=\frac{\mathbb{P}(X_n^2>a^2-||x_{-n}||^2)}{{(2\pi\sigma^2)}^{(n-1)/2}} \exp\left(-\frac{||x_{-n}||^2}{2\sigma^2}\right)\qquad\qquad\qquad\qquad\qquad$$ $$\qquad\qquad\qquad\times\frac{\mathbb{P}(X_n^2>a^2-||x_{-n}||^2)^{-1}}{{(2\pi\sigma^2)}^{1/2}} \exp\left(-\frac{x_n^2}{2\sigma^2}\right)\mathbb{I}_{x_n^2>a-||x_{-n}||^2}$$ który integruje się z $$f_{X_{-n}}(x_{-n}) \propto \frac{\mathbb{P}(X_n^2>a^2-||x_{-n}||^2)}{{(2\pi\sigma^2)}^{(n-1)/2}} \exp\left(-\frac{||x_{-n}||^2}{2\sigma^2}\right)$$ w $x_n$, pokazuje, że

1. Rozkład warunkowy $X_n$ biorąc pod uwagę inne elementy, $X_{-n}$, jest obciętym rozkładem normalnym;
2. Rozkład krańcowy pozostałych składników, $X_{-n}$, nie jest rozkładem normalnym z powodu dodatkowego terminu$\mathbb{P}(X_n^2>a^2-||x_{-n}||^2)$;

Jedynym sposobem, jaki widzę w korzystaniu z tej właściwości, jest uruchomienie samplera Gibbsa, jednego komponentu na raz, przy użyciu skróconych normalnych rozkładów warunkowych.

Pytanie pochodzi z pomysłu użycia - podstawowego rozkładu warunkowego rozkładów połączeń - w celu narysowania próbek wektorowych.

Pozwolić $X$ być wielowymiarowym gaussowskim z komponentami iid.

Pozwolić $\text{Prob}(||X||>a) \triangleq T$ i $Y\triangleq X.\mathbb{I}_{||X||>a}$

Algorytm, o którym mowa, zaproponowano w oparciu o następującą (całkowicie poprawną, ale wprowadzającą w błąd interpretację) faktoryzację warunkową:

$$f_Y(y) = \frac{1}{T}\frac{1}{{(2\pi\sigma^2)}^{n/2}} \exp\left(-\frac{||y||^2}{2\sigma^2}\right)\mathbb{I}_{||y||>a}\\ =\frac{1}{T}\frac{1}{{(2\pi\sigma^2)}^{n/2}} \exp\left(-\frac{y_1^2+\ldots+y_n^2}{2\sigma^2}\right)\mathbb{I}_{||y||>a}\\ =\left(\prod_{i=1}^{n-1} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{y_i^2}{2\sigma^2}\right)\right) \left(\frac{1}{T}\frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{y_n^2}{2\sigma^2}\right)\mathbb{I}_{||y||>a}\right)\\ =\underbrace{\left(\prod_{i=1}^{n-1} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{y_i^2}{2\sigma^2}\right)\right)}_{\text{Gaussians}} \underbrace{\left(\frac{1}{T}\frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{y_n^2}{2\sigma^2}\right)\mathbb{I}_{y_n^2>(a^2-y_1^2-\ldots y_{n-1}^2)}\right)}_{\text{Truncated Gaussian??}}$$

Najkrótsza odpowiedź jest taka, że ​​ten ostatni czynnik nie jest obciętym gaussowskim, (co ważniejsze) nawet rozkładem.

---

Oto szczegółowe wyjaśnienie, dlaczego sama powyższa faktoryzacja ma jakąś zasadniczą wadę. W jednym zdaniu: każda warunkowa faktoryzacja danego rozkładu połączeń musi spełniać pewne bardzo fundamentalne właściwości, a powyższa faktoryzacja ich nie spełnia (patrz poniżej).

Zasadniczo, jeśli kiedykolwiek uwzględnimy czynniki $f_{XY}(x,y)=f_X(x)\cdot f_{Y|X}(y|x)$ następnie $f_X(x)$ jest na marginesie $X$ i $f_{Y|X}(y|x)$ jest rozkładem warunkowym $Y$. Co znaczy:

1. Współczynnik $f(x,y)$ „zakładany jako” $f_X(x)$musi być dystrybucją. I,
2. Drugi czynnik „przyjęty jako” $f_{Y|X}(y|x)$musi być rozkładem dla każdego wyboru$x$

W powyższym przykładzie próbujemy warunkować jako $Y_n|(Y_1\ldots Y_{n-1})$. Oznacza to, że właściwość-1 powinna obejmować czynnik Gaussa, a właściwość-2 powinna utrzymywać się dobrze w drugiej części.

Oczywiste jest, że właściwość-1 dobrze trzyma się pierwszego czynnika. Ale problem dotyczy właściwości-2. Ostatni czynnik powyżej niestety nie jest wcale rozkładem (zapomnij o skróconym gaussowskim) dla prawie dowolnej wartości$(Y_1\ldots Y_{n-1})$!!

---

Taka propozycja algorytmu jest prawdopodobnie wynikiem następującego błędnego przekonania: Gdy rozkład naturalnie odłączy się od wspólnego rozkładu (takiego jak powyżej Gaussa), prowadzi to do faktoryzacji warunkowej. ---- Nie ma! ---- Drugi (drugi) czynnik musi być również dobry.

---

Uwaga: Istnieje świetna odpowiedź @whuber wcześniej, która faktycznie rozwiązuje problem generowania normalnego skróconego wielowymiarowego Gaussa. Przyjmuję jego odpowiedź. Ta odpowiedź ma jedynie na celu wyjaśnienie i podzielenie się moim własnym zrozumieniem i genezą pytania.