---
fontsize: 8pt
bibliography: 5_Referenzen.bib
citation_package: natbib
output:
  beamer_presentation:
    keep_tex: true
    includes:
      in_header: 5_header.tex
---


```{r, include = F}
source("5_R_common.R")
```

#  {.plain}
\center
```{r, echo = FALSE, out.width = "20%"}
knitr::include_graphics("5_Abbildungen/mvda_5_otto.png")
```

\vspace{2mm}

\Huge
Multivariate Datenanalyse
\vspace{6mm}

\large
MSc Psychologie WiSe 2021/22

\vspace{6mm}
\large
Prof. Dr. Dirk Ostwald

#  {.plain}

\vfill
\center
\huge
\textcolor{black}{(5) Faktorenanalyse}
\vfill

# Überblick
\small
\justifying
Generative Perspektive zu konfirmatorischer und exploratorischer Faktoranalyse

* "Generativ" bedeutet hier probabilistisch und modell-basiert.

Einführung zum Expectation-Maximization (EM) Algorithmus.

* Generischer Algorithmus zur Schätzung von Parametern in Modellen mit latenten Variablen.

* Moderne Sichtweise von EM als Evidende Lower Bound Maximierung.

* Ein Schritt in Richtung eines Verständnisses von Variational Inference.


\begin{itemize}
\justifying
\item[$\Rightarrow$] Integrative Perspektive von Inferenz und Lernen in Modellen mit
latenten Variablen mit einer natürlichen Generalisierung zu kontemporären Modellen
des maschinellen Lernens und der künstlicher Intelligenz (Variational Bayesian Filtering
oder ``Variational Autoencoders'' für die Generation Deep Learning)
\item[$\Rightarrow$]Ein Schritt zu einem Verständnis zeitgenössischer Theorien
zur Funktionsweise des Gehirns (Free Energy Principle, Active Inference,
Agent-based behavioral models)
\end{itemize}


# Überblick
\vfill
```{r, echo = FALSE, out.width = "100%"}
knitr::include_graphics("5_Abbildungen/mvda_5_wissenschaft.pdf")
```
\vfill

# Literaturhinweise

\small
\vfill
Wir folgen der Darstellung Linearer Normalverteilungsmodelle in @roweis_1999.
@dempster_1977 bietet eine inhaltsreiche Einführung zum EM Algorithmus. Die
Anwendung des EM Algorithmus im Rahmen der Faktorenanalyse geht auf @rubin_1982
zurück. Probabilistische Hauptkomponentenanalye wird in @tipping_1999
und @roweis_1998 diskutiert und geht auf Arbeiten von @lawley_1953 zurück. Ursprünglich
wurde die Hauptkomponentenanalyse von @pearson_1901 vorgeschlagen und insbesondere
von @hotelling_1933 verfeinert. Das Anwendungsbeispiel aus dem Gebiet der
kognitiven Fähigkeitsforschung geht auf das Beispiel zur konfirmatorischen
Faktorenanalyse in @rosseel_2012 basierend auf @joreskog_1969 und @holzinger_1939.
\vfill

#
\large
\setstretch{2}
\vfill
Grundlagen

Inferenz und Lernen

Faktorenanalyse

Hauptkomponentenanalyse Revisited

Anwendungsbeispiel

Selbstkontrollfragen

Appendix
\vfill

#
\large
\setstretch{2}
\vfill
**Grundlagen**

Inferenz und Lernen

Faktorenanalyse

Hauptkomponentenanalyse Revisited

Anwendungsbeispiel

Selbstkontrollfragen

Appendix
\vfill


# Grundlagen
\setstretch{1.3}
\small
\begin{definition}[Multivariate Normalverteilung]
\justifying
$X$ sei ein $n$-dimensionaler Zufallsvektor mit Ergebnisraum $\mathbb{R}^n$ und WDF
\begin{equation}
p : \mathbb{R}^n \to \mathbb{R}_{>0},\, x \mapsto p(x)
:= (2\pi)^{-\frac{n}{2}}|\Sigma|^{-\frac{1}{2}}\exp\left(-\frac{1}{2}(x-\mu)^T \Sigma^{-1} (x-\mu)\right).
\end{equation}
Dann sagen wird, dass $X$ einer \textit{multivariaten (oder $n$-dimensionalen)
Normalverteilung} mit \textit{Erwartungswertparameter} $\mu \in \mathbb{R}^n$ und
positive-definitem \textit{Kovarianzmatrixparameter} $\Sigma \in \mathbb{R}^{n \times n}$
unterliegt und nennen $X$ einen \textit{(multivariat) normalverteilten Zufallsvektor}.
Wir kürzen dies mit $X \sim N(\mu,\Sigma)$ ab. Die WDF eines multivariat
normalverteilten Zufallsvektors bezeichnen wir mit
\begin{equation}
N(x;\mu,\Sigma):= (2\pi)^{-\frac{n}{2}}|\Sigma|^{-\frac{1}{2}}\exp\left(-\frac{1}{2}(x-\mu)^T \Sigma^{-1} (x-\mu)\right).
\end{equation}
\end{definition}

\footnotesize
Bemerkungen

* Der Parameter $\mu \in \mathbb{R}^n$ entspricht dem Wert höchster Wahrscheinlichkeitsdichte
* Die Diagonalelemente von $\Sigma$ spezifizieren die Breite der WDF bezüglich $X_1,...,X_n$.
* Das $i,j$te Element von $\Sigma$ spezifiziert die Kovarianz on $X_i$ und $X_j$.
* Der Term $(2\pi)^{-\frac{n}{2}}|\Sigma|^{-\frac{1}{2}}$ ist die Normalsieriungskonstante für den Exponentialfunktionsterm.

# Grundlagen
\small
\textcolor{darkblue}{Zweidimensionale Normalverteilungen}
\vspace{2mm}

```{r, echo = FALSE, out.width = "100%"}
knitr::include_graphics("5_Abbildungen/mvda_5_normalverteilungen.pdf")
```

# Grundlagen
\small
\begin{definition}[Lineares Normalverteilungsmodell]
\justifying
$X$ sei ein kontinuierlicher nicht-beobachtbarer $k$-dimensionaler Zufallsvektor
und $Y$ sei ein kontinuierlicher beobachbarer $m$-dimensionaler Zufallsvektor.
$B \in \mathbb{R}^{m \times k}$ sei eine Matrix und $R \in \mathbb{R}^{m \times m}$
sei eine positiv-definite Matrix. Dann heißt ein probabilistisches Modell mit WDF
\begin{equation}
p(x,y) = p(y|x)p(x),
\end{equation}
wobei
\begin{equation}
p(y|x) := N(y;Bx,R)
\mbox{ und }
p(x) := N(x;0_k, I_k)
\end{equation}
gilt, ein \textit{lineares Normalverteilungsmodell (LNM)}. Die Parametermenge eines
LNMs ist $\theta := \{B,R\}$ und wir schreiben die WDFen von LNMen allgemein als
\begin{equation}
p_\theta(x,y) := N(y;Bx,R)N(x;0_k,I_k).
\end{equation}
\end{definition}


# Grundlagen
Bemerkungen

\small
* Der Zufallsvektor $X$ heißt auch *Zustandsvektor* oder *latenter Vektor*
* Der Zufallsvektor $Y$ heißt auch *Datenvektor*.
* In hierarchischer Form kann ein LNM geschrieben werden als
\begin{align}
\begin{split}
X & = \xi			\quad\quad\quad\quad\,	\xi  \sim N(0_k,I_k) \\
Y & = BX + \eta 	\quad 						\eta \sim N(0_m,R).
\end{split}
\end{align}

* $\xi$ heißt dabei *Zustandsrauschen*, $\eta$ heißt dabei *Beobachtungsrauschen/fehler*.
* Samplen eines LNMs resultiert in Realisierungen $(x^{(i)}, y^{(i)})$ mit $i = 1,...,n$.
* Die $x^{(i)} \in \mathbb{R}^k$ modellieren nicht beobachtbare/latente/virtuelle Daten.
* Die $y^{(i)} \in \mathbb{R}^m$ modellieren beobachtbare Daten.
* LNMe sind spezielle lineare normalverteilte Zustandsraummodelle.


# Grundlagen
\small
\begin{theorem}[LNM Datenverteilung]
\justifying
\normalfont
Die Datenverteilung des LNMs
\begin{equation}
p_\theta(x,y) = N(y;Bx,R)N(x;0_k,I_k)
\end{equation}
ist gegeben durch
\begin{equation}
p_\theta(y) = N(y;0_m, BB^T + R)
\end{equation}
und wird auch als \textit{marginale Datenverteilung} oder \textit{marginale Likelihood}
bezeichnet.
\end{theorem}


\footnotesize
\underline{Beweis}
\vspace{1mm}

Wir halten zunächst fest, dass mit dem Theorem zu Gemeinsamen Normalverteilungen
aus (3) Wahrscheinlichkeitstheorie direkt folgt, dass die WDF der gemeinsamen
Verteilung von $X$ und $Y$ gegeben ist durch
\begin{equation}
p_\theta(x,y) =
N\left(
\begin{pmatrix} x 		\\ y 	\end{pmatrix} ;
\begin{pmatrix} 0_k 	\\ 0_m 	\end{pmatrix},
\begin{pmatrix}
I_k		& B^T \\
B 		& BB^T + R
\end{pmatrix}
\right).
\end{equation}
Mit dem Theorem zu Marginalen Normalverteilungen aus (3) Wahrscheinlichkeitstheorie
folgt dann aber soforft, dass die WDF der marginalen Verteilung von $Y$ durch
\begin{equation}
p_\theta(y) = N\left(y; 0_m, BB^T + R \right).
\end{equation}
gegeben ist.
$\hfill \Box$


# Grundlagen
\setstretch{2}
\small
Bemerkungen

* LNMs modellieren zentrierte multivariat normalverteilte Datensätze mit Kovarianzmatrix.
\begin{equation}
\mathbb{C}(Y) =  BB^T + R.
\end{equation}
* "Modellieren" bedeutet hier insbesondere, die Datenkovarianzmatrix $\mathbb{C}(Y)$ zu erklären.
* Die Form von $\mathbb{C}(Y)$ resultiert dabei aus der Transformation einer latenten Normalverteilung.
* Die Matrizen $B$ und $R$ generieren/erklären $\mathbb{C}(Y)$ mechanistisch.
* $B$ und $R$ ermöglichen so oft eine kondensiertere Erklärung als $\mathbb{C}(Y)$ \textit{per se}.
* Verschiedene LNMe haben unterschiedliche Potentiale, Datenkovarianzmatrizen zu erklären.

# Grudlagen
```{r, echo = FALSE, out.width = "100%"}
knitr::include_graphics("5_Abbildungen/mvda_5_lnm.pdf")
```

\footnotesize
\setstretch{1.6}

* Der latente Zufallsvektor $X$ ist  in $\mathbb{R}^k$ (hier $k = 2$) spärisch verteilt.
* Durch $B$ wird diese Sphäre gedehnt, rotiert, und nach $\mathbb{R}^m$ (hier $m = 2$) transformiert.
* Bei $m < p$  sieht die Sphäre von $X$ zum Beispiel wie ein Pfannkuchen aus.
* Dieser Pfannkuchen wird dann noch mit der Kovarianz von des Beobachtungsrauschens $\eta$ konvolviert.
* Es mag helfen, sich diese Konvolution (Faltung) als "Addition von Rauschen" vorzustellen.


# Grundlagen
\setstretch{1.8}
\textcolor{darkblue}{Spezielle lineare Normalverteilungsmodelle}

\small
\begin{itemize}
\item Das Ziel der Datenmodellierung mit LNMe ist die Erklärung der Datenkovarianmatrixstruktur.
\item Die Datenkovarianmatrixstruktur kann durch Wahl von $B$ und $R$ erklärt werden
\item Spezielle LNMe entsprechen spezifischen Randbedingungen für $R$.
\item[] $\Rightarrow$ In der konfirmatorischen Faktorenanalyse wird $R$ als Diagonalmatrix vorausgesetzt.
\item[] $\Rightarrow$ In der probabilistischen PCA wird $R$ als sphärisch vorausgesetzt.
\item[] $\Rightarrow$ In der exploratorischen Faktorenanalyse (PCA) wird $R = 0_{mm}$ vorausgesetzt.
\end{itemize}

#
\large
\setstretch{2}
\vfill
Grundlagen

**Inferenz und Lernen**

Faktorenanalyse

Hauptkomponentenanalyse Revisited

Anwendungsbeispiel

Selbstkontrollfragen

Appendix

# Inferenz und Lernen
\setstretch{1.2}
\footnotesize
\begin{definition}[Verteilungen von LNM Datensätzen]
\justifying
\begin{equation}
Y := \begin{pmatrix}y^{(1)} & \cdots & y^{(n)} \end{pmatrix} \in \mathbb{R}^{m \times n}
\mbox{ und }
X := \begin{pmatrix}x^{(1)} & \cdots & x^{(n)} \end{pmatrix} \in \mathbb{R}^{k \times n}
\end{equation}
seien ein beobachteter Datensatz (Realisierungen des beobachtbaren Zufallsvektors)
und der assoziierte unbeobachtete Datensatz (Realisierungen des latenten Zufallsvektors).
Unter der Annahme unabhängier und identischer Verteilung der gemeinsamen Realisationen
$(x^{(i)},y^{(i)})$ für $i = 1,...,n$ ist die WDF der gemeinsamen Verteilung
eines LNM Datensatz durch
\begin{equation}
p_\theta(X,Y)
= \prod_{i=1}^n p_\theta\left(x^{(i)}, y^{(i)}\right)
= \prod_{i=1}^n N\left(y^{(i)}; Bx^{(i)}, R\right)N\left(x^{(i)};0_k,I_k\right)
\end{equation}
und die marginale WDF des beobachtbaren Datensatzes gegeben durch
\begin{equation}
p_\theta(Y)
= \prod_{i=1}^n N\left(y^{(i)}; 0_m, BB^T + R \right).
\end{equation}
\end{definition}

\small
Bemerkungen

* **$X$ und $Y$ bezeichnen ab jetzt keine Zufallsvektoren mehr.**
* **$X \in \mathbb{R}^{k \times n}$ und $Y \in \mathbb{R}^{m \times n}$ bezeichnen ab jetzt Matrizen.**

# Inferenz und Lernen
\justifying
\textcolor{darkblue}{Inferenz}

\small
Was ist die Verteilung der latenten Zufallsvektoren und was sind ihre
wahrscheinlichsten Werte für feste Werte der LNM Parameter $\theta := \{B,R\}$?

$\Rightarrow$ Die Antwort gibt das Theorem zu bedingten multivariaten Normalverteilungen.
\vspace{3mm}

\normalsize
\textcolor{darkblue}{Lernen}

\small
Welche Parameterwerte maximieren die Marginal-Likelihood Funktion
\begin{equation}
L : \Theta \to \mathbb{R}_{\ge 0}, \theta \mapsto L(\theta)
:= p_\theta(Y) = \int p_\theta(X,Y)\,dX,
\end{equation}
oder, äquivalent, die Log Marginal-Likelihood Funktion
\begin{equation}
\ell : \Theta \to \mathbb{R}, \theta \mapsto \ell(\theta)
:= \ln p_\theta(Y) =  \ln \int p_\theta(X,Y)\,dX \mbox{ ?}
\end{equation}
für einen festen beobachteten Datensatz $Y \in \mathbb{R}^{m \times n}$?

$\Rightarrow$ Die Antwort gibt der Expectation-Maximization Algorithmus.


# Inferenz und Lernen
\small
\begin{theorem}[LNM Inferenz]
\justifying
\normalfont
Es sei
\begin{equation*}
p_\theta(x,y) = N(y;Bx,R)N(x;0,I_k),
\end{equation*}
ein LNM. Dann ist die WDF der bedingten Verteilung des latenten Zufallsvektors
gegeben durch
\begin{equation}
p_\theta(x|y) = N\left(x; B^T(BB^T + R)^{-1}y, I_k - B^T(BB^T + R)^{-1}B\right).
\end{equation}
\end{theorem}
\footnotesize
Bemerkungen

* Die bedingte Verteilung des latenten Zufallsvektors ist eine Normalverteilung.
* Der wahrscheinlichste Wert des latenten ZVs gegeben eine Beobachtung des beobachtbaren ZVs ist also
\begin{equation}
\hat{x} := \mu_{x|y} = B^T(BB^T + R)^{-1}y.
\end{equation}
* Die mit diesem Wert assoziierte Unsicherheit ist $\Sigma_{x|y} := I_k - B^T(BB^T + R)^{-1}B$.

# Inferenz und Lernen
\footnotesize
\underline{Beweis}
\vspace{1mm}

Wir hatten oben bereits gesehen, dass die WDF der gemeinsamen Verteilung von
latentem und beobachtbarem Zufallsvektor gegeben ist durch
\begin{equation}
p_\theta(x,y) =
N\left(
\begin{pmatrix} x 	  \\ y \end{pmatrix} ;
\begin{pmatrix} 0_k 	\\ 0_m \end{pmatrix},
\begin{pmatrix}
I_k		& B^T \\
B 		& BB^T + R
\end{pmatrix}
\right)
\end{equation}
Mit dem Theorem zu Bedingten Normalverteilungen aus Einheit (3) Wahrscheinlichkeitstheorie
gilt dann mit Identifikation von
\begin{equation}
\mu_x := 0_k,
\mu_y := 0_m,
\Sigma_{xx} := I_k,
\Sigma_{xy} := B^T,
\Sigma_{yx} := B, \mbox{ und }
\Sigma_{yy} := BB^T + R,
\end{equation}
dass
\begin{equation}
p_\theta(x|y) = N\left(x; \mu_{x|y},\Sigma_{x|y}\right),
\end{equation}
wobei
\begin{equation}
\mu_{x|y}
= \mu_x + \Sigma_{xy}\Sigma_{yy}^{-1}(y - \mu_y)
= B^T(BB^T + R)^{-1}y,
\end{equation}
und wobei
\begin{equation}
\Sigma_{x|y}
= \Sigma_{xx} -  \Sigma_{xy}\Sigma_{yy}^{-1}\Sigma_{yx}
= I_k - B^T(BB^T + R)^{-1}B.
\end{equation}
ist.
$\hfill \Box$

# Inferenz und Lernen
\small
\begin{theorem}[LNM Datensatzinferenz]
\justifying
\normalfont
$p_\theta(X,Y)$ sei die gemeinsame Verteilung eines LNM Datenatzes. Dann ist die
WDF der bedingten Vertilung von $X$ gegeben $Y$ gegeben durch
\begin{equation}
p_\theta(X|Y) =
\prod_{i=1}^n N\left(x^{(i)}; B^T(BB^T + R)^{-1}y^{(i)}, I_k - B^T(BB^T + R)^{-1}B\right).
\end{equation}
\end{theorem}

\footnotesize
\underline{Beweis}
\vspace{1mm}

Mit den Ausdrücken für die WDFen der gemeinsamen und marginalen LNM Datensatzverteilungen gilt
\begin{align*}
p_\theta(X|Y)
= \frac{p_\theta(X,Y)}{p_\theta(Y)}
= \frac{\prod_{i=1}^n p_\theta(x^{(i)},y^{(i)})}{\prod_{i=1}^n p_\theta\left(y^{(i)} \right)}
= \prod_{i=1}^n\frac{p_\theta(x^{(i)},y^{(i)})}{p_\theta(y^{(i)})}
= \prod_{i=1}^n p_\theta\left(x^{(i)}|y^{(i)}\right).
\end{align*}
Das Theorem folgt dann direkt mit dem Theorem zur LNM Inferenz.

$\hfill \Box$

# Inferenz und Lernen
\footnotesize
\begin{theorem}[Evidence Lower Bound]
\justifying
\normalfont
Für einen Datensatz $Y\in \mathbb{R}^{m\times n}$ sei $\ln p_\theta(Y)$ die WDF
der Log Marginal-Likelihood Verteilung eines LNMs. Dann gilt für jede WDF $q(X)$,
dass
\begin{align*}
\ln p_\theta(Y)
\ge  \int q(X) \ln p_\theta(X,Y) \,dX - \int q(X) \ln q(X) \,dX
=: \mbox{ELBO}(q(X),\theta).
\end{align*}
$\mbox{ELBO}(q(X),\theta)$ heißt \textit{Evidence Lower Bound}.
\end{theorem}

\setstretch{1.0}
\underline{Beweis}

Mit der Jensenschen Ungleichung (Appendix) gilt
\begin{equation*}
\ln p_\theta(Y)
:= \ln \int p_\theta(X,Y)\,dX
 = \ln \int q(X) \left(\frac{p_\theta(X,Y)}{q(X)}\right) \,dX
\ge \int q(X) \ln \left(\frac{p_\theta(X,Y)}{q(X)} \right)\,dX.
\end{equation*}
Damit aber folgt
\begin{align*}
\ln p_\theta(Y)
 \ge 	\int q(X) \ln \left(\frac{p_\theta(X,Y)}{q(X)} \right)\,dX
& =  	\int q(X) \left(\ln p_\theta(X,Y) - \ln q(X) \right)\,dX  	\\
& = 	\int q(X) \ln p_\theta(X,Y) \,dX - \int q(X) \ln q(X) \,dX.
\end{align*}


# Inferenz und Lernen
\setstretch{1.8}
Bemerkungen
\small

* Für einen festen Datensatz $Y$ ist $\mbox{ELBO}(q(X),\theta)$ eine Funktion von $q(X)$ und $\theta$.
* Die Bezeichnung Evidence Lower Bound geht auf den Term "Evidence" für $p_\theta(Y)$ zurück.
* In den kognitiven Neurowissenschaften ist die ELBO als "Freie Energie" bekannt.
* Die Signifikanz der ELBO geht weit über LNMe hinaus:
  * \small \justifying Die ELBO ist für die Variational Inference zentral.
  * Variational Inference ist für moderne Theorien zur Funktion des Gehirns zentral.
* Für Einführungen zur Variational Inference siehe zum Beispiel
  * \small @ostwald_2014, @starke_2017, @blei_2017.


# Inferenz und Lernen
\small
\begin{definition}[Expectation-Maximization Algorithmus]
\justifying
Die iterative koordinatenweise Maximierung der ELBO hinsichtlich $q(X)$ und $\theta$
heißt \textit{Expectation-Maximization (EM) Algorithmus}. Der Algorithmus
hat die allgemeine Form
\vspace{2mm}

\underline{EM Algorithmus}
\vspace{1mm}

0. Initialisierung von $q^{(0)}(X)$ und $\theta^{(0)}$
\vspace{1mm}

Für $j = 1,2,...$
\begin{itemize}
\item[1.] E Schritt: Setze  $q^{(j)}(X) := \argmax_{q(X)} \mbox{ELBO}\left(q(X),\theta^{(j-1)}\right)$
\item[2.] M Schritt: Setze  $\,\,\,\,\,\,\,\,\theta^{(j)} := \argmax_{\theta} \mbox{ELBO}\left(q^{(j)}(X),\theta\right)$
\end{itemize}
Nach Konvergenz, nutze $\hat{\theta} := \theta^{(j)}$ als Parameterschätzer.
\end{definition}

\footnotesize
Bemerkungen

* "Expectation Schritt" ist eine Fehlbezeichnung, es handelt sich auch um einen Maximization Schritt...
* ... allerdings ergibt die Bezeichnung im sogenannten "exakten" EM Algorithmus Sinn.

# Inferenz und Lernen
\small
\begin{theorem}[Exakter Expectation-Maximization Algorithmus]
\justifying
\normalfont
Das Setzen von
\begin{equation}
q^{(j)}(X) := p_{\theta^{(j-1)}}(X|Y) \mbox{ für alle } j = 1,2,...
\end{equation}
im E Schritt des EM Algorithmus maximiert die ELBO hinsichtlich $q(X)$ und
heißt \textit{exakter E Schritt}. Der Algorithmus hat dann die folgende Form

\vspace{2mm}

\underline{Exakter EM Algorithmus}
\vspace{1mm}

0. Initialisierung von  $\theta^{(0)}$
\vspace{1mm}

Für $j = 1,2,...$
\begin{itemize}
\item[1.] E Step
$q^{(j)}(X) := p_{\theta^{(j-1)}}(X|Y)$
\item[2.] M Step
$\,\,\,\,\,\,\,\,\theta^{(j)} := \argmax_{\theta} \int p_{\theta^{(j-1)}}(X|Y) \ln p_\theta(X,Y) \,dX$
\end{itemize}
Nach Konvergenz, nutze $\hat{\theta} := \theta^{(j)}$ als Parameterschätzer.
\end{theorem}

\footnotesize
Bemerkungen

* Für LNMe kann $p_{\theta^{(j-1)}}(X|Y)$ analytisch evaluiert werden $\Rightarrow$ Inferenz.
* Der M Schritt des exakten EM Algorithmus für LNM Parameterschätzung $\Rightarrow$ Lernen.


# Inferenz und Lernen
\footnotesize
\setstretch{1.1}
\underline{Beweis}

Wir zeigen zunächst, dass die ELBO $q^{(j)}(X) := p_{\theta^{(j-1)}}(X|Y)$ den maximalen Wert $\ln p_{\theta^{(j-1)}}(Y)$ annimmt:
\begin{align*}
\mbox{ELBO}(p_{\theta^{(j-1)}}(X|Y),\theta)
& = \int p_{\theta^{(j-1)}}(X|Y) \ln \left(\frac{p_{\theta^{(j-1)}}(X,Y)}{p_{\theta^{(j-1)}}(X|Y)}\right) \,dX 		\\
& = \int p_{\theta^{(j-1)}}(X|Y) \ln \left(\frac{p_{\theta^{(j-1)}}(Y)p_{\theta^{(j-1)}}(X|Y)}{p_{\theta^{(j-1)}}(X|Y)}\right) \,dX 	\\
& = \int p_{\theta^{(j-1)}}(X|Y) \ln p_{\theta^{(j-1)}}(Y) \,dX 										\\
& = \ln  p_{\theta^{(j-1)}}(Y) \int p_{\theta^{(j-1)}}(X|Y) \,dX \\
&  = \ln  p_{\theta^{(j-1)}}(Y).
\end{align*}
Der M Schritt hat dementsprechend die Form
\begin{align*}
\theta^{(j)}
& = \argmax_{\theta} \mbox{ELBO}\left(p_{\theta^{(j-1)}}(X|Y),\theta\right) 	\\
& = \argmax_{\theta} \int p_{\theta^{(j-1)}}(X|Y) \ln p_\theta(X,Y) \,dX - \int p_{\theta^{(j-1)}}(X|Y) \ln p_{\theta^{(j-1)}}(X|Y) \,dX.
\end{align*}
Der M Schritt des exakten EM Algorithmus folgt dann damit, dass der zweite Integralterm hier nicht von $\theta$ abhängt.


# Inferenz und Lernen
\vspace{2mm}
\textcolor{darkblue}{Visuelle Intuition}
\vspace{2mm}
```{r, echo = FALSE, out.width = "80%"}
knitr::include_graphics("5_Abbildungen/mvda_5_exakter_em_algorithmus.pdf")
```

# Inferenz und Lernen
Weitere Bemerkungen

\small
* \justifying Der M Step der $i$ten Iteration des exakten EM Algorithmus entspricht der
Maximierung des Erwartungswertes der logarithmierten WDF der gemeinsamen
Datenverteilung  $p_\theta(X,Y)$ hinsichtlich $\theta$, wobei der Erwartungswert
hinsichtlihch der WDF der bedingten Datenverteilung von $X$ gegeben $Y$
basierend auf der Parameterschätzung $\theta^{(j-1)}$, die in der $(j-1)$ten
Iteration des exakten EM Algorithmus gewonnen wurde.
\vspace{1mm}
* Überraschenderweise garantiert aufgrund der inherenten Logik des EM Algrithmus die
Maximierung des Erwartungswertes
\begin{equation}
\mathbb{E}_{p_{\theta^{(j-1)}}(X|Y)} \left(\ln p_\theta(X,Y)\right)
= \int p_{\theta^{(j-1)}}(X|Y) \ln p_\theta(X,Y) \,dX
\end{equation}
auch die Maximierung der  tatsächlichen Funktion von Interesse, $\ln p_\theta(Y)$.
* Für konkrete Algorithmen und für spezifische LNMe muss obiger Erwartungswert
analytisch als Funktion von $\theta^{(j-1)}$ ausgewertet werden und dann hinsichtlich $\theta$
entweder analytisch oder numerisch maximiert werden um die Parameterschätzung $\theta^{(j)}$.
zu erhalten.

#
\large
\setstretch{2}
\vfill
Grundlagen

Inferenz und Lernen

**Faktorenanalyse**

Hauptkomponentenanalyse Revisited

Anwendungsbeispiel

Selbstkontrollfragen

Appendix

# Faktorenanalyse
\small
\begin{definition}[Faktorenanalysemodell]
Es sei
\begin{equation}
p_\theta(x,y) = N(y; Bx, R)N(x; 0, I_k)
\end{equation}
ein LGM mit diagonaler Beobachtungsrauschen Kovarianzmatrix
\begin{equation}
R = \mbox{diag}(\sigma_1^2,\sigma_2^2,...,\sigma_m^2) \in \mathbb{R}^{m\times m}, \sigma_i^2 > 0,\,i = 1,...,m.
\end{equation}
Dann heißt $p_\theta(x,y)$ \textit{Faktorenanalysemodell}.
\end{definition}

Bemerkungen

* Das Modell heißt auch **Modell der konfirmatorischen Faktorenanalyse**.
* Die Komponenten (Zufallsvariablen) des latenten Zufallsvektors werden **Faktoren** genannt.
* Die Matrix $B$ des Faktorenanalysemodells wird **Faktorenladungsmatrix** genannt.
* Die Diagonalelemente von $R$ werden **Uniquenesses** genannt.

# Faktorenanalyse
\textcolor{darkblue}{Visuelle Intuition}
\vspace{4mm}

```{r, echo = FALSE, out.width = "100%"}
knitr::include_graphics("5_Abbildungen/mvda_5_faktorenanalysemodell.pdf")
```

# Faktorenanalyse
\setstretch{1.8}
\small
\justifying
Bemerkungen

\small
* Faktorenenalysemodelle erklären die Struktur Datenkovarianzmatrizen dadurch, dass
  * \footnotesize \small alle Korrelationen zwischen Datendimensionen durch $B$ erklärt werden,
  * alle Varianzen dimensionsspezifisch durch $R$ erklärt werden und
  * die Komponenten des beobachtbaren ZVs als bedingt unabhängig angenommen werden.
* Faktorenanalysemodelle behandeln Datenkovarianzen und Varianzen nicht identisch.
* Die Datenkovarianzstruktur wird als bedeutsam angesehen

  $\quad\quad\quad\quad\quad\quad\Leftrightarrow$ Das Beobachtungsrauschen wird als unkorreliert angenommen.
* Exploratorische und konfirmatorische FA entsprechen unterschiedlichen Bedingungen an $R$.
* Die Parameter des Modells können mit dem exakten EM Algorithmus geschätzt werden.


# Faktorenanalyse
\footnotesize
\setstretch{1.2}
\begin{theorem}[Exakter EM Algorithmus für Faktorenanalysemodelle]
\justifying
\normalfont
\vspace{1mm}

0. Initialisiere $B^{(0)}$ und $R^{(0)}$.
\vspace{1mm}

Für $i = 1,2,...$

\begin{enumerate}
\item E Schritt. Setze $\tilde{B} := B^{(j-1)}$ und $\tilde{R} := R^{(j-1)}$ und
\begin{equation}
\begin{split}
q^{(j)}(X)
:= \prod_{j=1}^n N\left(x^{(i)}; \hat{x}^{(i)}, \hat{\Sigma}^{(i)} \right),
\end{split}
\end{equation}
wobei
\begin{equation}
\hat{x}^{(i)} 		:= \tilde{B}^T(\tilde{B}\tilde{B}^T  + \tilde{R})^{-1}y^{(i)} \mbox{ und }
\hat{\Sigma}^{(i)}  := I_k - \tilde{B}^T(\tilde{B}\tilde{B}^T  + \tilde{R})^{-1}\tilde{B}.
\end{equation}
\item M Schritt. Setze
\begin{equation}
B^{(j)}
:= \sum_{i=1}^n y^{(i)}\hat{x}^{(i)^T}\left(\sum_{i=1}^n \hat{x}^{(i)}\hat{x}^{(i)^T} + \hat{\Sigma}^{(i)}\right)^{-1}
\end{equation}
und
\begin{equation}
R^{(j)}
:= \frac{1}{n}\mbox{diag}\left(\sum_{i=1}^n y^{(i)}y^{(i)^T} - \sum_{i=1}^n y^{(i)}\hat{x}^{(i)^T} B^{(i)^T}\right).
\end{equation}
\end{enumerate}
\end{theorem}

Siehe Appendix für einen Beweis.


# Faktorenanalyse
\textcolor{darkblue}{Simulationsbeispiel}

\small
Datengeneration

\vspace{1mm}
\footnotesize
\setstretch{.9}
```{r}
# mvrnorm() Paket
library(MASS)

# Paramterspezifikation
k     = 3                                # Dimension des latenten Zufallsvektors
m     = 9                                # Dimension des beobachtbaren Zufallsvektors
B     = matrix(c( 1,0,0,                 # Faktorenladungsmatrix
                  1,0,0,
                  1,0,0,
                  0,1,0,
                  0,1,0,
                  0,1,0,
                  0,0,1,
                  0,0,1,
                  0,0,1),
                 nrow  = m,
                 byrow = TRUE)
u      = rep(1,m)                        # Uniquenesses
R      = diag(u)                         # Beobachtungsrauschenkovarianmatrix
mu     = rep(0,k+m)                      # Erwartungswertparameter  p_\theta(x,y)
Sigma  = rbind(cbind(diag(k), t(B)),     # Kovarianzmatrixparameter p_\theta(x,y)
         cbind(B, B %*% t(B) + R))

# Datengeneration
n      = 1e3                             # Beobachtungsanzahl
XY     = t(mvrnorm(n,mu,Sigma))          # p_theta(x,y) sampling und z-score normalization
X      = XY[1:k,]                        # virtuelle Daten    X \in R^{k x n}
Y      = XY[(k+1):nrow(XY),]             # Beobachtete Daten  Y \in R^{m x n}
```


# Faktorenanalyse
\textcolor{darkblue}{Simulationsbeispiel}

\small
Wahre, aber unbekannte, Parameter

```{r, eval = F, echo = F}
library(plot.matrix)
library(RColorBrewer)
graphics.off()
fdir        =  file.path(getwd(), "5_Abbildungen")
dev.new()
par(
family      = "sans",
mfcol       = c(1,2),
pty         = c("s", "s"),
bty         = "l",
lwd         = 1,
las         = 1,
mgp         = c(2,1,0),
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = 1,
cex.main    = 3,
mar         = c(1,2,3,1)) # bottom, left, top, right margin
plot(
B,
breaks      = seq(-1.2,1.2,len = 11),
col         = rev(brewer.pal(n = 11, name = "RdBu")),
digits      = 2,
key         = NULL,
cex         = 0.8,
polygon.key = NULL,
axis.key    = NULL,
xlab        = "",
ylab        = "",
main        = TeX("$B$"))
plot(
R,
digits      = 2,
breaks      = seq(-1.2,1.2,len = 11),
col         = rev(brewer.pal(n = 11, name = "RdBu")),
key         = NULL,
cex         = 0.8,
polygon.key = NULL,
axis.key    = NULL,
xlab        = "",
ylab        = "",
main        = TeX("$R$"))
dev.copy2pdf(
file        = file.path(fdir, "mvda_5_simulation_B_R.pdf"),
width       = 12,
height      = 8)
```


```{r, echo = FALSE, out.width = "80%"}
knitr::include_graphics("5_Abbildungen/mvda_5_simulation_B_R.pdf")
```

# Faktorenanalyse
\textcolor{darkblue}{Simulationsbeispiel}

\small
Marginale Kovarianzmatrix und marginale Korrelationsmatrix

```{r, eval = F, echo = F}
# Marginale Kovarianzmatrix und marginale Korrelationsmatrix
library(plot.matrix)
library(RColorBrewer)
Sigma_y = B %*% t(B) + R
D_rho   = diag(1/sqrt(diag(Sigma_y)))
rho_y   = D_rho %*% Sigma_y %*% D_rho
graphics.off()
fdir        =  file.path(getwd(), "5_Abbildungen")
dev.new()
par(
family      = "sans",
mfcol       = c(1,2),
pty         = "s",
bty         = "l",
lwd         = 1,
las         = 1,
mgp         = c(3,1,0),
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = 1,
cex.main    = 3,
mar         = c(2,2,3,2)) # bottom, left, top, right margin
plot(
Sigma_y,
breaks      = seq(-2.3,2.3,len = 11),
col         = rev(brewer.pal(n = 11, name = "RdBu")),
digits      = 1,
key         = NULL,
cex         = 0.8,
polygon.key = NULL,
axis.key    = NULL,
xlab        = "",
ylab        = "",
main        = TeX("$\\Sigma_y$"))
plot(
rho_y,
digits      = 1,
breaks      = seq(-1.2,1.2,len = 11),
col         = rev(brewer.pal(n = 11, name = "RdBu")),
key         = NULL,
cex         = 0.8,
polygon.key = NULL,
axis.key    = NULL,
xlab        = "",
ylab        = "",
main        = TeX("$\\rho_y$"))
dev.copy2pdf(
file        = file.path(fdir, "mvda_5_simulation_Sigma_y_rho_y.pdf"),
width       = 12,
height      =  8)
```

```{r, echo = FALSE, out.width = "80%"}
knitr::include_graphics("5_Abbildungen/mvda_5_simulation_Sigma_y_rho_y.pdf")
```

# Faktorenanalyse
\textcolor{darkblue}{Simulationsbeispiel}

\small
Vollständiger Datensatz (n = 100)

```{r, eval = F, echo = F}
# Vollständiger Datensatz (n = 100)
library(plot.matrix)
library(RColorBrewer)
graphics.off()
fdir        =  file.path(getwd(), "5_Abbildungen")
dev.new()
par(
family       = "sans",
mfcol        = c(1,1),
pty          = "m",
bty          = "l",
lwd          = 1,
las          = 1,
mgp          = c(3,1,0),
xaxs         = "i",
yaxs         = "i",
font.main    = 1,
cex          = 1.5,
cex.main     = 1.5,
mar          = c(2,4,3,4)) # bottom, left, top, right margin
varnames     = c("x1",
                 "x2",
                 "x3",
                 "y1",
                 "y2",
                 "y3",
                 "y4",
                 "y5",
                 "y6",
                 "y7",
                 "y8",
                 "y9")
rownames(XY) = varnames
plot(
XY,
border       = NA,
breaks       = seq(-6,6,len = 11),
col          = rev(brewer.pal(n = 11, name = "RdBu")),
key          = list(side     = 4,
                    font     = 1,
                    cex.axis = 1),
fmt.key      = "%.2f",
polygon.key  = NULL,
axis.key     = NULL,
spacing.key  = c(3,2,2),
cex          = 0.8,
xlab         = "Realisierungen",
ylab         = "",
main         = TeX("$(X,Y) \\in R^{(k + m) \\times n}$"))
dev.copy2pdf(
file         = file.path(fdir, "mvda_5_simulation_XY.pdf"),
width        = 12,
height       =  8)
```

```{r, echo = FALSE, out.width = "80%"}
knitr::include_graphics("5_Abbildungen/mvda_5_simulation_XY.pdf")
```


# Faktorenanalyse
\textcolor{darkblue}{Simulationsbeispiel}

\small
Stichprobenkovarianzmatrix und Stichprobenkorrelationsmatrix

```{r, eval = F, echo = F}

# Stichprobenkovarianzmatrix und Stichprobenkorrelationsmatrix
library(plot.matrix)
library(RColorBrewer)
I_n         = diag(n)                                     # Einheitsmatrix I_n
J_n         = matrix(rep(1,n^2), nrow = n)                # 1_{nn}
C_y         = (1/(n-1))*(Y %*% (I_n-(1/n)*J_n) %*% t(Y))  # Stichprobenkovarianzmatrix
D_y         = diag(1/sqrt(diag(C_y)))                     # Kov-Korr-Transformationsmatrix
R_y         = D_y %*% C_y %*% D_y                         # Stichprobenkorrelationsmatrix

graphics.off()
fdir        =  file.path(getwd(), "5_Abbildungen")
dev.new()
par(
family      = "sans",
mfcol       = c(1,2),
pty         = "s",
bty         = "l",
lwd         = 1,
las         = 1,
mgp         = c(3,1,0),
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = 1,
cex.main    = 3,
mar         = c(2,2,3,2)) # bottom, left, top, right margin
plot(
C_y,
breaks      = seq(-3,3,len = 11),
col         = rev(brewer.pal(n = 11, name = "RdBu")),
digits      = 1,
key         = NULL,
cex         = 0.8,
polygon.key = NULL,
axis.key    = NULL,
xlab        = "",
ylab        = "",
main        = TeX("$\\S_y$"))
plot(
R_y,
digits      = 1,
breaks      = seq(-1.2,1.2,len = 11),
col         = rev(brewer.pal(n = 11, name = "RdBu")),
key         = NULL,
cex         = 0.8,
polygon.key = NULL,
axis.key    = NULL,
xlab        = "",
ylab        = "",
main        = TeX("$\\R_y$"))
dev.copy2pdf(
file        = file.path(fdir, "mvda_5_simulation_S_y_R_y.pdf"),
width       = 12,
height      =  8)
```

```{r, echo = FALSE, out.width = "80%"}
knitr::include_graphics("5_Abbildungen/mvda_5_simulation_S_y_R_y.pdf")
```


# Faktorenanalyse

\textcolor{darkblue}{Simulationsbeispiel}

\small
Parameterschätzung mit R's `factanal()` Funktion

\footnotesize

```{r}
# Lawley & Maxwell (1971) Schätzverfahren für Faktorenanalyse
fa              = factanal(t(Y), factors = 3, rotation = "none")  # Parameterschätzung
B_hat_fa        = fa$loadings[1:9,1:3]                            # \hat{B}
R_hat_fa        = diag(fa$uniquenesses)                           # \hat{R}

# Parametschätzer-basierte Kovarianz und Korrelationsmatrix
Sigma_hat_y_fa  = B_hat_fa %*% t(B_hat_fa) + R_hat_fa
D_hat_y_fa      = diag(1/sqrt(diag(Sigma_hat_y_fa)))
rho_hat_y_fa    = D_hat_y_fa %*% Sigma_hat_y_fa %*% D_hat_y_fa
```



# Faktorenanalyse
\textcolor{darkblue}{Simulationsbeispiel}

\small
`factanal()` basierte Parameterschätzer

```{r, eval = F, echo = F}
library(plot.matrix)
library(RColorBrewer)
graphics.off()
fdir        =  file.path(getwd(), "5_Abbildungen")
dev.new()
par(
family      = "sans",
mfcol       = c(1,2),
pty         = c("s", "s"),
bty         = "l",
lwd         = 1,
las         = 1,
mgp         = c(2,1,0),
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = 1,
cex.main    = 3,
mar         = c(1,2,4,1)) # bottom, left, top, right margin
plot(
B_hat_fa,
breaks      = seq(-1.2,1.2,len = 11),
col         = rev(brewer.pal(n = 11, name = "RdBu")),
digits      = 2,
key         = NULL,
cex         = 0.8,
polygon.key = NULL,
axis.key    = NULL,
xlab        = "",
ylab        = "",
main        = TeX("$\\hat{B}_{FA}$"))
plot(
R_hat_fa,
digits      = 2,
breaks      = seq(-1.2,1.2,len = 11),
col         = rev(brewer.pal(n = 11, name = "RdBu")),
key         = NULL,
cex         = 0.8,
polygon.key = NULL,
axis.key    = NULL,
xlab        = "",
ylab        = "",
main        = TeX("$\\hat{R}_{FA}$"))
dev.copy2pdf(
file        = file.path(fdir, "mvda_5_simulation_B_hat_R_hat_fa.pdf"),
width       = 12,
height      = 8)
```


```{r, echo = FALSE, out.width = "80%"}
knitr::include_graphics("5_Abbildungen/mvda_5_simulation_B_hat_R_hat_fa.pdf")
```


# Faktorenanalyse
\textcolor{darkblue}{Simulationsbeispiel}

\small
`factanal()` basierte marginale Kovarianmatrix und Korrelationsmatrix

```{r, eval = F, echo = F}

# factanal() basierte marginale Kovarianmatrix und Korrelationsmatrix
library(plot.matrix)
graphics.off()
fdir        =  file.path(getwd(), "5_Abbildungen")
dev.new()
par(
family      = "sans",
mfcol       = c(1,2),
pty         = "s",
bty         = "l",
lwd         = 1,
las         = 1,
mgp         = c(3,1,0),
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = 1,
cex.main    = 2,
mar         = c(2,2,3,2)) # bottom, left, top, right margin
plot(
Sigma_hat_y_fa,
breaks      = seq(-1.2,1.2,len = 11),
col         = rev(brewer.pal(n = 11, name = "RdBu")),
digits      = 1,
key         = NULL,
cex         = 0.8,
polygon.key = NULL,
axis.key    = NULL,
xlab        = "",
ylab        = "",
main        = TeX("$\\hat{\\Sigma}_y$"))
plot(
rho_hat_y_fa,
digits      = 1,
breaks      = seq(-1.2,1.2,len = 11),
col         = rev(brewer.pal(n = 11, name = "RdBu")),
key         = NULL,
cex         = 0.8,
polygon.key = NULL,
axis.key    = NULL,
xlab        = "",
ylab        = "",
main        = TeX("$\\hat{\\rho}_y$"))
dev.copy2pdf(
file        = file.path(fdir, "mvda_5_simulation_Sigma_hat_y_rho_hat_y_fa.pdf"),
width       = 12,
height      =  8)
```
```{r, echo = FALSE, out.width = "80%"}
knitr::include_graphics("5_Abbildungen/mvda_5_simulation_Sigma_hat_y_rho_hat_y_fa.pdf")
```

# Faktorenanalyse
\setstretch{1.2}
\textcolor{darkblue}{Simulationsbeispiel}

\small
EM Algorithmus

\tiny
\setstretch{.9}
```{r, eval = F}
library(matlib)                                   # Matrizenrechnungspaket
max_j = 2^2                                       # Maximale Anzahl Iterationen

# EM Initialisierung
B_j = matrix(runif(m*k), nrow = m)
R_j = diag(runif(0,m))

# Iterations
for(j in 1:max_j){

  # E Schritt
  X_hat_j      = t(B_j) %*% inv(B_j %*% t(B_j) + R) %*% Y
  Sigma_hat_j  = diag(k) - (t(B_j) %*% inv(B_j %*% t(B_j) + R) %*% B_j)

  # M Schritt
  yxT = Y %*% t(X_hat_j)
  xxT = X_hat_j %*% t(X_hat_j)
  yyT = Y %*% t(Y)
  B_j = yxT %*% inv(xxT + Sigma_hat_j)
  R_j = (1/n) * diag(diag((yyT - (yxT %*% t(B_j)))))
}
print(yyT)
print(Y %*% t(Y))
# Parameterschätzer
B_hat_em = B_j
R_hat_em = R_j

# Parametschätzer-basierte Kovarianz und Korrelationsmatrix
Sigma_hat_y_em  = B_hat_em %*% t(B_hat_em) + R_hat_em
D_hat_y_em      = diag(1/sqrt(diag(Sigma_hat_y_em)))
rho_hat_y_em    = D_hat_y_em %*% Sigma_hat_y_em %*% D_hat_y_em

```

# Faktorenanalyse
\textcolor{darkblue}{Simulationsbeispiel}

\small
EM Algorithmus-basierte Parameterschätzer

```{r, eval = F, echo = F}
library(plot.matrix)
library(RColorBrewer)
graphics.off()
fdir        =  file.path(getwd(), "5_Abbildungen")
dev.new()
par(
family      = "sans",
mfcol       = c(1,2),
pty         = c("s", "s"),
bty         = "l",
lwd         = 1,
las         = 1,
mgp         = c(2,2,0),
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = 1,
cex.main    = 2.5,
mar         = c(1,2,4,1)) # bottom, left, top, right margin
plot(
B_hat_em,
breaks      = seq(-3.2,3.2,len = 11),
col         = rev(brewer.pal(n = 11, name = "RdBu")),
digits      = 2,
key         = NULL,
cex         = 0.8,
polygon.key = NULL,
axis.key    = NULL,
xlab        = "",
ylab        = "",
main        = TeX("$\\hat{B}_{EM}$"))
plot(
R_hat_em,
digits      = 2,
breaks      = seq(-1.2,1.2,len = 11),
col         = rev(brewer.pal(n = 11, name = "RdBu")),
key         = NULL,
cex         = 0.8,
polygon.key = NULL,
axis.key    = NULL,
xlab        = "",
ylab        = "",
main        = TeX("$\\hat{R}_{EM}$"))
dev.copy2pdf(
file        = file.path(fdir, "mvda_5_simulation_B_hat_R_hat_em.pdf"),
width       = 12,
height      = 8)
```


```{r, echo = FALSE, out.width = "80%"}
knitr::include_graphics("5_Abbildungen/mvda_5_simulation_B_hat_R_hat_em.pdf")
```

# Faktorenanalyse
\textcolor{darkblue}{Simulationsbeispiel}

\small
EM Algorithmus basierte marginale Kovarianmatrix und Korrelationsmatrix

```{r, eval = F, echo = F}
# EM Algorithmus basierte marginale Kovarianmatrix und Korrelationsmatrix
library(plot.matrix)
graphics.off()
fdir        =  file.path(getwd(), "5_Abbildungen")
dev.new()
par(
family      = "sans",
mfcol       = c(1,2),
pty         = "s",
bty         = "l",
lwd         = 1,
las         = 1,
mgp         = c(3,1,0),
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = 1,
cex.main    = 2,
mar         = c(2,2,3,2)) # bottom, left, top, right margin
plot(
Sigma_hat_y_em,
breaks      = seq(-20,20,len = 11),
col         = rev(brewer.pal(n = 11, name = "RdBu")),
digits      = 1,
key         = NULL,
cex         = 0.8,
polygon.key = NULL,
axis.key    = NULL,
xlab        = "",
ylab        = "",
main        = TeX("$\\hat{\\Sigma}_y^{EM}$"))
plot(
rho_hat_y_em,
digits      = 1,
breaks      = seq(-1.2,1.2,len = 11),
col         = rev(brewer.pal(n = 11, name = "RdBu")),
key         = NULL,
cex         = 0.8,
polygon.key = NULL,
axis.key    = NULL,
xlab        = "",
ylab        = "",
main        = TeX("$\\hat{\\rho}_y^{EM}$"))
dev.copy2pdf(
file        = file.path(fdir, "mvda_5_simulation_Sigma_hat_y_rho_hat_y_em.pdf"),
width       = 12,
height      =  8)
```

```{r, echo = FALSE, out.width = "80%"}
knitr::include_graphics("5_Abbildungen/mvda_5_simulation_Sigma_hat_y_rho_hat_y_em.pdf")
```


#
\large
\setstretch{2}
\vfill
Grundlagen

Inferenz und Lernen

Faktorenanalyse

**Hauptkomponentenanalyse Revisited**

Selbstkontrollfragen

Appendix


# Hauptkomponentenanalyse Revisited
\small
\begin{definition}[Modell der probabilistischen Hauptkomponentenanalyse]
Es sei
\begin{equation}
p_\theta(x,y) = N(y; Bx, R)N(x; 0_k, I_k)
\end{equation}
ein LNM mit der sphärischen Beobachtungsrauschen Kovarianzmatrix
\begin{equation}
R := \sigma^2 I_m.
\end{equation}
Dann heißt $p_\theta(x,y)$ \textit{probabilistisches Hauptkomponentanalysemodell}.
\end{definition}

\footnotesize
Bemerkungen

* Die Matrix $B$ bedingt die Beziehung zur klassischen Hauptkomponentenanalyse.
* $\sigma^2$ heißt *globales Rauschlevel*.
* $B$ and $\sigma^2$ können durch den EM Algorithmus geschätzt werden.
* $B$ and $\sigma^2$ können auch durch direkte Maximierung der marginalen Likelihood Funktion geschätzt werden.

# Hauptkomponentenanalyse Revisited
\textcolor{darkblue}{Visuelle Intuition zur probabilistischen Hauptkomponentenanalyse}
\vspace{6mm}

```{r, echo = FALSE, out.width = "100%"}
knitr::include_graphics("5_Abbildungen/mvda_5_ppcamodell.pdf")
```


# Hauptkomponentenanalyse Revisited
\small
\begin{theorem}[Parameter der probabilistischen Hauptkomponentnanalyse]
\justifying
\normalfont
Es sei
\begin{equation}
p_\theta(x,y) = N(y; Bx, \sigma^2 I_m)N(x; 0, I_k)
\end{equation}
ein probabilitisches Hauptkomponentenanalysemodell und es sei
\begin{equation}
\mathbb{C}(Y) = Q\Lambda Q^T
\end{equation}
die Orthogonalzerlegung der Kovarianzmatrix des beobachtbaren Zufallsvektors.
Dann kann der Parameter $B$ des probabilistischen Hauptkomponentnanalysemodells
geschrieben werden als
\begin{equation}
B = Q\left(\Lambda - \sigma^2I_m\right)^{1/2}.
\end{equation}
\end{theorem}

\footnotesize
Bemerkungen

* Wir haben $\mathbb{C}(Y) = Q\Lambda Q^T$ Hauptkomponentenanalyse von $\mathbb{C}(Y)$ genannt.
* $Q$ besteht aus den Eigenvektoren von $\mathbb{C}(y)$, $\Lambda$ aus den assoziierten Eigenwerten.
* Wir haben die Eigenvektoren von $\mathbb{C}(Y)$ die Hauptkomponenten von $\mathbb{C}(Y)$ genannt.
* Die Spalten von $B$ sind die Hauptkomponenten gewichtet mit $(\Lambda - \sigma^2I_m)^{1/2}$.


# Hauptkomponentenanalyse Revisited
\footnotesize
\underline{Beweis}

Wir halten zunächst fest, dass die marginale Datenverteilung eines probabilistischen
Hauptkomponentenanalysemodells gegeben ist durch
\begin{equation}
p_\theta(y) = N(y;0_m, BB^T + \sigma^2 I_m)
\mbox{ and thus }
\mathbb{C}(Y) =  BB^T + \sigma^2 I_m.
\end{equation}
Einsetzen von $B = Q(\Lambda - \sigma^2I_m)^{1/2}$ ergibt dann
\begin{align}
\begin{split}
\mathbb{C}(Y)
& =
BB^T
+ \sigma^2 I_m
\\
& = Q(\Lambda - \sigma^2I_m)^{1/2}(Q(\Lambda - \sigma^2I_m)^{1/2})^T + \sigma^2 I_m
\\
& = Q(\Lambda - \sigma^2I_m)^{1/2}((\Lambda - \sigma^2I_m)^{1/2})^TQ^T + \sigma^2 I_m \\
& = Q(\Lambda - \sigma^2I_m)^{1/2}(\Lambda - \sigma^2I_m)^{1/2}Q^T + \sigma^2 I_m \\
& = Q(\Lambda - \sigma^2I_m)Q^T + \sigma^2 I_m \\
& = (Q\Lambda - \sigma^2QI_m)Q^T + \sigma^2 I_m \\
& = Q\Lambda Q^T - \sigma^2QQ^T + \sigma^2 I_m \\
& = Q\Lambda Q^T - \sigma^2 I_m + \sigma^2 I_m \\
& = Q\Lambda Q^T.
\end{split}
\end{align}
Also ergibt sich die Äquivalenz
\begin{equation}
\mathbb{C}(y) = Q\Lambda Q^T
\Leftrightarrow
B = Q(\Lambda - \sigma^2I_m)^{\frac{1}{2}}.
\end{equation}



# Hauptkomponentenanalyse Revisited
\footnotesize
\begin{theorem}[Direkte Marginal-Maximum Likelihood Schätzung]
\normalfont
\justifying
Es sei
\begin{equation}
p_\theta(x,y) = N(y;Bx,\sigma^2 I_m)N(x;0,I_k)
\end{equation}
ein probabilistisches Hauptkomponentenanalyse model und $Y \in \mathbb{R}^{m \times n}$
sei ein von diesem Modell generierter Datensatz. Weiterhin sei
\begin{equation}
C = \frac{1}{n}YY^T
\mbox{ und }
C = Q \Lambda Q^T
\end{equation}
die unkorrigierte Stichprobenkovarianzmatrix und ihre Orhtogonazerlegung, respektive.
Schließlich seien $Q_q \in \mathbb{R}^{m \times q}$ und $\Lambda_q \in \mathbb{R}^{q \times q}$
die Matrizen die durch die spaltenweise Konkatenation der höchsten Eigenwerte und
ihrer assoziierten Eigenvektoren entstehen. Dann sind Maximum Marginal Likelihood
Schätzer von $B$ und $\sigma^{2}$ durch
\begin{equation}
\hat{B} = Q_q(\Lambda_q - \sigma^2 I_m)^{1/2}
\mbox{ and }
\hat{\sigma}^2 = \frac{1}{m - l} \sum_{j = l + 1}^m \lambda_j,
\end{equation}
gegeben.
\end{theorem}
Siehe Appendix für einen Beweis.


# Hauptkomponentenanalyse Revisited
\small
\begin{definition}[Modell der Hauptkomponentenanalyse]
Es sei
\begin{equation}
p_\theta(x,y) = N(y; Bx, R)N(x; 0, I_k)
\end{equation}
ein LNM mit der Beobachtungsrauschenkovarianzmatrix
\begin{equation}
R := \lim_{\sigma^2  \to 0} \sigma^2  I_m \in \mathbb{R}^{m\times m}.
\end{equation}
Dann heißt $p_\theta(x,y)$ \textit{Modell der Hauptkomponentenanalyse}.
\end{definition}
Bemerkungen
\begin{itemize}
\item $B$ enkodiert die Beziehung zur klassischen Hauptkomponentenanalyse.
\item Das Beobachtungsrauschen wird als nicht-existent angenommen.
\end{itemize}


# Hauptkomponentenanalyse Revisited
\textcolor{darkblue}{Visuelle Intuition zur Hauptkomponentenanalyse}
\vspace{6mm}

```{r, echo = FALSE, out.width = "100%"}
knitr::include_graphics("5_Abbildungen/mvda_5_pcamodell.pdf")
```

# Hauptkomponentenanalyse Revisited
\small
\begin{theorem}[Parameter der Hauptkomponentenanalyse]
\justifying
\normalfont
Es sei
\begin{equation}
p_\theta(x,y) = N\left(y; Bx, \lim_{\sigma^2  \to 0} \sigma^2 I_m\right)N(x; 0_k, I_k)
\end{equation}
ein Hauptkomponentenanalysemodell und es sei
\begin{equation}
\mathbb{C}(Y) = Q^T \Lambda Q
\end{equation}
die Hauptkomponentenanalyse der zugehörigen marginalen Datenverteilung. Dann gilt
\begin{equation}
B = Q\Lambda^{\frac{1}{2}}.
\end{equation}
\end{theorem}

\footnotesize
\underline{Beweis}

Wir halten zunächst fest, dass die marginale Datenverteilung des
Hauptkomponentenanalysemodells im Limit $\sigma^2 \to 0$ gegeben ist durch
\begin{equation}
p_\theta(y) = N(y;0_m, BB^T)
\mbox{ and thus }
\mathbb{C}(y) = BB^T.
\end{equation}
Es ergibt sich also
\begin{equation}
\mathbb{C}(y) = BB^T \Leftrightarrow Q\Lambda Q^T = BB^T \Leftrightarrow B = Q \Lambda^{\frac{1}{2}}.
\end{equation}
$\hfill \Box$


#
\large
\setstretch{2}
\vfill
Grundlagen

Inferenz und Lernen

Faktorenanalyse

Hauptkomponentenanalyse Revisited

**Anwendungsbeispiel**

Selbstkontrollfragen

Appendix
\vfill


# Anwendungsbeispiel
\textcolor{darkblue}{Nine mental ability tests | Intelligenzforschungsdatensatz} 

@holzinger_1939

\small
Visualization

1. Visual Perception
2. Cubes
3. Lozenges

Verbal intelligence

4. Paragraph Comprehension
5. Sentence Completion
6. Word Meaning

Speed

7. Addition
8. Counting dots
9. Straight-Curved Capitals


# Anwendungsbeispiel
\textcolor{darkblue}{Beobachteter Datensatz (n = 301)}
\small

301 Proband:innen | 11 - 16 Jahre

Probandin:innen 1 - 10

\footnotesize
```{r, echo = F, message = F}
# Datensatz
library(lavaan)                                      # Lavaan SEM Paket
data(HolzingerSwineford1939)                         # Datensatz
Y             = t(HolzingerSwineford1939[,7:15])     # Datenmatrix
m             = nrow(Y)                              # Anzahl Tests/Variablen
n             = ncol(Y)                              # Anzahl Proband:innen
tests         = c("Perception",                      # Variablennamen
                  "Cubes",
                  "Lozenges",
                  "Comprehension",
                  "Completion",
                  "Word Meaning",
                  "Addition",
                  "Counting",
                  "Capitals")
rownames(Y)   = tests
colnames(Y)   = 1:ncol(Y)

# Tabellevisualisierung
knitr::kable(Y[,1:10], digits = 2, "pipe")
```
# Anwendungsbeispiel

\textcolor{darkblue}{Beobachteter Datensatz (n = 301)}

```{r, eval = F, echo = F}
# Beobachteter Datensatz (n = 301)
library(plot.matrix)
library(RColorBrewer)
graphics.off()
fdir        =  file.path(getwd(), "5_Abbildungen")
dev.new()
par(
family       = "sans",
mfcol        = c(1,1),
pty          = "m",
bty          = "l",
lwd          = 1,
las          = 1,
mgp          = c(3,1,0),
xaxs         = "i",
yaxs         = "i",
font.main    = 1,
cex          = 1.2,
cex.main     = 1.7,
mar          = c(2,7,3,4)) # bottom, left, top, right margin
plot(
Y,
border       = NA,
breaks       = seq(0,9,len = 11),
col          = rev(brewer.pal(n = 11, name = "RdBu")),
key          = list(side     = 4,
                    font     = 1,
                    cex.axis = 1),
fmt.key      = "%.2f",
polygon.key  = NULL,
axis.key     = NULL,
spacing.key  = c(3,2,2),
cex          = 0.8,
xlab         = "Realisierungen",
ylab         = "",
main         = TeX("$(X,Y) \\in R^{(k + m) \\times n}$"))
dev.copy2pdf(
file         = file.path(fdir, "mvda_5_anwendung_Y.pdf"),
width        = 12,
height       =  8)
```

```{r, echo = FALSE, out.width = "80%"}
knitr::include_graphics("5_Abbildungen/mvda_5_anwendung_Y.pdf")
```

# Anwendungsbeispiel

\textcolor{darkblue}{Stichprobenkovarianzmatrix und Stichprobenkorrelationsmatrix}

```{r, eval = F, echo = F}

# Stichprobenkovarianzmatrix und Stichprobenkorrelationsmatrix
library(plot.matrix)
library(RColorBrewer)
I_n         = diag(n)                                     # Einheitsmatrix I_n
J_n         = matrix(rep(1,n^2), nrow = n)                # 1_{nn}
C_y         = (1/(n-1))*(Y %*% (I_n-(1/n)*J_n) %*% t(Y))  # Stichprobenkovarianzmatrix
D_y         = diag(1/sqrt(diag(C_y)))                     # Kov-Korr-Transformationsmatrix
R_y         = D_y %*% C_y %*% D_y                         # Stichprobenkorrelationsmatrix

graphics.off()
fdir        =  file.path(getwd(), "5_Abbildungen")
dev.new()
par(
family      = "sans",
mfcol       = c(1,2),
pty         = "s",
bty         = "l",
lwd         = 1,
las         = 2,
mgp         = c(3,1,0),
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = 1,
cex.main    = 2,
mar         = c(2,8,3,2)) # bottom, left, top, right margin
plot(
C_y,
breaks      = seq(-2.1,2.1,len = 11),
col         = rev(brewer.pal(n = 11, name = "RdBu")),
digits      = 1,
key         = NULL,
cex         = 0.8,
polygon.key = NULL,
axis.key    = NULL,
xlab        = "",
ylab        = "",
main        = TeX("$\\S_y$"))
rownames(R_y) = row.names(C_y)
colnames(R_y) = row.names(C_y)
plot(
R_y,
digits      = 1,
breaks      = seq(-1.2,1.2,len = 11),
col         = rev(brewer.pal(n = 11, name = "RdBu")),
key         = NULL,
cex         = 0.8,
polygon.key = NULL,
axis.key    = NULL,
xlab        = "",
ylab        = "",
main        = TeX("$\\R_y$"))
dev.copy2pdf(
file        = file.path(fdir, "mvda_5_anwendung_S_y_R_y.pdf"),
width       = 12,
height      =  8)
```

```{r, echo = FALSE, out.width = "90%"}
knitr::include_graphics("5_Abbildungen/mvda_5_anwendung_S_y_R_y.pdf")
```


# Anwendungsbeispiel

\textcolor{darkblue}{Modellschätzung}
\small

Die Art der Tests motiviert ein 3-Faktor-Modell mit Faktoren

* Visualisierungsvermögen
* Verbale Intelligenz
* Schnelligkeit

\vspace{2mm}
\footnotesize
```{r}
# Lawley & Maxwell (1971) Schätzverfahren für Faktorenanalyse
fa              = factanal(t(Y), factors = 3, rotation = "none")  # Parameterschätzung
B_hat_fa        = fa$loadings[1:9,1:3]                            # \hat{B}
R_hat_fa        = diag(fa$uniquenesses)                           # \hat{R}

# Parametschätzer-basierte Kovarianz und Korrelationsmatrix
Sigma_hat_y_fa  = B_hat_fa %*% t(B_hat_fa) + R_hat_fa
D_hat_y_fa      = diag(1/sqrt(diag(Sigma_hat_y_fa)))
rho_hat_y_fa    = D_hat_y_fa %*% Sigma_hat_y_fa %*% D_hat_y_fa
```

# Anwendungsbeispiel
\textcolor{darkblue}{Modellschätzung}

```{r, eval = F, echo = F}
library(plot.matrix)
library(RColorBrewer)
graphics.off()
fdir        =  file.path(getwd(), "5_Abbildungen")
dev.new()
par(
family      = "sans",
mfcol       = c(1,2),
pty         = c("s", "s"),
bty         = "l",
lwd         = 1,
las         = 1,
mgp         = c(2,1,0),
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = 1,
cex.main    = 2,
mar         = c(1,2,4,1)) # bottom, left, top, right margin
plot(
B_hat_fa,
breaks      = seq(-1.2,1.2,len = 11),
col         = rev(brewer.pal(n = 11, name = "RdBu")),
digits      = 2,
key         = NULL,
cex         = 0.8,
polygon.key = NULL,
axis.key    = NULL,
xlab        = "",
ylab        = "",
main        = TeX("$\\hat{B}_{FA}$"))
plot(
R_hat_fa,
digits      = 2,
breaks      = seq(-1.2,1.2,len = 11),
col         = rev(brewer.pal(n = 11, name = "RdBu")),
key         = NULL,
cex         = 0.8,
polygon.key = NULL,
axis.key    = NULL,
xlab        = "",
ylab        = "",
main        = TeX("$\\hat{R}_{FA}$"))
dev.copy2pdf(
file        = file.path(fdir, "mvda_5_anwendung_B_hat_R_hat_fa.pdf"),
width       = 12,
height      = 8)
```


```{r, echo = FALSE, out.width = "80%"}
knitr::include_graphics("5_Abbildungen/mvda_5_anwendung_B_hat_R_hat_fa.pdf")
```


# Anwendungsbeispiel
\textcolor{darkblue}{Marginale Kovarianmatrix und Korrelationsmatrix}


```{r, eval = F, echo = F}

# factanal() basierte marginale Kovarianmatrix und Korrelationsmatrix
rownames(Sigma_hat_y_fa) = row.names(C_y)
colnames(Sigma_hat_y_fa) = row.names(C_y)
rownames(rho_hat_y_fa)   = row.names(C_y)
colnames(rho_hat_y_fa)   = row.names(C_y)

library(plot.matrix)
graphics.off()
fdir        =  file.path(getwd(), "5_Abbildungen")
dev.new()
par(
family      = "sans",
mfcol       = c(1,2),
pty         = "s",
bty         = "l",
lwd         = 1,
las         = 2,
mgp         = c(3,2,0),
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = 1,
cex.main    = 2,
mar         = c(2,8,3,2)) # bottom, left, top, right margin
plot(
Sigma_hat_y_fa,
breaks      = seq(-1.2,1.2,len = 11),
col         = rev(brewer.pal(n = 11, name = "RdBu")),
digits      = 1,
key         = NULL,
cex         = 0.8,
polygon.key = NULL,
axis.key    = NULL,
xlab        = "",
ylab        = "",
main        = TeX("$\\hat{\\Sigma}_y$"))
plot(
rho_hat_y_fa,
digits      = 1,
breaks      = seq(-1.2,1.2,len = 11),
col         = rev(brewer.pal(n = 11, name = "RdBu")),
key         = NULL,
cex         = 0.8,
polygon.key = NULL,
axis.key    = NULL,
xlab        = "",
ylab        = "",
main        = TeX("$\\hat{\\rho}_y$"))
dev.copy2pdf(
file        = file.path(fdir, "mvda_5_anwendung_Sigma_hat_y_rho_hat_y_fa.pdf"),
width       = 12,
height      =  8)
```

```{r, echo = FALSE, out.width = "90%"}
knitr::include_graphics("5_Abbildungen/mvda_5_anwendung_Sigma_hat_y_rho_hat_y_fa.pdf")
```

\normalsize
\begin{center}
$\Rightarrow$ Modellierung des Datensatzes mit einem 3-Faktor Modell ist möglich.
\end{center}

# Anwendungsbeispiel

\textcolor{darkblue}{Modellvergleich | Welche Anzahl von Faktoren zwischen 1 und 5 ist am besten?}

\small
Bayesian Information Criterion (cf. @horvath_2021, @schwarz_1978)
\normalsize
\begin{equation}
\mbox{BIC} := \ln p_{\hat{\theta}}(Y) - \frac{l}{2} \ln n
\end{equation}

$\ln p_{\hat{\theta}}(Y)$  

\small
* Logarithmierte marginale Datenwahrscheinlichkeit(sdichte) bei optimierten Parametern
* Maß für die Passungsgüte des Modells

\normalsize
$\frac{l}{2} \ln n$

\small
* Stichprobengröße gewichtete Anzahl an Parametern $l$ 
* Maß für die Komplexität des Modells

\normalsize
\begin{equation}
\mbox{BIC} = \mbox{Passungsgüte} - \mbox{Komplexität}
\end{equation}

# Anwendungsbeispiel

\textcolor{darkblue}{Modellvergleich | Welche Anzahl von Faktoren zwischen 1 und 5 ist am besten?}

\vfill
\small
\begin{theorem}[BIC für Faktorenanalysemodelle]
\normalfont
\justifying
$p_\theta(X,Y)$ sei die WDF der gemeinsamen Verteilung eines Faktorenanalysemodelldatensatzes
mit $X \in \mathbb{R}^{k \times n}$ und $Y \in \mathbb{R}^{m \times n}$.  Weiterhin 
seien $\hat{\theta} := \{\hat{B},\hat{R}\}$ (marginale) Maximum Likelihood Schätzer 
von $\theta := \{B,R\}$ und es sei
\begin{equation}
\hat{\Sigma} := \hat{B}\hat{B}^T + \hat{R}.
\end{equation}
Dann ergibt sich das Bayesian Information Criterion zu  
\begin{equation}
\mbox{BIC} =  -\frac{nm}{2}\ln 2\pi - \frac{n}{2} \ln |\hat{\Sigma}| - \frac{1}{2} \mbox{tr}\left(\hat{\Sigma}^{-1}YY^T\right)  - \frac{mk + m}{2}\ln n
\end{equation}
\end{theorem}
\vfill

# Anwendungsbeispiel

\textcolor{darkblue}{Modellvergleich | Welche Anzahl von Faktoren zwischen 1 und 5 ist am besten?}
\footnotesize

\setstretch{1}
\underline{Beweis}

Nach Definition der Verteilung von LNM Datensätzen gilt
\begin{equation}
p_{\hat{\theta}}(Y) = \prod_{i=1}^n N\left(y^{(i)}; 0_m, \hat{B}\hat{B}^T + \hat{R} \right).
\end{equation}
Die Anzahl der Parameter eines Faktorenanalysemodells ergibt sich als $l = mk + m$, 
wobei $mk$ die Anzahl der Einträge in $B$ und $m$ die Anzahl der Einträge in $R$ sind.
Mit der Eigenschaft 
\begin{equation}
x^T Ax = \mbox{tr}(Axx^T) 
\end{equation}
der Matrix Spur
\begin{equation}
\mbox{tr}(A) := \sum_{i=1}^n a_ii \mbox{ für } A := (a_{ij})_{1\le i,j\le n} \in \mathbb{R}^{n \times n}
\end{equation}
und der Definition von 
\begin{equation}
\hat{\Sigma} := \hat{B}\hat{B}^T + \hat{R}
\end{equation}
ergibt sich dann
\begin{align}
\begin{split}
\mbox{BIC} = \ln p_{\hat{\theta}}(Y) - \frac{l}{2} \ln n = \ln \left(\prod_{i=1}^n N\left(y^{(i)}; 0_m, \hat{\Sigma} \right) \right) - \frac{mk + m}{2}\ln n
\end{split}
\end{align}


# Anwendungsbeispiel
\textcolor{darkblue}{Modellvergleich | Welche Anzahl von Faktoren zwischen 1 und 5 ist am besten?}
\footnotesize
\vspace{1mm}
\underline{Beweis (fortgeführt)}

Es ergibt sich also
\setstretch{.9}
\vspace{-2mm}
\tiny
\begin{align}
\begin{split}
\mbox{BIC}
& = \sum_{i=1}^n \ln N\left(y^{(i)}; 0_m, \hat{\Sigma}  \right)  - \frac{mk + m}{2}\ln n \\
& = \sum_{i=1}^n \ln \left((2\pi)^{-\frac{m}{2}} |\hat{\Sigma}|^{-\frac{1}{2}} \exp\left(-\frac{1}{2}y^{(i)^T}\Sigma^{-1}y^{(i)} \right) \right)  - \frac{mk + m}{2}\ln n \\
& = \sum_{i=1}^n \left( -\frac{m}{2}\ln 2\pi - \frac{1}{2} \ln |\hat{\Sigma}| -\frac{1}{2}y^{(i)^T}\Sigma^{-1}y^{(i)}\right) - \frac{mk + m}{2}\ln n \\
& = -\frac{nm}{2}\ln 2\pi - \frac{n}{2} \ln |\hat{\Sigma}| - \frac{1}{2} \sum_{i=1}^n y^{(i)^T}\Sigma^{-1}y^{(i)}  - \frac{mk + m}{2}\ln n \\
& = -\frac{nm}{2}\ln 2\pi - \frac{n}{2} \ln |\hat{\Sigma}| - \frac{1}{2} \mbox{tr}\left(\Sigma^{-1}\sum_{i=1}^n y^{(i)}y^{{(i)}^T}\right)  - \frac{mk + m}{2}\ln n \\
& = -\frac{nm}{2}\ln 2\pi - \frac{n}{2} \ln |\hat{\Sigma}| - \frac{1}{2} \mbox{tr}\left(\hat{\Sigma}^{-1}YY^T\right)  - \frac{mk + m}{2}\ln n \\
\end{split}
\end{align}

# Anwendungsbeispiel
\textcolor{darkblue}{Modellvergleich | Welche Anzahl von Faktoren zwischen 1 und 5 ist am besten?}

\small 
Simulationsbasierte Validierung des BIC bei $k = 3, n = 301$ und $R := uI_m$ mit $u > 0$.

\vspace{5mm}
\footnotesize
\setstretch{.9}
```{r, echo = F, eval = F}
# mvrnorm() und inv() Pakete
library(MASS)
library(matlib)

# Paramterspezifikation
k_tbu = 3                                # Dimension des latenten Zufallsvektors
m     = 9                                # Dimension des beobachtbaren Zufallsvektors
n     = 301                              # Beobachtungsanzahl
sim   = 1e1                              # Anzahl Datensimulationen
mod   = 5                                # Anzahl Modelle mit 1,...,mod Faktoren
B     = matrix(c( 1,0,0,                 # Faktorenladungsmatrix
                  1,0,0,
                  1,0,0,
                  0,1,0,
                  0,1,0,
                  0,1,0,
                  0,0,1,
                  0,0,1,
                  0,0,1),
                 nrow  = m,
                 byrow = TRUE)
mu    = rep(0,k_tbu+m)                      # Erwartungswertparameter p_\theta(x,y)
u     = c(1e-3,1e-2,1e-1,1e0,1e1)               # u > 0
BIC   = matrix(rep(NaN, length(u)*mod),      # Simulationsgemittelte BIC Werte
               nrow = mod)


# Iterationen über Parameter
for(i in 1:length(u)){

  # Parameterspezifikation
  R      = diag(rep(u[i],m))            # Beobachtungsrauschenkovarianmatrix
  Sigma  = rbind(cbind(diag(k_tbu), t(B)),  # Kovarianzmatrixparameter p_\theta(x,y)
           cbind(B, B %*% t(B) + R))

  # Iterationen über Datensimulationen
  BICs   = matrix(rep(NaN, sim*mod),          # Parameter-spezifische BIC Werte
                  nrow = mod)

  for(s in 1:sim){
    XY  = t(mvrnorm(n,mu,Sigma))              # p_theta(x,y) sampling und z-score normalization
    Y   = XY[(k_tbu+1):nrow(XY),]               # Beobachtete Daten  Y \in R^{m x n}

    # Iterationen über Datenanalysemodelle
    for(k in 1:mod){

      # Modellschätzung und Parameterschätzer
      fa        = factanal(t(Y), factors = k, rotation = "none")
      B_hat     = fa$loadings[,1:k]
      R_hat     = diag(fa$uniquenesses)
      Sigma_hat = B_hat %*% t(B_hat) + R_hat
      BICs[k,s] = (-((n*m)/2)*log(2*pi)
                   -(n/2*log(det(Sigma_hat)))
                   -(1/2)*sum(diag(inv(Sigma_hat)%*%Y%*%t(Y)))
                   -(((m*k)+m)/2)*log(n))
    }
  }
  BIC[,i] = rowMeans(BICs)
}
```

```{r, eval = F, echo = F}
graphics.off()
fdir        =  file.path(getwd(), "5_Abbildungen")
dev.new()
par(
family      = "sans",
mfcol       = c(1,5),
pty         = "m",
bty         = "l",
lwd         = .5,
las         = 1,
mgp         = c(3,2,0),
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex.main    = 1.5)
labels = c(TeX("$u = 0.001"),
           TeX("$u = 0.01"),
           TeX("$u = 0.1"),
           TeX("$u = 1"),
           TeX("$u = 10"))
for(i in 1:5){
  barplot(
  BIC[,i] - BIC[3,i],
  names.arg = c("1", "2", "3", "4", "5"),
  xlab      = "k",
  main      = labels[i])
}


dev.copy2pdf(
file        = file.path(fdir, "mvda_5_bic_simulation.pdf"),
width       = 12,
height      =  3)
```
```{r, echo = FALSE, out.width = "100%"}
knitr::include_graphics("5_Abbildungen/mvda_5_bic_simulation.pdf")
```

\begin{center}
$\Rightarrow$ BIC erlaubt Model Recovery in Szenarien mit niedrigem bis mittlerem Beobachtungsrauschen.
\end{center}

# Anwendungsbeispiel
\textcolor{darkblue}{Modellvergleich | Welche Anzahl von Faktoren zwischen 1 und 5 ist am besten?}

\small 
BIC Modellevaluation
\vspace{1mm}

\footnotesize
\setstretch{1}
```{r, eval = F}

# Datensatz
library(lavaan)                                      # Lavaan SEM Paket
data(HolzingerSwineford1939)                         # Datensatz
Y             = t(HolzingerSwineford1939[,7:15])     # Datenmatrix
m             = nrow(Y)                              # Anzahl Tests/Variablen
n             = ncol(Y)                              # Anzahl Proband:innen
max_k         = 5                                    # Maximale Faktorenanzahl
BIC           = rep(NaN,max_k)                       # BIC

# Modell Iterationen
for(k in 1:max_k){
  
  # Modellformulierung und Modellschätzung 
  fa        = factanal(t(Y), factors = k, rotation = "none")
  B_hat     = fa$loadings[,1:k]
  R_hat     = diag(fa$uniquenesses)
  
  # Modellschätzer
  Sigma_hat = B_hat %*% t(B_hat) + R_hat
  
  # Modellevaluation
  BIC[k]     = (-((n*m)/2)*log(2*pi)
               -(n/2*log(det(Sigma_hat)))
               -(1/2)*sum(diag(inv(Sigma_hat)%*%Y%*%t(Y)))
               -(((m*k)+m)/2)*log(n))
}

```

```{r, eval = F, echo = F}
graphics.off()
fdir        =  file.path(getwd(), "5_Abbildungen")
dev.new()
par(
family      = "sans",
mfcol       = c(1,1),
pty         = "s",
bty         = "l",
lwd         = .5,
las         = 1,
mgp         = c(2,1,0),
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex.main    = 1.4)
barplot(
BIC,
names.arg = c("1", "2", "3", "4", "5"),
xlab      = "Anzahl Faktoren",
main      = "BIC")
dev.copy2pdf(
file        = file.path(fdir, "mvda_5_bic_anwendung.pdf"),
width       = 4,
height      = 4)
```
# Anwendungsbeispiel
\textcolor{darkblue}{Modellvergleich | Welche Anzahl von Faktoren zwischen 1 und 5 ist am besten?}

\small 
BIC Modellevaluation
\vspace{3mm}
```{r, echo = FALSE, out.width = "50%"}
knitr::include_graphics("5_Abbildungen/mvda_5_bic_anwendung.pdf")
```
\begin{center}
$\Rightarrow$ Das BIC Modellvergleichskriterium legt ein Modell mit 2 Faktoren nahe.
\end{center}

#
\large
\setstretch{2}
\vfill
Grundlagen

Inferenz und Lernen

Faktorenanalyse

Hauptkomponentenanalyse Revisited

Anwendungsbeispiel

**Selbstkontrollfragen**

Appendix
\vfill

# Selbstkontrollfragen

\footnotesize
\begin{enumerate}
\item Definieren Sie den Begriff des Linearen Normalverteilungsmodells (LNMs).
\item Erläutern Sie die hierarchische Form eines LNMs.
\item Geben Sie das Theorem zur LNM Datenverteilung wieder.
\item Was modellieren/erklären LNMs?
\item Definieren Sie die WDFen der gemeinsamen und marginalen Datenverteilungen eines LNM Datensatzes.
\item Erläutern Sie die Fragestellungen der Inferenz und des Lernens bei der Schätzung von LNMen.
\item Geben Sie das Theorem zum Exakten Expectation-Maximization Algorithmus wieder.
\item Erläutern Sie das Theorem zum Exakten Expectation-Maximization Algorithmus.
\item Geben Sie die Definition eines Faktorenanalysemodells wieder.
\item Erläutern Sie das Verhältnis von Datenkomponentenkovarianzen und -varianzen im Faktorenanalysemodell.
\item Geben Sie die Definition eines Hauptkomponentenanalysemodells wieder.
\item Geben Sie das Theorem zum Parameter des Hauptkomponentenanalysemodells wieder.
\item Erläutern Sie den Zusammenhang von modell-freier PCA (Einheit (4)) und modell-basierter PCA  (Einheit (5)).
\item Definieren Sie das Bayesian Information Criterion (BIC).
\item Erläutern Sie die Intuition zum BIC im Kontext von Modellvergleichen.
\end{enumerate}

#
\large
\setstretch{2}
\vfill
Grundlagen

Inferenz und Lernen

Faktorenanalyse

Hauptkomponentenanalyse Revisited

Selbstkontrollfragen

**Appendix**
\vfill

# Appendix
\footnotesize
\setstretch{1}
\begin{theorem}[Jensen's inequality]
\justifying
\normalfont
Let $x$ be a random variable and $g$ be a convex function, i.e.
\begin{equation}
g(\lambda x_1 + (1-\lambda)x_2) \le \lambda g(x_1) + (1-\lambda)g(x_2).
\end{equation}
Then
\begin{equation}
\mathbb{E}(g(x)) \ge g(\mathbb{E}(x)).
\end{equation}
Conversely, let $g$ be a concave function, i.e.
\begin{equation}
g(\lambda x_1 + (1-\lambda)x_2) \ge \lambda g(x_1) + (1-\lambda)g(x_2).
\end{equation}
Then
\begin{equation}
\mathbb{E}(g(x)) \le g(\mathbb{E}(x)).
\end{equation}
\end{theorem}
\underline{Proof}

By adapting the proof of Theorem 4.7.8 in @casella_2012 , we show the inequality for the concave case.
Let $f$ be a tangent line at the point $g(\mathbb{E}(x))$, i.e. is a linear-affine
function of the form $f(x) := ax + b$ for some $a,b\in \mathbb{R}$ with
$f(\mathbb{E}(x)) = g(\mathbb{E}(x))$. Because $g$ is concave, we have $g(x) \le ax + b$
for all $x \in \mathbb{R}$ and thus als $g(x) \le ax + b$. Hence,
\begin{equation}
\mathbb{E}(g(x)) \le \mathbb{E}(ax + b) = a\mathbb{E}(x) + b = f(\mathbb{E}(x)) = g(\mathbb{E}(x)).
\end{equation}

Remarks
\begin{itemize}
\item For convex $g$ the function's graph lies below the straight line $g(x_1)$ to $g(x_2)$.
\item For concave $g$ the function's graph lies above the straight line $g(x_1)$ to $g(x_2)$.
\item The logarithm is a concave function, hence $\mathbb{E}(\ln x) \le \ln \mathbb{E}(x)$.
\end{itemize}

# Appendix
\footnotesize
\textbf{\underline{Beweis des Theorems zum exakten EM Algorithmus eines Faktorenanalysemodells}}
\setstretch{1.2}

The E Step of the algorithm follows directly with the LGM data set inference theorem.
We thus focus on the derivation of the M Step. To this end, recall that the M Step
of the exact EM algorithm takes the form
\begin{equation}
\theta^{(j)} := \argmax_{\theta} \int p_{\theta^{(j-1)}}(X|Y) \ln p_\theta(X,Y) \,dX
\end{equation}
For LGMs, the maximization of the expected joint likelihood with respect to $\theta$
can be carried out analytically and in the sense of the necessary condition for a maximum.

To ease notation, we will write $\tilde{\theta} := \theta^{(j-1)}$ in the following
and denote the expectation of a function $f$ of the unobserved data $X$ under the conditional distribution $p_{\tilde{\theta}}(X|Y)$ as a conditional expectation:
\begin{equation}
\mathbb{E}_{\tilde{\theta}}(f(X)|Y)
:= \mathbb{E}_{p_{\tilde{\theta}}(X|Y)}(f(X))
= \int f(X)\,p_{\tilde{\theta}}(X|Y)\,dX.
\end{equation}

\setstretch{0.5}
With these simplifications, we thus aim to evaluate
\begin{equation}
\frac{\partial}{\partial B} \mathbb{E}_{p_{\tilde{\theta}}(X|Y)}\left(\ln p_\theta(X,Y)\right)
\mbox{ and }
\frac{\partial}{\partial R} \mathbb{E}_{p_{\tilde{\theta}}(X|Y)}\left(\ln p_\theta(X,Y)\right),
\end{equation}
\noindent set the results to zero, and solve for update equations for $B^{(j)}$ and $R^{(j)}$.

We proceed in four steps. (1) We first use the IID data assumption and the linearity
of conditional expectations and derivatives to simplify the problem of evaluating the
expected data set joint likelihood and its partial derivatives for a single data point
$(x^{(i)}, y^{(i)})$. We then (2) evaluate the conditional expectation and (3)
evaluate the respective partial derivatives. By capitalizing on the results form
the first step, we then (4) evaluate and simplify the ensuing parameter update equations.


# Appendix
\setstretch{1}
\footnotesize
\noindent (1) Expected joint likelihood partial derivatives under IID data assumptions
\tiny
\begin{align}
\begin{split}
\frac{\partial}{\partial B} \mathbb{E}_{p_{\tilde{\theta}}(X|Y)}\left(\ln p_\theta(X,Y)\right)
& =
\frac{\partial}{\partial B}
\left(
\mathbb{E}_{p_{\tilde{\theta}}(X|Y)}
\left(\ln \prod_{i=1}^n p_\theta\left(x^{(i)},y^{(i)}\right)\right)
\right) \\
& =
\frac{\partial}{\partial B}
\left(
\mathbb{E}_{p_{\tilde{\theta}}(X|Y)}\left(\sum_{i=1}^n \ln p_\theta\left(x^{(i)},y^{(i)}\right)\right)
\right) \\
& =
\frac{\partial}{\partial B}
\sum_{i=1}^n \mathbb{E}_{\prod_{i=1}^n p_{\tilde{\theta}}\left(x^{(i)},y^{(i)}\right)}
\left(\ln p_\theta\left(x^{(i)},y^{(i)}\right)\right)
\\
& =
\frac{\partial}{\partial B}
\sum_{i=1}^n \mathbb{E}_{p_{\tilde{\theta}}\left(x^{(i)},y^{(i)}\right)}
\left(\ln p_\theta\left(x^{(i)},y^{(i)}\right)\right) \\
& =
\sum_{i=1}^n
\frac{\partial}{\partial B}
\mathbb{E}_{p_{\tilde{\theta}}\left(x^{(i)},y^{(i)}\right)}
\left(\ln p_\theta\left(x^{(i)},y^{(i)}\right)\right),
\end{split}
\end{align}
\begin{equation}
\frac{\partial}{\partial R} \mathbb{E}_{p_{\tilde{\theta}}(X|Y)}\left(\ln p_\theta(X,Y)\right)
= \sum_{i=1}^n
\frac{\partial}{\partial R}
\mathbb{E}_{p_{\tilde{\theta}}\left(x^{(i)},y^{(i)}\right)}
\left(\ln p_\theta\left(x^{(i)},y^{(i)}\right)\right).
\end{equation}


# Appendix
\setstretch{1}
\footnotesize
\noindent (2) Expected joint likelihood for a single data point

For ease of notation, we omit the $(i)$ superscript indexing the data realizations in this step.
\vspace{-2mm}
\begin{tiny}
\begin{align*}
\begin{split}
& \mathbb{E}_{p_{\tilde{\theta}}(x|y)}\left(\ln p_\theta(x,y)\right) 		\\
& = \mathbb{E}_{\tilde{\theta}}\left(\ln p_\theta(x,y)| y \right)					\\
& = \mathbb{E}_{\tilde{\theta}}\left(\ln\left(N(y;Bx,R)N(x;0,I_k)\right)|y \right) 	\\
& = \mathbb{E}_{\tilde{\theta}}\left( \ln N(y;Bx,R) + \ln N(x;0,I_k)|y \right) 		\\
& =
\mathbb{E}\left(
\ln \left((2\pi)^{-\frac{m}{2}}|R  |^{-\frac{1}{2}}\exp\left(-\frac{1}{2}(y - Bx)^T R^{-1}(y-Bx) \right) \right) +
\ln \left((2\pi)^{-\frac{k}{2}}|I_k|^{-\frac{1}{2}}\exp\left(-\frac{1}{2} x^Tx \right)  \right)
|y\right)
\\
& =
\mathbb{E}_{\tilde{\theta}}\left(
- \frac{m+k}{2} \ln 2\pi - \frac{1}{2} \ln |R|
- \frac{1}{2}\left(y^TR^{-1}y - 2y^TR^{-1}Bx + x^TB^TR^{-1}Bx \right)
- \frac{1}{2} \ln 1 - \frac{1}{2}x^Tx
|y\right)
\\
& =
\mathbb{E}_{\tilde{\theta}}\left(
- \frac{m+k}{2} \ln 2\pi - \frac{1}{2} \ln |R|
- \frac{1}{2}y^TR^{-1}y + y^TR^{-1}Bx -\frac{1}{2} x^TB^TR^{-1}Bx - \frac{1}{2}x^Tx
|y\right)
\\
& =
\mathbb{E}_{\tilde{\theta}}\left(
- \frac{m+k}{2} \ln 2\pi - \frac{1}{2} \ln |R|
- \frac{1}{2}y^TR^{-1}y + y^TR^{-1}Bx -\frac{1}{2} \mbox{tr}\left(B^TR^{-1}Bxx^T\right) - \frac{1}{2}x^Tx
|y\right)
\\
& =
- \frac{m+k}{2} \ln 2\pi
- \frac{1}{2} \ln |R|
- \frac{1}{2}y^TR^{-1}y
+ y^TR^{-1}B\mathbb{E}_{\tilde{\theta}}(x|y)
- \frac{1}{2}\mbox{tr}\left(B^TR^{-1}B\mathbb{E}_{\tilde{\theta}}\left(xx^T|y\right)\right)
- \frac{1}{2}\mathbb{E}_{\tilde{\theta}}\left(x^Tx|y\right)
\\
\end{split}
\end{align*}
\end{tiny}
where in the 7th equality, we made use of the fact that $x^TAx = \mbox{tr}(Axx^T)$ for $x \in \mathbb{R}^n, A \in \mathbb{R}^{n \times n}$.



# Appendix
\footnotesize
\setstretch{1}
\noindent (3) Partial derivatives

To evaluate the partial derivatives of the conditional expected joint likelihood
with respect to the matrices $B$ and $R$, we require the  following identities from matrix calculus (z.B. @petersen_2012):
\begin{tiny}
\begin{align}
\begin{split}
\frac{\partial}{\partial X}A^TXB = AB^T,
\frac{\partial}{\partial X}\mbox{tr}(XA) = A^T,
\frac{\partial}{\partial X}\mbox{tr}(X^TAXB) = AXB + A^TXB^T,
\frac{\partial}{\partial X}\ln |X|				  = \left(X^{-1}\right)^T.
\end{split}
\end{align}
\end{tiny}
We then have
\begin{tiny}
\begin{align}
\begin{split}
\frac{\partial}{\partial B}\mathbb{E}_{p_{\tilde{\theta}}(x|y)}\left(\ln p_\theta(x,y)\right)
& =
\frac{\partial}{\partial B}
\left(
y^TR^{-1}B\mathbb{E}_{\tilde{\theta}}(x|y)
- \frac{1}{2}\mbox{tr}\left(B^TR^{-1}B\mathbb{E}_{\tilde{\theta}}\left(xx^T|y\right)\right)
\right)
\\
&
= \frac{\partial}{\partial B}\left(y^TR^{-1}B\mathbb{E}_{\tilde{\theta}}(x|y)\right)
- \frac{1}{2}\frac{\partial}{\partial B}\mbox{tr}\left(B^TR^{-1}B\mathbb{E}_{\tilde{\theta}}\left(xx^T|y\right)\right)
\\
&
= \left(y^TR^{-1}\right)^T\mathbb{E}_{\tilde{\theta}}(x|y)^T
- \frac{1}{2}R^{-1}B\mathbb{E}_{\tilde{\theta}}\left(xx^T|y\right)
- \frac{1}{2} \left(R^{-1}\right)^TB \mathbb{E}_{\tilde{\theta}}\left(xx^T|y\right)^T
\\
&
= R^{-1}y\mathbb{E}_{\tilde{\theta}}(x|y)^T
- \frac{1}{2}R^{-1}B\mathbb{E}_{\tilde{\theta}}\left(xx^T|y\right)
- \frac{1}{2}R^{-1}B\mathbb{E}_{\tilde{\theta}}\left(xx^T|y\right)
\\
&
=  R^{-1}y\mathbb{E}_{\tilde{\theta}}(x|y)^T
-  R^{-1}B\mathbb{E}_{\tilde{\theta}}\left(xx^T|y\right).
\end{split}
\end{align}
\end{tiny}

# Appendix
\setstretch{1}
\footnotesize
Similarly, we have
\begin{tiny}
\begin{align*}
\begin{split}
& \frac{\partial}{\partial R^{-1}}
\mathbb{E}_{p_{\tilde{\theta}}(x|y)}\left(\ln p_\theta(x,y)\right) \\
& =
\frac{\partial}{\partial R^{-1}}
\left(
- \frac{1}{2} \ln |R|
- \frac{1}{2}y^TR^{-1}y
+ y^TR^{-1}B\mathbb{E}_{\tilde{\theta}}(x|y)
- \frac{1}{2}\mbox{tr}\left(B^TR^{-1}B\mathbb{E}_{\tilde{\theta}}\left(xx^T|y\right)\right)
\right)
\\
& =
- \frac{1}{2} \frac{\partial}{\partial R^{-1}} \ln |R|
- \frac{1}{2} \frac{\partial}{\partial R^{-1}} y^TR^{-1}y
+ \frac{\partial}{\partial R^{-1}} y^TR^{-1}C\mathbb{E}_{\tilde{\theta}}(x|y)
- \frac{1}{2}\frac{\partial}{\partial R^{-1}}\mbox{tr}\left(R^{-1} C \mathbb{E}_{\tilde{\theta}}\left(xx^T|y\right)B^T \right) \\
& =
\frac{1}{2} R
- \frac{1}{2} yy^T
+ y\left(B\mathbb{E}_{\tilde{\theta}}(x|y)\right)^T
- \frac{1}{2}\left(B \mathbb{E}_{\tilde{\theta}}\left(xx^T|y\right)B^T \right)^T \\
& =
\frac{1}{2} R
- \frac{1}{2} yy^T
+ y \mathbb{E}_{\tilde{\theta}}(x|y)^T B^T
- \frac{1}{2}B \mathbb{E}_{\tilde{\theta}}\left(xx^T|y\right)B^T.
\end{split}
\end{align*}
\end{tiny}


# Appendix
\setstretch{1}
\footnotesize
\noindent (4) Parameter update equations

Re-substitution then yields for the partial derivative with respect to $B$
\begin{tiny}
\begin{align}
\begin{split}
\frac{\partial}{\partial B} \mathbb{E}_{p_{\tilde{\theta}}(X|Y)}\left(\ln p_\theta(X,Y)\right)
& =
\sum_{i=1}^n
\frac{\partial}{\partial B}
\mathbb{E}_{p_{\tilde{\theta}}\left(x^{(i)},y^{(i)}\right)}
\left(\ln p_\theta\left(x^{(i)},y^{(i)}\right)\right) \\
& =
\sum_{i=1}^n R^{-1}y^{(i)}\mathbb{E}_{\tilde{\theta}}\left(x^{(i)}|y^{(i)}\right)^T
- R^{-1}B\mathbb{E}_{\tilde{\theta}}\left(x^{(i)}x^{{(i)}^T}|y^{(i)}\right)  \\
& =
R^{-1}\sum_{i=1}^n y^{(i)}\mathbb{E}_{\tilde{\theta}}\left(x^{(i)}|y^{(i)}\right)^T
- R^{-1}B \sum_{i=1}^n \mathbb{E}_{\tilde{\theta}}\left(x^{(i)}x^{{(i)}^T}|y^{(i)}\right).
\end{split}
\end{align}
\end{tiny}

# Appendix
\setstretch{1}
\footnotesize
Setting to zero and solving for $B^{(j)}$ then yields
\begin{tiny}
\begin{align}
\begin{split}
R^{-1}B^{(j)} \sum_{i=1}^n \mathbb{E}_{\tilde{\theta}}\left(x^{(i)}x^{{(i)}^T}|y^{(i)}\right)
& = R^{-1}\sum_{i=1}^n y^{(i)}\mathbb{E}_{\tilde{\theta}}\left(x^{(i)}|y^{(i)}\right)^T.
\\
\Leftrightarrow
B^{(j)}
& = \sum_{i=1}^n y^{(i)}\mathbb{E}_{\tilde{\theta}}\left(x^{(i)}|y^{(i)}\right)^T
\left(\sum_{i=1}^n \mathbb{E}_{\tilde{\theta}}\left(x^{(i)}x^{{(i)}^T}|y^{(i)}\right)\right)^{-1}.
\end{split}
\end{align}
\end{tiny}

# Appendix
\setstretch{1}
\footnotesize
Similarly, re-substitution yields for the partial derivative with respect to $R$
\begin{tiny}
\begin{align}
\begin{split}
& \frac{\partial}{\partial R^{-1}}
\mathbb{E}_{p_{\tilde{\theta}}(X|Y)}\left(\ln p_\theta(X,Y)\right) \\
& = \sum_{i=1}^n
\frac{\partial}{\partial R^{-1}}
\mathbb{E}_{p_{\tilde{\theta}}\left(x^{(i)},y^{(i)}\right)}
\left(\ln p_\theta\left(x^{(i)},y^{(i)}\right)\right)  \\
& = \sum_{i=1}^n
\frac{1}{2} R
- \frac{1}{2} y^{(i)}y^{(i)^T}
+ y^{(i)}\mathbb{E}_{\tilde{\theta}}\left(x^{(i)}|y^{(i)}\right)^T B^T
- \frac{1}{2}B \mathbb{E}_{\tilde{\theta}}\left(x^{(i)}x^{(i)^T}|y^{(i)}\right)B^T \\
& =
\frac{n}{2} R
- \frac{1}{2} \sum_{i=1}^n y^{(i)}y^{(i)^T}
+ \sum_{i=1}^n y^{(i)}\mathbb{E}_{\tilde{\theta}}\left(x^{(i)}|y^{(i)}\right)^T B^T
- \frac{1}{2}B \sum_{i=1}^n \mathbb{E}_{\tilde{\theta}}\left(x^{(i)}x^{(i)^T}|y^{(i)}\right)B^T.
\end{split}
\end{align}
\end{tiny}

# Appendix
\footnotesize
\underline{Beweis (fortgeführt)}

Setting to zero and solving for $R^{(j)}$ then yields
\begin{tiny}
\begin{align*}
\begin{split}
\frac{n}{2} R^{(j)}
&
= \frac{1}{2} \sum_{i=1}^n y^{(i)}y^{(i)^T}
- \sum_{i=1}^n y^{(i)}\mathbb{E}_{\tilde{\theta}}\left(x^{(i)}|y^{(i)}\right)^T B^T
+ \frac{1}{2}B \sum_{i=1}^n \mathbb{E}_{\tilde{\theta}}\left(x^{(i)}x^{(i)^T}|y^{(i)}\right)B^T \\
R^{(j)}
&
= \frac{1}{n} \sum_{i=1}^n y^{(i)}y^{(i)^T}
- \frac{2}{n}\sum_{i=1}^n y^{(i)}\mathbb{E}_{\tilde{\theta}}\left(x^{(i)}|y^{(i)}\right)^T B^T
+ \frac{1}{n}B \sum_{i=1}^n \mathbb{E}_{\tilde{\theta}}\left(x^{(i)}x^{(i)^T}|y^{(i)}\right)B^T
\end{split}
\end{align*}
\end{tiny}

# Appendix
\setstretch{1}
\footnotesize
Substitution of the update equation for $B$ further yields
\begin{tiny}
\begin{align*}
\begin{split}
R^{(j)}
=
&    \frac{1}{n} \sum_{i=1}^n y^{(i)}y^{(i)^T}
  - \frac{2}{n}\sum_{i=1}^n y^{(i)}\mathbb{E}_{\tilde{\theta}}\left(x^{(i)}|y^{(i)}\right)^T B^{(j)^T} \\
& + \frac{1}{n}\sum_{i=1}^n y^{(i)}\mathbb{E}_{\tilde{\theta}}\left(x^{(i)}|y^{(i)}\right)^T
    \left(\sum_{i=1}^n \mathbb{E}_{\tilde{\theta}}\left(x^{(i)}x^{{(i)}^T}|y^{(i)}\right)\right)^{-1}
    \sum_{i=1}^n \mathbb{E}_{\tilde{\theta}}\left(x^{(i)}x^{(i)^T}|y^{(i)}\right)B^{(j)^T} \\
=
&
   \frac{1}{n} \left(
    \sum_{i=1}^n y^{(i)}y^{(i)^T}
 - 2\sum_{i=1}^n y^{(i)}\mathbb{E}_{\tilde{\theta}}\left(x^{(i)}|y^{(i)}\right)^T B^{(j)^T}
 + \sum_{i=1}^n y^{(i)}\mathbb{E}_{\tilde{\theta}}\left(x^{(i)}|y^{(i)}\right)^T B^{(j)^T}
 \right) \\
=
&
\frac{1}{n}
\left(
\sum_{i=1}^n y^{(i)}y^{(i)^T}
- \sum_{i=1}^n y^{(i)}\mathbb{E}_{\tilde{\theta}}\left(x^{(i)}|y^{(i)}\right)^T B^{(j)^T}
\right)
\end{split}
\end{align*}
\end{tiny}

# Appendix
\footnotesize
\setstretch{1}
We thus obtained the parameter update equations
\begin{tiny}
\begin{align*}
\begin{split}
B^{(j)}
& = \sum_{i=1}^n y^{(i)}\mathbb{E}_{\tilde{\theta}}\left(x^{(i)}|y^{(i)}\right)^T
\left(\sum_{i=1}^n \mathbb{E}_{\tilde{\theta}}\left(x^{(i)}x^{{(i)}^T}|y^{(i)}\right)\right)^{-1} \\
R^{(j)}
& =
\frac{1}{n}
\left(
\sum_{i=1}^n y^{(i)}y^{(i)^T}
- \sum_{i=1}^n y^{(i)}\mathbb{E}_{\tilde{\theta}}\left(x^{(i)}|y^{(i)}\right)^T B^{(j)^T}
\right)
\end{split}
\end{align*}
\end{tiny}
\noindent The computational forms of these update equations can be further simplified by noting that with
\begin{equation}
\mathbb{E}\left(xx^T|y\right) = \mu_{x|y}\mu_{x|y}^T + \Sigma_{x|y}
\end{equation}
and the LGM inference theorem it holds that
\begin{tiny}
\begin{align}
\mathbb{E}_{\tilde{\theta}}\left(x^{(i)}|y^{(i)}\right)
= \hat{x}^{(i)} \mbox{ and }
\mathbb{E}_{\tilde{\theta}}\left(x^{(i)}x^{{(i)}^T}|y^{(i)}\right)
= \hat{x}^{(i)}\hat{x}^{(i)^T} + \hat{\Sigma}^{(i)}.
\end{align}
\end{tiny}

# Appendix
\setstretch{1}
\footnotesize
We thus obtain
\begin{tiny}
\begin{align*}
\begin{split}
B^{(j)} & = \sum_{i=1}^n y^{(i)}\hat{x}^{(i)^T}\left(\sum_{i=1}^n \hat{x}^{(i)}\hat{x}^{(i)^T} + \hat{\Sigma}^{(i)}\right)^{-1} \\
R^{(j)}
& =
\frac{1}{n}\left(\sum_{i=1}^n y^{(i)}y^{(i)^T} - \sum_{i=1}^n y^{(i)}\hat{x}^{(i)^T} B^{(j)^T}\right).
\end{split}
\end{align*}
\end{tiny}
Finally, enforcing the diagonality constraint on $R$ can be achieved by setting
\begin{tiny}
\begin{equation}
R^{(j)}
:= - \frac{1}{n}\mbox{diag}\left(\sum_{i=1}^n y^{(i)}y^{(i)^T} + \sum_{i=1}^n y^{(i)}\hat{x}^{(i)^T} B^{(j)^T}\right)
\end{equation}
\end{tiny}



# Appendix
\footnotesize
\setstretch{1.2}
\textbf{\underline{Beweis des Theorems zur direkten Maximum Marginal Likelihood Schätzung des probabilistischen PCA Modells}}

We only show that $\hat{B}$ as defined in the theorem corresponds to maximum
of the log marginal likelihood function. To this end, we closely follow the
respective proof in @tipping_1999, and proceed in three steps:
(1) We first rewrite the marginal data set log likelihood function of the PPCA
model in a suitable manner. (2) We then evaluate its gradient with respect to $B$
and (3) finally evaluate the resulting maximum marginal likelihood estimator.

\noindent (1) Log likelihood function

We first rewrite the marginal data set log likelihood function of the PPCA model
\begin{equation}
\ell : \Theta \to \mathbb{R}, \theta \mapsto \ell(\theta)
:= \ln p_\theta(Y) = \sum_{i=1}^n  \ln N\left(y^{(i)}; 0_m, BB^T + \sigma^2 I_m \right)
\end{equation}
with $\theta := \{B, \sigma^2\}$ in a way more amenable to direct maximization.

With the definitions of
\begin{equation}
\Sigma := BB^T + \sigma^2 I_m \mbox{ and }
C := \frac{1}{n} \sum_{i=1}^n y^{(i)} y^{(i)^T}
\end{equation}
and the trace operator properties
\begin{equation}
x^TAx = \mbox{tr}(Axx^T)
\mbox{ and } \mbox{tr}(A) + \mbox{tr}(B) = \mbox{tr}(A + B),
\end{equation}
we have



# Appendix
\tiny
\begin{align}
\begin{split}
\ln p_\theta(Y)
& = \sum_{i=1}^n  \ln N\left(y^{(i)}; 0_m, \Sigma \right) \\
& = \sum_{i=1}^n  \ln \left((2\pi)^{-\frac{m}{2}}|\Sigma|^{-\frac{1}{2}}
                  \exp\left(-\frac{1}{2}y^{(i)^T} \Sigma^{-1} y^{(i)}\right)\right) \\
& = -\frac{n}{2}\left(m\ln 2 \pi + \ln |\Sigma| + \sum_{i=1}^n \frac{1}{n} y^{(i)^T} \Sigma^{-1} y^{(i)} \right) \\
& = -\frac{n}{2}\left(m\ln 2 \pi + \ln |\Sigma| + \sum_{i=1}^n \mbox{tr}\left(\Sigma^{-1} \frac{1}{n}y^{(i)} y^{(i)^T} \right)\right) \\
& = -\frac{n}{2}\left(m\ln 2 \pi + \ln |\Sigma| + \mbox{tr}\left(\sum_{i=1}^n \Sigma^{-1} \frac{1}{n} y^{(i)} y^{(i)^T} \right)\right) \\
& = -\frac{n}{2}\left(m\ln 2 \pi + \ln |\Sigma| + \mbox{tr}\left(\Sigma^{-1} \frac{1}{n}\sum_{i=1}^n   y^{(i)} y^{(i)^T} \right)\right) \\
& = -\frac{n}{2}\left(m\ln 2 \pi + \ln |\Sigma| + \mbox{tr}\left(\Sigma^{-1} C\right) \right).
\end{split}
\end{align}



# Appendix
\footnotesize
\noindent (2) Gradient of the log marginal likelihood function

With
\begin{tiny}
\begin{equation}
\frac{\partial}{\partial X} \ln |X| = \left(X^{-1}\right)^T,
\frac{\partial}{\partial X} XX^T = 2X, \mbox{ and }
\frac{\partial}{\partial X} \mbox{tr}(XA) = A^T,
\end{equation}
\end{tiny}
\noindent the gradient of the log marginal likelihood function with respect to $B$ evaluates to
\begin{tiny}
\begin{align}
\begin{split}
\frac{\partial}{\partial B} \ell(\theta)
& = -\frac{n}{2}\frac{\partial}{\partial B}
\left(m \ln 2 \pi + \ln |BB^T + \sigma^2 I_m| + \mbox{tr}\left(\left(BB^T + \sigma^2 I_m \right)^{-1} C\right)\right) \\
& =
-\frac{n}{2}\frac{\partial}{\partial B} \ln |BB^T + \sigma^2 I_m|
-\frac{n}{2}\frac{\partial}{\partial B} \mbox{tr}\left(\left(BB^T + \sigma^2 I_m \right)^{-1} C\right) \\
& =
-\frac{n}{2} 2 \left(\left(BB^T + \sigma^2 I_m \right)^{-1} B \right)^T
+\frac{n}{2} 2\left(BB^T + \sigma^2 I_m \right)^{-1} C \left(BB^T + \sigma^2 I_m \right)^{-1} B \\
& = n\left(- \Sigma^{-1}B + \Sigma^{-1}C\Sigma^{-1}B \right)\\
& = n\left(\Sigma^{-1}C\Sigma^{-1}B - \Sigma^{-1}B \right).\\
\end{split}
\end{align}
\end{tiny}

# Appendix
\footnotesize
\noindent (3) Maximum marginal likelihood estimator evaluation
\setstretch{1}

Setting the gradient of $\ell$ with respect to $B$ to zero then yields
\begin{equation}
\Sigma^{-1}C\Sigma^{-1}\hat{B} - \Sigma^{-1}\hat{B} = 0
\Leftrightarrow
\Sigma^{-1}\hat{B} = \Sigma^{-1}C\Sigma^{-1}\hat{B}
\Leftrightarrow
\hat{B} = C\Sigma^{-1}\hat{B} .
\end{equation}
We consider solutions of this necessary condition for a stationary point of the
log marginal likelihood function with $B \neq 0$ and $\Sigma \neq C$. To find these,
we first express $\hat{B}$ in terms of its singular value decomposition
\begin{equation}
\hat{B} = ULV^T
\end{equation}
where $U = (u_1, u_2, ..., u_l)$ is an $m \times q$ matrix of orthonormal column
vectors, $L = \mbox{diag}(l_1,l_2,...,l_q)$ is a $q \times q$ diagonal matrix of
singular values, and $V$ is a $q \times q$ orthogonal matrix. Substitution in
the necessary condition for a stationary point then yields
\begin{equation}\label{eq:cul}
CUL = U(\sigma^2 I_m + L^2)L.
\end{equation}
For $l_j \neq 0$, eq. \eqref{eq:cul} implies that
\begin{equation}
Cu_j = (\sigma^2 + l_j^2)u_j
\end{equation}
Hence, each column of $U$ must be an eigenvector of $C$ with corresponding
eigenvalue $\lambda_j = \sigma^2 + l_j$, and thus
\begin{equation}
\lambda_j = \sigma^2 + l_j^2
\Leftrightarrow
l_j^2  = \lambda_j - \sigma^2
\Leftrightarrow
l_j = (\lambda_j - \sigma^2)^{\frac{1}{2}}.
\end{equation}
For $l_j = 0$, $u_j$ is arbitrary. Under the assumption that $l_j \neq 0$ for
$j = 1,...,m$, all potential solutions for $\hat{B}$ can thus be written in the form
\begin{equation}
\hat{B} = U_q\left(\Lambda_q - \sigma^2 I_m\right)^{\frac{1}{2}}R,
\end{equation}
where $U_q$ is a $m \times q$ matrix whose $q$ columns are the eigenvectors of
$C$, $\Lambda_q$ is the diagonal matrix of the corresponding eigenvalues, and $R$
is an arbitrary $q \times q$ orthogonal matrix, for example, $R = I_q$.

$\hfill \Box$

# References {.allowframebreaks}
\footnotesize

