Ograniczone maksymalne prawdopodobieństwo z mniej niż pełną pozycją kolumny

To pytanie dotyczy oszacowania ograniczonego maksymalnego prawdopodobieństwa (REML) w określonej wersji modelu liniowego, a mianowicie:

Y = X (α) β + ϵ, ϵ \sim N_{n} (0, Σ (α)),

$Y = X(\alpha)\beta + \epsilon, \\ \epsilon\sim N_n(0, \Sigma(\alpha)),$

gdzie $X(\alpha)$ jest macierzą ( $n \times p$ ) sparametryzowaną przez $\alpha \in \mathbb R^k$ , podobnie jak $\Sigma(\alpha)$ . $\beta$ jest nieznanym wektorem parametrów uciążliwych; interesuje nas oszacowanie $\alpha$ , a mamy $k\leq p\ll n$ . Oszacowanie modelu według maksymalnego prawdopodobieństwa nie stanowi problemu, ale chcę użyć REML. Jest dobrze znane, patrz np. LaMotte , że prawdopodobieństwo $A'Y$ , gdzie $A$ jest dowolną półprostokątną macierzą taką, że można zapisać $A'X=0$

L_{REML} (α ∣ Y) \propto | X^{'} X |^{1 / 2} | Σ |^{- 1 / 2} | X^{'} Σ^{- 1} X |^{- 1 / 2} \exp {- \frac{1}{2} r^{'} Σ^{- 1} r}, r = (I - X (X^{'} Σ^{- 1} X)^{+} X^{'} Σ^{- 1}) Y,

$L_{\text{REML}}(\alpha\mid Y) \propto\vert X'X\vert^{1/2} \vert \Sigma\vert^{-1/2}\vert X'\Sigma^{-1}X\vert^{-1/2}\exp\left\{-\frac{1}{2} r'\Sigma^{-1}r \right\}, \\ r = (I - X(X'\Sigma^{-1}X)^+X'\Sigma^{-1})Y,$

gdy jest pełną pozycją kolumny $X$ .

Moim problemem jest to, że dla niektórych całkowicie rozsądnych i interesujących naukowo macierz nie ma pełnej rangi kolumny. Wszystkie pochodne, które widziałem powyżej o ograniczonym prawdopodobieństwie, wykorzystują wyznaczniki równości, które nie mają zastosowania, gdy , czyli zakładają pełną rangę kolumny . Oznacza to, że powyższe ograniczone prawdopodobieństwo jest poprawne tylko dla mojego ustawienia na częściach przestrzeni parametrów, a zatem nie jest tym, co chcę zoptymalizować. $\alpha$ $X(\alpha)$ $\vert X'X\vert=0$ $X$

Pytanie: Czy istnieją bardziej ogólne ograniczone prawdopodobieństwa, wyprowadzone w literaturze statystycznej lub gdzie indziej, bez założenia, że będzie pełną pozycją kolumny? Jeśli tak, to jak wyglądają? $X$

Niektóre spostrzeżenia:

Wyprowadzenie części wykładniczej nie stanowi problemu dla żadnego i można je zapisać w kategoriach odwrotności Moore-Penrose'a jak powyżej $X(\alpha)$
Kolumny stanowią (dowolną) ortonormalną podstawę dla $A$ $C(X)^\bot$
Dla znanego prawdopodobieństwo można łatwo zapisać dla każdego , ale oczywiście liczba wektorów bazowych, tj. Kolumn, w zależy od rangi kolumny $A$ $A'Y$ $\alpha$ $A$ $X$

Jeśli ktoś zainteresowany tym pytaniem uważa, że dokładna parametryzacja , daj mi znać, a ja je zanotuję. W tym momencie jednak najbardziej interesuje mnie REML dla ogólnego prawidłowych wymiarów. $X,\Sigma$ $X$

Bardziej szczegółowy opis modelu znajduje się tutaj. Niech będzie wymiarową autoregresją wektorową pierwszego rzędu [VAR (1)] gdzie $y_t = \mu + Ay_{t - 1} + v_t, t = 1, \dots, T$ $r$ . Załóżmy, że proces rozpoczyna się od pewnej stałej wartości w czasie . $v_t \overset{iid}{\sim}N(0, \Omega)$ $y_0$ $t = 0$

Zdefiniuj . Model można zapisać w postaci modelu liniowego przy użyciu następujących definicji i notacji: $Y = [y_1', \dots, y_T']'$ $Y = X\beta + \varepsilon$

\begin{aligned} X & = [1_{T} \otimes I_{r}, C^{- 1} B] \\ β & = [μ^{'}, y_{0}^{'} - μ^{'}]^{'} \\ v a r (ε)^{- 1} & = C^{'} (I_{T} \otimes Ω^{- 1}) C \\ C & = [\begin{matrix} I_{r} & 0 & 0 & \dots \\ - A & I_{r} & 0 & \dots \\ 0 & - A & I_{r} & \dots \\ ⋮ & ⋮ & ⋮ & ⋱ \end{matrix}] \\ B & = e_{1, T} \otimes A, \end{aligned}

$\begin{align} X &= [1_T \otimes I_r, C^{-1}B] \\ \beta &= [\mu', y_0' - \mu']' \\ \mathrm{var}(\varepsilon)^{-1} &= C'(I_T \otimes \Omega^{-1})C \\ C &= \begin{bmatrix} I_r & 0 & 0 & \cdots \\ -A & I_r & 0 & \cdots \\ 0 & -A & I_r & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{bmatrix} \\ B &= e_{1, T} \otimes A, \end{align}$

gdzie oznacza wymiarowy wektor zer i pierwszą standardową podstawę wektorowych . $1_T$ $T-$ $e_{1,T}$ $\mathbb R^T$

Oznacz . Zauważ, że jeśli nie jest pełną rangą, nie jest pełną rangą kolumny. Obejmuje to na przykład przypadki, w których jeden ze składników nie zależy od przeszłości. $\alpha = \mathrm{vec}(A)$ $A$ $X(\alpha)$ $y_t$

Pomysł oszacowania VAR przy użyciu REML jest dobrze znany, na przykład, w literaturze dotyczącej regresji predykcyjnej (patrz np. Phillips i Chen i odnośniki tam zawarte).

Warto wyjaśnić, że macierz nie jest macierzą projektową w zwykłym sensie, po prostu wypada z modelu i jeśli nie jest a priori wiedza na temat , o ile mogę powiedzieć, nie ma sposobu na ponowną parametryzację ma być pełna ranga. $X$ $A$

Na stronie mat.stackexchange zamieściłem pytanie, które jest z tym związane w tym sensie, że odpowiedź na pytanie matematyczne może pomóc w ustaleniu prawdopodobieństwa udzielenia odpowiedzi na to pytanie.

mixed-model maximum-likelihood linear-model optimization reml ekvall
źródło

Być może jednym ze sposobów rozwiązania tego pytania jest pytanie, co dzieje się w liniowych modelach mieszanych, gdy macierz modeli nie ma pełnej rangi kolumn?

Greenparker

Dzięki za nagrodę @Greenparker. I tak, gdyby można było zapisać ograniczone prawdopodobieństwo dla liniowego modelu mieszanego z matrycą projektową o ustalonych efektach mniejszą niż pełna kolumna, to by pomogło.

ekvall

Wyprowadzenie części wykładniczej nie stanowi problemu dla żadnego X (α) X (α) i można je zapisać w kategoriach odwrotności Moore-Penrose'a jak wyżej

Mam wątpliwości, czy ta obserwacja jest poprawna. Uogólniona inwersja faktycznie nakłada dodatkowe liniowe ograniczenie na twoje estymatory [Rao i Mitra], dlatego powinniśmy rozważyć wspólne prawdopodobieństwo jako całość zamiast zgadywać „Odwrotność Moore-Penrose'a będzie działała dla części wykładniczej”. Wydaje się to formalnie poprawne, ale prawdopodobnie nie rozumiesz poprawnie modelu mieszanego.

(1) Jak poprawnie myśleć o modelach z efektami mieszanymi? $\blacksquare$

Musisz pomyśleć model mieszanego efektu w inny sposób, zanim spróbujesz podłączyć odwrotność g (odwrotność Moore'a-Penrose'a, która jest specjalnym rodzajem zwrotnego odwrotności g [Rao i Mitra]) mechanicznie do formuły podanej przez RMLE (Restricted Estymator maksymalnego prawdopodobieństwa, to samo poniżej.).

X = (\begin{array}{cc} f i x e d e f f e c t \\ r a n d o m e f f e c t \end{array})

$\boldsymbol{X}=\left(\begin{array}{cc} fixed\quad effect\\ & random\quad effect \end{array}\right)$

Powszechnym sposobem myślenia efektu mieszanego jest to, że część losowego efektu w macierzy projektowej jest wprowadzana przez błąd pomiaru, który nosi inną nazwę „predyktor stochastyczny”, jeśli bardziej zależy nam na przewidywaniu niż na szacowaniu. Jest to także jedna historyczna motywacja badań macierzy stochastycznych w ustalaniu statystyk.

Mój problem polega na tym, że dla niektórych całkowicie rozsądnych i interesujących naukowo α α macierz X (α) X (α) nie ma pełnej rangi kolumny.

Biorąc pod uwagę ten sposób myślenia, prawdopodobieństwo, że nie jest pełnej rangi, wynosi zero. Wynika to z faktu, że funkcja wyznaczająca jest ciągła we wpisach macierzy, a rozkład normalny jest ciągłym rozkładem, który przypisuje zerowe prawdopodobieństwo pojedynczemu punktowi. Prawdopodobieństwo uszkodzenia rangi jest dodatnie, jeśli sparametryzowano ją w sposób patologiczny np. $X(\alpha)$ $X(\alpha)$ . $\left(\begin{array}{ccc} \alpha & \alpha\\ \alpha & \alpha\\ & & random\quad effect \end{array}\right)$

Więc rozwiązanie twojego pytania jest również dość proste, po prostu zaburzasz swoją macierz projektową (zaburza tylko część z efektem stałym) i używasz zaburzonej macierzy (która jest pełna ranga), aby przeprowadzić wszystkie pochodne. O ile twój model nie ma skomplikowanych hierarchii lub sam nie jest prawie osobliwy, nie widzę poważnego problemu, gdy weźmiesz w wyniku końcowym, ponieważ funkcja wyznacznika jest ciągła i możemy wziąć limit wewnątrz funkcji wyznacznika. $X_\epsilon(\alpha)=X(\alpha)+\epsilon\left(\begin{array}{cc} I & 0\\ 0 & 0 \end{array}\right)$ $X$ $\epsilon\rightarrow 0$ . W formie perturbacji odwrotność można uzyskać według twierdzenia Shermana-Morrisiona-Woodbury'ego. A wyznacznik macierzypodano w standardowej książce algebry liniowej, takiej jak [Horn & Johnson]. Oczywiście możemy zapisać wyznacznik w kategoriach każdego wpisu macierzy, ale zaburzenie jest zawsze preferowane [Horn & Johnson]. $lim_{\epsilon\rightarrow 0}|X_\epsilon|=|lim_{\epsilon\rightarrow 0}X_\epsilon|$ $X_\epsilon$ $I+X$

(2) Jak powinniśmy radzić sobie z uciążliwymi parametrami w modelu? $\blacksquare$

Jak widzisz, aby poradzić sobie z częścią efektu losowego w modelu, powinniśmy traktować ją jako rodzaj „uciążliwego parametru”. Problem polega na tym: czy RMLE jest najbardziej odpowiednim sposobem na wyeliminowanie uciążliwego parametru? Nawet w modelach GLM i efektach mieszanych RMLE nie jest jedynym wyborem. [Basu] wskazał, że istnieje wiele innych sposobów eliminacji parametrów przy ustalaniu oszacowania. Dzisiaj ludzie wybierają między RMLE a modelowaniem bayesowskim, ponieważ odpowiadają one dwóm popularnym rozwiązaniom komputerowym: odpowiednio EM i MCMC.

Moim zdaniem zdecydowanie bardziej odpowiednie jest wprowadzenie ołtarza w sytuacji wadliwej rangi w części ze stałym efektem. Możesz też sparametryzować swój model, aby uzyskać pełną rangę.

Ponadto, jeśli twój ustalony efekt nie ma pełnej rangi, możesz martwić się o źle określoną strukturę kowariancji, ponieważ stopnie swobody w ustalonych efektach powinny przejść do części błędu. Aby zobaczyć ten punkt jaśniej, może warto rozważyć również MLE (LSE) dla GLS (General najmniej gdzie jest struktura kowariancji warunek błędu, w przypadku, gdy nie jest pełnej rangi. $\hat{\beta}=(X\Sigma^{-1} X')^{-1}\Sigma^{-1}y$ $\Sigma$ $X(\alpha)$

(3) Dalsze komentarze $\blacksquare$

Problem nie polega na tym, jak zmodyfikujesz RMLE, aby działał w przypadku, gdy część matrycy o ustalonym efekcie nie ma pełnej rangi; problemem jest to, że w takim przypadku sam model może być problematyczny, jeśli przypadek niepełnej rangi ma dodatnie prawdopodobieństwo.

Jednym z istotnych przypadków, z którymi się spotkałem, jest to, że w przypadku przestrzennym ludzie mogą chcieć zmniejszyć rangę części o stałym efekcie ze względu na względy obliczeniowe [Wikle].

Nie widziałem żadnego „interesującego naukowo” przypadku w takiej sytuacji, czy mógłbyś wskazać literaturę, w której sprawa nie w pełni rangi budzi poważne obawy? Chciałbym wiedzieć i dyskutować dalej, dzięki.

Odniesienie $\blacksquare$

[Rao i Mitra] Rao, Calyampudi Radhakrishna i Sujit Kumar Mitra. Uogólniona odwrotność macierzy i jej zastosowania. Vol. 7. New York: Wiley, 1971.

[Basu] Basu, Debabrata. „W sprawie eliminacji uciążliwych parametrów”. Journal of American Statistics Association 72.358 (1977): 355-366.

[Horn & Johnson] Horn, Roger A. i Charles R. Johnson. Analiza macierzowa. Prasa uniwersytecka Cambridge, 2012.

[Wikle] Wikle, Christopher K. „Reprezentacje niskiej rangi dla procesów przestrzennych”. Handbook of Spatial Statistics (2010): 107-118.

Henry.L
źródło

Thanks for your interest and very thought through answer, + 1 for effort. I will read it in more detail and come back with some clarifications. I think a first thing that I will have to clarify is that there are no random effects in this model, and the matrix

X

$X$ is not a design matrix at all, except perhaps by name fr lack of a better word; it's a highly non-linear function (deterministic) of the parameter

α

$\alpha$ which consists of (the vectorization of) the coefficient matrix in a vector autoregressive process, so the concept of probability of being low-rank is not meaningful.

ekvall

@Student001 Yes, feel free to make any clarification since I also feel it more like a GLM instead of mixed model. I will try to answer again if I can:)

Henry.L

@Student001 If you can, do write the whole model and I would like to study such case, possibly AR(1) in spatial setting I guess.

Henry.L

"Given this way of thinking the likelihood, the probability that

X (α)

$X(\alpha)$ is not of full rank is zero." Right answer, wrong problem. The probability that it will be numerically not of full rank in finite precision is non-zero.

Mark L. Stone

@ MarkL.Stone Już podałem zaburzenie jako rozwiązanie, jeśli dokładnie czytasz wiersze, co jest standardowym rozwiązaniem osobliwości liczbowej. I OP powiedział, że zaktualizuje opis, więc myślę, że dojdziemy do konsensusu w sprawie poprawnie sformułowanego problemu.

Henry.L,

Ograniczone maksymalne prawdopodobieństwo z mniej niż pełną pozycją kolumny

Odpowiedzi: