Szacowanie niepewności w problemach wnioskowania wielowymiarowego bez próbkowania?

Pracuję nad problemem wnioskowania o dużych wymiarach (około 2000 parametrów modelu), dla którego jesteśmy w stanie solidnie przeprowadzić oszacowanie MAP poprzez znalezienie globalnego maksimum log-tylnego przy użyciu kombinacji optymalizacji opartej na gradiencie i algorytmu genetycznego.

Bardzo chciałbym mieć możliwość oszacowania niepewności parametrów modelu oprócz znalezienia oszacowania MAP.

Jesteśmy w stanie efektywnie obliczyć gradient log-tylny w odniesieniu do parametrów, więc długoterminowo zamierzamy użyć Hamiltonian MCMC do wykonania próbkowania, ale na razie jestem zainteresowany szacunkami nieopartymi na próbkach.

Jedyne podejście, jakie znam, to obliczenie odwrotności Hesji w trybie w celu przybliżenia tylnej jako normalnej wielowymiarowej, ale nawet to wydaje się niemożliwe dla tak dużego układu, ponieważ nawet jeśli elementów Hesji jestem pewien, że nie mogliśmy znaleźć jego odwrotności. $\sim 4\times10^{6}$

Czy ktoś może zasugerować, jakie metody są zwykle stosowane w takich przypadkach?

Dzięki!

EDYCJA - dodatkowe informacje o problemie

Tło
To jest odwrotny problem związany z dużym eksperymentem fizyki. Mamy trójkątną siatkę 2D, która opisuje niektóre pola fizyczne, a naszymi parametrami modelu są fizyczne wartości tych pól na każdym wierzchołku siatki. Siatka ma około 650 wierzchołków, a my modelujemy 3 pola, stąd nasze parametry 2000 modeli.

Nasze dane eksperymentalne pochodzą z instrumentów, które nie mierzą bezpośrednio tych pól, ale wielkości, które są skomplikowanymi nieliniowymi funkcjami pól. Dla każdego z różnych instrumentów mamy model do przodu, który odwzorowuje parametry modelu na przewidywania danych eksperymentalnych, a porównanie prognozy i pomiaru daje logarytmiczne prawdopodobieństwo.

Następnie sumujemy prawdopodobieństwa dziennika z wszystkich tych różnych instrumentów, a także dodajemy wartości wcześniejszego dziennika, które nakładają pewne fizyczne ograniczenia na pola.

W związku z tym wątpię, by ten „model” należał do kategorii - nie mamy wyboru, czym jest model, jest on podyktowany tym, jak działają rzeczywiste instrumenty, które zbierają nasze dane eksperymentalne.

Zestaw
danych Zestaw danych składa się z 500 x 500 zdjęć, a dla każdej kamery jest jeden obraz, więc całkowita liczba punktów danych wynosi 500 x 500 x = . $10^6$

Model
błędów Obecnie bierzemy wszystkie błędy w tym problemie za Gaussa. W pewnym momencie mógłbym spróbować przejść do modelu błędu ucznia-t tylko dla pewnej dodatkowej elastyczności, ale wydaje się, że nadal działa dobrze tylko z Gaussianami.

Przykład prawdopodobieństwa
Jest to eksperyment z fizyką plazmy, a ogromna większość naszych danych pochodzi z kamer skierowanych na plazmę z określonymi filtrami przed obiektywami, aby patrzeć tylko na określone części spektrum światła.

Aby odtworzyć dane, należy wykonać dwa kroki; najpierw musimy zamodelować światło pochodzące z plazmy na siatce, a następnie modelować to światło z powrotem do obrazu z kamery.

Modelowanie światła pochodzącego z plazmy niestety zależy od tego, jakie są efektywne współczynniki szybkości, które mówią, ile światła jest emitowane przez różne procesy w danych polach. Stawki te są przewidywane przez niektóre drogie modele numeryczne, więc musimy przechowywać ich wyniki na siatkach, a następnie interpolować, aby wyszukać wartości. Dane funkcji szybkości obliczane są zawsze tylko raz - przechowujemy je, a następnie budujemy z nich splajn po uruchomieniu kodu, a następnie ten splajn jest wykorzystywany do wszystkich ocen funkcji.

Załóżmy, że i są funkcjami szybkości (które oceniamy przez interpolację), a następnie emisja w -tym wierzchołku siatki jest podana przez gdzie to 3 pola, które modelujemy na siatce. Przeniesienie wektora emisji na obraz z kamery jest łatwe, wystarczy pomnożenie przez macierz która koduje części siatki, przez które patrzy każdy piksel kamery. $R_1$ $R_2$ $i$ $\mathcal{E}_i$

E_{i} = R_{1} (x_{i}, y_{i}) + z_{i} R_{2} (x_{i}, y_{i})

$\mathcal{E}_i = R_1(x_i, y_i) + z_i R_2(x_i, y_i)$

(x, y, z)

$(x,y,z)$

G

$\mathbf{G}$

Ponieważ błędy są gaussowskie, prawdopodobieństwo dziennika dla tego konkretnego aparatu wynosi zatem

L = - \frac{1}{2} (G \vec{E} - \vec{d})^{⊤} Σ^{- 1} (G \vec{E} - \vec{d})

$\mathcal{L} = -\frac{1}{2} (\mathbf{G}\vec{\mathcal{E}} - \vec{d})^{\top}\mathbf{\Sigma}^{-1} (\mathbf{G}\vec{\mathcal{E}} - \vec{d})$

gdzie to dane kamery. Całkowite prawdopodobieństwo logarytmu jest sumą 4 powyższych wyrażeń, ale dla różnych kamer, z których wszystkie mają różne wersje funkcji szybkości ponieważ patrzą na różne części spektrum światła. $\vec{d}$ $R_1, R_2$

Wcześniejszy przykład
Mamy różne priorytety, które skutecznie po prostu ustalają pewne górne i dolne granice dla różnych wielkości, ale te zwykle nie działają zbyt silnie na problem. Mamy jeden wcześniej działający silnie, który skutecznie stosuje wygładzanie typu Laplaciana na polach. Przybiera również postać gaussowską:

log-prior = - \frac{1}{2} {\vec{x}}^{⊤} S \vec{x} - \frac{1}{2} {\vec{y}}^{⊤} S \vec{y} - \frac{1}{2} {\vec{z}}^{⊤} S \vec{z}

$\text{log-prior} = -\frac{1}{2}\vec{x}^{\top}\mathbf{S}\vec{x} -\frac{1}{2}\vec{y}^{\top}\mathbf{S}\vec{y} -\frac{1}{2}\vec{z}^{\top}\mathbf{S}\vec{z}$

— CBowman
źródło

Jaki model pasujesz? Regresja liniowa? GP? Hierarchiczny model zliczania? Bayesowska kalibracja modelu komputerowego? Dodaj więcej szczegółów na temat rozwiązanego problemu, a ja napiszę odpowiedź za i przeciw VI.

— DeltaIV

@DeltaIV Zaktualizowałem pytanie o więcej informacji - być może nie opracowałem dokładnie tego, czego szukałeś. Jeśli tak, daj mi znać, a ja dokonam kolejnej edycji, dzięki!

— CBowman

@DeltaIV Jeszcze raz dziękuję! Dodano więcej informacji, daj mi znać, jeśli mogę coś jeszcze dodać.

— CBowman

@DeltaIV obrazy danych mają wymiary 500 x 500 i jest jeden dla każdej kamery, więc łączna liczba punktów danych wynosi 500 x 500 x = . Dane funkcji szybkości obliczane są zawsze tylko raz - przechowujemy je, a następnie budujemy z nich splajn po uruchomieniu kodu, a następnie ten splajn jest wykorzystywany do wszystkich ocen funkcji.

10^{6}

$10^6$

— CBowman

Nie mam odniesienia, ale istnieje wiele przybliżeń niskiej rangi do obliczania macierzy odwrotnej. np znaleźć największe eigenvalues załóżmy, pozostając są równe i używać szorstkiej Przybliżony wektorów własnych odpowiadających niskiej wartości własnej. Jestem pewien, że istnieją również przybliżone / iteracyjne dekompozycje Cholesky'ego, które są zbieżne z dokładną wartością. po prostu zakończ iteracje po tym, jak czekasz na maksymalny czas

k

$k$

2000 - k

$2000-k$

— probabilityislogic

Odpowiedzi:

Po pierwsze, myślę, że twój model statystyczny jest zły. Zmieniam więc waszą notację na bardziej znaną statystykom, niech więc

d = y = (y_{1}, \dots, y_{N}), N = 10^{6}

$\mathbf{d}=\mathbf{y}=(y_1,\dots,y_N),\ N=10^6$

być wektorem obserwacji (danych) i

\begin{aligned} x & = θ = (θ_{1}, \dots, θ_{p}) \\ y & = ϕ = (ϕ_{1}, \dots, ϕ_{p}) \\ z & = ρ = (ρ_{1}, \dots, ρ_{p}), p \approx 650 \end{aligned}

$\begin{align} \mathbf{x}&=\boldsymbol{\theta}=(\theta_1,\dots,\theta_p) \\ \mathbf{y}&=\boldsymbol{\phi}=(\phi_1,\dots,\phi_p) \\ \mathbf{z}&=\boldsymbol{\rho}=(\rho_1,\dots,\rho_p), \ p \approx 650 \\ \end{align}$

twoje wektory parametrów o wymiarze całkowitym . . Następnie, jeśli dobrze zrozumiałem, zakładasz model $d=3p \approx 2000$

y = G r_{1} (θ, ϕ) + ρ G r_{2} (θ, ϕ)) + ϵ, ϵ \sim N (0, I_{N})

$\mathbf{y} = \mathbf{G}\mathbf{r_1}(\boldsymbol{\theta}, \boldsymbol{\phi})+\boldsymbol{\rho}\mathbf{G}\mathbf{r_2}(\boldsymbol{\theta}, \boldsymbol{\phi}))+\boldsymbol{\epsilon},\ \boldsymbol{\epsilon}\sim\mathcal{N}(0,I_N)$

gdzie jest macierzą interpolacji splajnu . $\mathbf{G}$ $N\times d$

To jest oczywiście złe. Nie ma możliwości, aby błędy w różnych punktach obrazu z tej samej kamery, aw tym samym punkcie na obrazach z różnych kamer, były niezależne. Powinieneś przyjrzeć się statystykom przestrzennym i modelom, takim jak uogólnione najmniejsze kwadraty, estymacja semiwariogramu, kriging, procesy Gaussa itp.

To powiedziawszy, ponieważ twoje pytanie nie dotyczy tego, czy model jest dobrym przybliżeniem faktycznego procesu generowania danych, ale jak oszacować taki model, pokażę ci kilka opcji, aby to zrobić.

HMC

2000 parametrów nie jest bardzo dużym modelem, chyba że trenujesz to na laptopie. Zestaw danych jest większy ( punktów danych), ale mimo to, jeśli masz dostęp do wystąpień w chmurze lub maszyn z procesorami graficznymi, frameworki takie jak Pyro lub Tensorflow Prawdopodobieństwo szybko rozwiążą taki problem. Tak więc możesz po prostu użyć zasilanego GPU Hamiltoniana Monte Carlo. $10^6$

Plusy : wnioskowanie „dokładne” w granicach nieskończonej liczby próbek z łańcucha.

Wady : brak ścisłego ograniczenia błędu estymacji, istnieje wiele wskaźników diagnostycznych zbieżności, ale żaden nie jest idealny.

Przybliżenie dużej próbki

Z nadużywaniem notacji, niech Oznaczmy przez wektora otrzymanego przez złączenie swoje trzy wektory parametrów. Następnie, korzystając z centralnego twierdzenia Bayesa o ograniczeniach (Bernstein-von Mises), możesz aproksymować pomocą , gdzie jest wartością parametru „true”, jest oszacowaniem MLE dla a to matryca informacji Fishera oceniana przy . Oczywiście, jest nieznany, użyjemy $\theta$ $p(\theta\vert \mathbf{y})$ $\mathcal{N}(\hat{\theta_0}_n,I_n^{-1}(\theta_0))$ $\theta_0$ $\hat{\theta_0}_n$ $\theta_0$ $I_n^{-1}(\theta_0)$ $\theta_0$ $\theta_0$ $I_n^{-1}(\hat{\theta_0}_n)$ zamiast. Ważność twierdzenia Bernsteina-von Misesa zależy od kilku hipotez, które można znaleźć, np. Tutaj : w twoim przypadku, zakładając, że są płynne i zróżnicowane, twierdzenie jest ważne, ponieważ poparcie Gaussa Prior to cała przestrzeń parametrów. Lub, lepiej, byłoby prawidłowe, gdyby twoje dane były rzeczywiście takie, jak zakładasz, ale nie sądzę, że są, jak wyjaśniłem na początku. $R_1,R_2$

Zalet : szczególnie użyteczne w przypadków. Gwarantujemy zbieżność z właściwą odpowiedzią, w ustawieniu iid, gdy prawdopodobieństwo jest płynne i zróżnicowane, a przeor jest niezerowy w sąsiedztwie . $p<<N$ $\theta_0$

Minusy : największym problemem, jak zauważyłeś, jest potrzeba odwrócenia matrycy informacji Fishera. Poza tym nie wiedziałbym, jak empirycznie ocenić dokładność aproksymacji, poza użyciem próbnika MCMC do pobierania próbek z . Oczywiście, to w pierwszej kolejności zniweczyłoby użyteczność B-vM. $p(\theta\vert \mathbf{y})$

Wnioskowanie wariacyjne

W tym przypadku, zamiast znalezienie dokładnego (co wymaga obliczenia wyniku idimensional całkowitej), zdecydowaliśmy się w przybliżeniu z , gdzie należy do rodziny parametrycznej indeksowanej przez wektor parametrów . Szukamy st . Zminimalizowano pewną miarę rozbieżności między i . Wybierając tę miarę jako dywergencję KL, otrzymujemy metodę wnioskowania wariacyjnego: $p(\theta\vert \mathbf{y})$ $d-$ $p$ $q_{\phi}(\theta)$ $q$ $\mathcal{Q}_{\phi}$ $\phi$ $\phi^*$ $q$ $p$

ϕ^{*} = \underset{ϕ \in Φ}{a r g m i n} D_{K L} (q_{ϕ} (θ) | | p (θ | y))

$\DeclareMathOperator*{\argmin}{arg\,min} \phi^*=\argmin_{\phi\in\Phi}D_{KL}(q_{\phi}(\theta)||p(\theta\vert\mathbf{y}))$

Wymagania dotyczące : $q_{\phi}(\theta)$

powinno być możliwe do rozróżnienia względem , abyśmy mogli zastosować metody optymalizacji na dużą skalę, takie jak Stochastic Gradient Descent, aby rozwiązać problem minimalizacji. $\phi$
powinien być wystarczająco elastyczny, aby mógł dokładnie przybliżać dla pewnej wartości , ale także wystarczająco prosty, aby łatwo było z niego próbkować. Jest tak, ponieważ oszacowanie dywergencji KL (nasz cel optymalizacji) wymaga oszacowania oczekiwań wr . $p(\theta\vert\mathbf{y})$ $\phi$ $q$

Możesz wybrać aby być w pełni podzielonym na czynniki, tj. Iloczyn jednoznacznych rozkładów prawdopodobieństwa: $q_{\phi}(\theta)$ $d$

q_{ϕ} (θ) = \prod_{i = 1}^{d} q_{ϕ_{i}} (θ_{i})

$q_{\phi}(\theta)=\prod_{i=1}^d q_{\phi_i}(\theta_i)$

jest to tak zwana metoda Bayesa wariacyjnego o średnim polu . Można udowodnić (patrz np. Rozdział 10 tej książki ), że optymalnym rozwiązaniem dla każdego z czynników jest $q_{\phi_j}(\theta_j)$

\log q_{j}^{*} (θ_{j}) = E_{i \neq j} [\log p (y, θ)] + const.

$\log{q_j^*(\theta_j)} = \mathbb{E}_{i\neq j}[\log{p(\mathbf{y},\theta)}] + \text{const.}$

gdzie jest wspólnym rozkładem parametrów i danych (w twoim przypadku jest to iloczyn twojego prawdopodobieństwa Gaussa i gaussowskiego pierwszeństwa względem parametrów), a oczekiwanie dotyczy innej wariacji dystrybucje jednowymiarowe . Oczywiście, ponieważ rozwiązanie jednego z czynników zależy od wszystkich innych czynników, musimy zastosować procedurę iteracyjną, inicjując wszystkie dystrybucje celu wstępnego odgadnięcia, a następnie iteracyjnie aktualizując je jeden jednocześnie z powyższym równaniem. Zauważ, że zamiast obliczać powyższe oczekiwania jako $p(\mathbf{y},\theta)$ $q_1^*(\theta_1),\dots,q_{j-1}^*(\theta_{j-1}),q_{j+1}^*(\theta_{j+1}),\dots,q_{d}^*(\theta_{d})$ $q_{i}(\theta_{i})$ $(d-1)-$ całka wymiarowa, która byłaby zaporowa w twoim przypadku, gdy priorytety i prawdopodobieństwo nie są sprzężone, możesz użyć oszacowania Monte Carlo, aby oszacować oczekiwanie.

Algorytm zmiennych pól zmiennych w polu średniej nie jest jedynym możliwym algorytmem VI, którego można użyć: zmiennokształtny autoenkoder zaprezentowany w Kingma & Welling, 2014, „Automatyczne kodowanie zmiennych wariacyjnych” jest interesującą alternatywą, w której zamiast przyjąć w pełni złożoną formę dla , a następnie wyprowadzając wyrażenie w postaci zamkniętej dla , zakłada się , że jest wielowymiarowym gaussowskim, ale z prawdopodobnie różnymi parametrami w każdym z punktów danych. Aby zamortyzować koszt wnioskowania, sieć neuronowa służy do mapowania przestrzeni wejściowej na przestrzeń parametrów wariacyjnych. Szczegółowy opis algorytmu znajduje się w artykule: implementacje VAE są ponownie dostępne we wszystkich głównych platformach głębokiego uczenia. $q$ $q_i$ $q$ $N$

— DeltaIV
źródło

ten model niezależności VB może być okropnym podejściem do pomiarów dokładności . Zwykle jest to przybliżenie typu wtyczki bez regulacji. proste przykłady nie są za pomocą „stopni swobody” korekt w ty i stosując zamiast normalnego rozkładu t. szczególnie problem hiper parametrów

s^{2}

$s^2$

— prawdopodobieństwo

@DeltaIV Model statystyczny jest ogólnie całkiem niezły, błędy między różnymi kamerami są bardzo bardzo niezależne, a różne piksele w tej samej kamerze będą w zasadzie niezależne, chyba że dosłownie sąsiadują ze sobą. Możemy zakodować pewną korelację przestrzenną w sąsiednich pikselach, używając prawdopodobieństwa procesu Gaussa, ale wymagałoby to albo bezpośredniego odwrócenia macierzy kowariancji, albo rozwiązania rzadkiego układu liniowego za każdym razem, gdy chcemy ocenić prawdopodobieństwo, co jest znacznie więcej drogie (choć nie wykluczone).

— CBowman

możesz chcieć sprawdzić niektóre oprogramowanie „bayesX”, a być może także oprogramowanie „inla”. oba mogą mieć kilka pomysłów, które możesz wypróbować. wygoogluj to

oba polegają bardzo mocno na wykorzystaniu rzadkości w parametryzacji macierzy dokładności (tj. niezależność warunkowa, model typu Markowa) - i mają zaprojektowane do tego algorytmy inwersji. większość przykładów opiera się na wielopoziomowych lub auto-regresywnych modelach guassian. powinien być dość podobny do opublikowanego przykładu

— prawdopodobieństwo prawdopodobieństwa
źródło