# A Contemporary Overview of Probabilistic Latent Variable Models

Rick Farouni<EMAIL_ADDRESS>
Department of Molecular Pathology,
Massachusetts General Hospital,
Harvard Medical School

###### Abstract

In this paper we provide a conceptual overview of latent variable models
within a probabilistic modeling framework, an overview that emphasizes the
compositional nature and the interconnectedness of the seemingly disparate
models commonly encountered in statistical practice.

Keywords: Latent Variable Models, Deep Generative Models, Probabilistic Models

##  1 Introduction

In this paper, we will be concerned with building probabilistic latent
variable models for data that populate an $N\times P$ matrix $\mathbf{Y}$
whose every element $y_{np}\in\operatorname{\mathbb{R}}$ is real-valued
measurement.

$\displaystyle\mathbf{Y}=\begin{bmatrix}\mathbf{y}_{1}^{\rm T}\\\ \vdots\\\
\mathbf{y}_{n}^{\rm T}\\\ \vdots\\\ \mathbf{y}_{N}^{\rm
T}\end{bmatrix}=\begin{bmatrix}y_{1,1}&\dots&y_{1,P}\\\
\vdots&\ddots&\vdots\\\ y_{n,1}&\dots&y_{n,P}\\\ \vdots&\ddots&\vdots\\\
y_{N,1}&\dots&y_{N,P}\end{bmatrix}$ (1)

In general, we will be limiting our focus mainly to the the multivariate
statistics setting in which we treat an observation as a multivariate random
vector $\mathbf{y}_{n}=\left(y_{n,1},\ldots,y_{n,P}\right)$ consisting of $P$
features and the data as a set of $N$ observations
$\mathbf{y}=\\{\mathbf{y}_{1},\cdots,\mathbf{y}_{n},\cdots,\mathbf{y}_{N}\\}$.
Accordingly, we can think of the data matrix $\mathbf{Y}$ as consisting of $N$
observations, each a $P$-dimensional vector, that populate the rows of the
matrix.

##  2 Basic Definitions

In order to give a rigorous constructive account of probabilistic latent
variable models, we begin with defining several foundational concepts that
will later serve as our building blocks. In what follows, we attempt to convey
the intuition at the expense of mathematical rigor. More rigorous definitions
can be found in Durrett (2010) and Athreya and Lahiri (2006).

###  2.1 What is a Probability Model?

Our first building block is a probability model. A _probability model_ is
defined as the triplet consisting of a set, a set of sets, and a function.
More precisely given by

$\left(\mathcal{Y},\ \mathcal{B}(\mathcal{Y}),\ p\right)$

where

  1. 1.

$\mathcal{Y}$ is the _sample space_

  2. 2.

$\mathcal{B}\left(\mathcal{Y}\right)$ is the _Borel $\sigma$-algebra_

  3. 3.

$p$ is a _probability measure_

Whereas the sample space is a relatively easy concept to understand - it is a
set of all possible _outcomes_ \- the notion of a Borel $\sigma$-algebra tends
to be more involved. So what exactly is a Borel $\sigma$-algebra? The short
answer is that it is the set of all _events_ of interest; technically, the set
of subsets of $\mathcal{Y}$. For the long answer, we first need to define an
_algebra_. An algebra in this context is a set along with a collection of
operations. Here is an example. First take our sample space, the set
$\mathcal{Y}$, and equip it with two binary operations
$\left(\cup,\cap\right)$, the union and intersection. Now let us consider two
events only - the empty set $\emptyset$ denoting the event that nothing
happens and the sample space $\mathcal{Y}$ denoting the event that all
outcomes enumerated in the sample space occur. Starting with these two events
only, we can, using the two operations, construct the set of all possible
subsets of $\mathcal{Y}$, which is called the powerset

$\big{(}\mathcal{P}\left(\mathcal{Y}\right),\cap,\cup,\neg,\emptyset,\mathcal{Y}\big{)}$

Since every possible subset of $\mathcal{Y}$ is an event, accordingly the
powerset consists of all possible events. With a finite number of outcomes,
the number of all possible events is finite as well. Therefore, the powerset
is an algebra that is closed under the union and intersection of finitely many
subsets. However, when our sample space is not finite, but rather countably
infinite, we would instead need a specialized set of subsets to enumerate all
possible events. That set is called a $\sigma$-algebra. It is simply a sub-
algebra of the powerset, completed to include countably infinite operations.
Lastly, when the events are subsets of the real line $\mathbb{R}$, we would
need to resort to a more specialized algebra, called the Borel
$\sigma$-algebra. It is defined as the $\sigma$-algebra generated by the open
sets and which can be constructed from open or closed sets by repeatedly
taking countable unions and intersections. In notation, we define it as such.

$\mathcal{B}(\mathbb{R}):=\sigma(\mathcal{C});\quad\mathcal{C}=\\{(a,b],-\infty\leq
a\leq b\leq\infty\\}$

Whereas the powerset is the largest possible $\sigma$-algebra, the collection
of subsets $\\{\emptyset,\mathcal{Y}\\}$ is the smallest. The Borel
$\sigma$-algebra lies between the two extremes

$\\{\emptyset,\mathcal{Y}\\}\subset\mathcal{B}(\mathcal{Y})\subset\mathcal{P}(\mathcal{Y})$

and it is the smallest $\sigma$-algebra containing all open sets in
$\operatorname{\mathbb{R}}$ (hyper-rectangles in
$\operatorname{\mathbb{R}}^{P}$).

The third component of a probability model involves the notion of a measure.
What is a probability measure? It is simply a function
$\mathit{p}:\mathcal{B}\left(\mathcal{Y}\right)\to[0,1]$ that maps events to a
real valued number between zero and one, thus assigning each event a finite
probability of occurring. The notion of a probability measure is tightly
intertwined with the concept of a random variable. To illustrate, our
multivariate observation $\mathbf{y}=\left(y_{1},\ldots,y_{P}\right)$ is a
real-valued continuous random vector with distribution $\mu$ that takes values
in the _measurable space_ $\big{(}\mathcal{Y},\
\mathcal{B}(\mathcal{Y})\big{)}$. Given the probability triplet, we can define
a random vector (with a slight abuse of notation) as a measurable function
from the sample space to a $P$-dimensional real space
$\mathbf{y}:\mathcal{Y}\mapsto\operatorname{\mathbb{R}}^{P}$ such that

$\\{\mathbf{y}\leq\textrm{y}\\}\in B;\ where\quad
B\in\mathcal{B}\left(\operatorname{\mathbb{R}}^{P}\right)$

The important point is that the distribution of $\mathbf{y}$ induces a
probability measure on the measurable space $\big{(}\mathcal{Y},\
\mathcal{B}(\mathcal{Y})\big{)}$ and is given by

$\mu(B)=p\left(\mathbf{y}\in B\right);\ where\quad
B\in\mathcal{B}\left(\operatorname{\mathbb{R}}^{P}\right)$

where $p\left(\mathbf{y}\in B\right)$ denotes the probability that the event
$\\{\mathbf{y}\leq\textrm{y}\\}$ lies in a Borel set $B$. Note that a _Borel
set_ $B$ is an element of a Borel $\sigma$-algebra.

###  2.2 What is a Statistical Model?

We build a probabilistic model starting with the joint distribution of all
observable and unobservable variables. Since we have not introduced any
unobservable quantities yet, our entry point is the following joint
probability distribution for all $N$ observations,

$\displaystyle
p\left(\mathbf{y}\right)=p\left(\mathbf{y}_{1},\mathbf{y}_{2},\cdots,\mathbf{y}_{N}\right)$
(2)

where $p$ is a joint probability density function corresponding to a specific,
known probability measure, or equivalently, to a known random variable.

Technically, the joint density given by Equation 2 denotes a probability
model, not a statistical model. To illustrate, let $\mathcal{M}(\mathcal{Y})$
denote the space of all probability measures on the measurable space
$\big{(}\mathcal{Y},\ \mathcal{B}(\mathcal{Y})\big{)}$. As can be seen, there
are infinitely many possible probability measures we can choose from. That is,
$\mathcal{M}(\mathcal{Y})$ can be too large a space for many inference
problems with limited data. To make inference a less daunting task we can
restrict ourselves to a indexed subset of the space by introducing an
unobservable variable, a parameter $\theta\in\Theta$ that represents the
pattern that explains the data, where the parameter space $\Theta$ is the set
of all possible values of $\theta$. As a result, the probability measures
$P_{\theta}$ are now elements of $\mathcal{PM}(\mathcal{Y})$, the space of all
probability measures on $\Theta$ with elements
$P_{\theta}\in\mathcal{PM}(\mathcal{Y})$ indexed by a parameter
$\theta\in\Theta$.

Accordingly, a statistical model $\mathcal{P}$ can be defined as a subset of
the space of all probability measures
$\mathcal{P}\subset\mathcal{M}(\mathcal{Y})$ such that
$\mathcal{P}=\\{P_{\theta}\mid\theta\in\Theta\\}$ where $\theta\to P_{\theta}$
is a bijective and measurable assignment. We say that the model $\mathcal{P}$
is parametric statistical model if $\Theta\subset\mathbb{R}^{d}$ for some
$d\in\mathbb{N}$. The subset $\mathcal{P}\subset\mathcal{PM}(\mathcal{Y})$ can
be even restricted further. We can specify a family of parametric models
$\mathcal{G}=\\{G_{\theta}\mid\theta\in\Theta\\}$ where $\theta\to G_{\theta}$
is smooth. For example, $\mathcal{G}=\\{N(\theta,1):\theta\in\Theta\\}$
specifies the one-dimensional normal location family of models. If on the
other hand we choose the parameter space $\Theta$ to be infinite dimensional,
then we say that $\mathcal{P}$ is a nonparametric statistical model such that
$\Theta$ is equivalent to $\mathcal{PM}(\mathcal{Y})$, the space of all
probability measures on $\Big{(}\mathcal{Y},\mathcal{B}(\mathcal{Y})\Big{)}$.
As a result, we obtain the following hierarchy of probability model spaces.

$\mathcal{G}\subset\mathcal{P}\subseteq\mathcal{PM}\subset\mathcal{M}$

Now, by introducing the parameter $\boldsymbol{\theta}$, we can now specify
the joint density of a parametric statistical model for our multivariate
observations as such

$\displaystyle
p\left(\mathbf{y};\boldsymbol{\theta}\right)=p\left(\mathbf{y}_{1},\mathbf{y}_{2},\cdots,\mathbf{y}_{N};\boldsymbol{\theta}\right)$
(3)

where $\boldsymbol{\theta}\in\Theta\subset\mathbb{R}^{d}$. In the case of a
parametric statistical model, the superscript $d$ is a finite integer, whereas
in the case of a nonparametric statistical model, it is an infinite integer.
Note that we use the semicolon to denote that the parameter
$\boldsymbol{\theta}$ is not a random variable.

###  2.3 What is a Bayesian Model?

In frequentist statistics, we treat a parameter $\boldsymbol{\theta}$ as an
unobservable unknown quantity and estimation involves finding
$\boldsymbol{\theta}^{\ast}\in\Theta$ that governs the distribution of the
data we observe if the true generating model indeed lies inside the family of
probability models we specify. In the probabilistic perspective, we treat a
parameter as an unobservable random quantity whose prior distributions
quantifies our initial uncertainty. Learning involves computing the posterior
distribution of $\boldsymbol{\theta}$ that quantifies the residual uncertainty
that results after we observe the data.

Accordingly, the joint probability density of all observable and unobservable
random variables can now be given by

$\displaystyle
p\left(\mathbf{y},\boldsymbol{\theta}\right)=p\left(\mathbf{y}_{1},\mathbf{y}_{2},\cdots,\mathbf{y}_{N},\boldsymbol{\theta}\right)$
(4)

Note that Equation 4 defines a general probabilistic model, not a Bayesian
model. To reduce it to a Bayesian model, we need to factorize the joint into
the product of a conditional data distribution and a marginal distribution
prior distribution as such.

$\displaystyle
p\left(\mathbf{y},\theta\right)=\prod_{n=1}^{N}p\left(\mathbf{y}_{n}\mid\boldsymbol{\theta}\right)p\left(\boldsymbol{\theta}\right)$
(5)

This factorization entails a particular set of assumptions. To illustrate
them, we need to define a Bayesian model first.

A parametric Bayesian statistical model ($\mathcal{P},\Pi$) consists of a
model $\mathcal{P}$, the data distribution, and a prior distribution $\Pi$ on
$\Theta$ such that $\boldsymbol{\theta}$ is a random variable taking values in
$\Theta$ given that
$\Pi\Big{(}\\{P_{\boldsymbol{\theta}}:\boldsymbol{\theta}\in\Theta\\}\Big{)}=1$.
A Bayesian model is a generative model in which the data is generated
hierarchically in two stages:

$\displaystyle\begin{split}\boldsymbol{\theta}&\sim\Pi\\\
\mathbf{y}_{n}\mid\boldsymbol{\theta}&\overset{i.i.d.}{\sim}P_{\boldsymbol{\theta}}\qquad\boldsymbol{\theta}\in\Theta\subset\mathbb{R}^{d}\end{split}$
(6)

A Bayesian model can be thought of as a random mixture model where we first
sample from a mixing measure $\boldsymbol{\theta}\sim\Pi$, then sample from a
component $\mathbf{y}_{n}\mid\boldsymbol{\theta}\sim P_{\boldsymbol{\theta}}$.
After we observe the data, we then update the prior to the posterior
$\Pi(.\mid\mathbf{y}_{1},\mathbf{y}_{2},\cdots)$.

In contrast, a nonparametric Bayesian model is a Bayesian model whose prior
$\Pi$ is defined on an infinite dimensional parameter space $\Theta$. The
corresponding two-stage hierarchical model can be given as such.

$\displaystyle\begin{split}P&\sim\Pi\\\ \mathbf{y}_{n}\mid P&\sim_{iid}P\qquad
P\in\mathcal{P}\end{split}$ (7)

A prior distribution on an infinite dimensional space is a stochastic process.
Defining an infinite dimensional prior distributions is not straightforward,
but one way to construct a prior distribution $\Pi$ on $\Theta$ is through De
Finetti’s Theorem.

######  Theorem 1 (De Finetti’s Theorem)

A sequence of random variables $\\{\mathbf{y}_{n}\\}_{n=1}^{\infty}$ with
values on $\mathcal{Y}$ is exchangeable if and only if there is a unique
measure $\Pi$ on $\Theta$ such that for all $N$

$\displaystyle
p(\mathbf{y}_{1},\mathbf{y}_{2},...,\mathbf{y}_{N})=\int_{\Theta}\left(\prod_{n=1}^{N}p(\mathbf{y}_{n}\mid\theta)\right)p(\theta)d\theta$
(8)

De Finetti’s gives us a infinite mixture representation of the joint
probability of the observations (Ghosh and Ramamoorthi, 2006). More
importantly, the theorem implies that there exists a random probability
measure $P$ such that the sequence of exchangeable observations are $i.i.d$
conditional on $P$. Implicit in the resulting mixture representation is a
data-generating process that gives rise to the observations.

What is exchangeability? We say a sequence of random variables is exchangeable
when their distribution is invariant under permutation of the indices. When
building statistical models, we can use exchangeability and De Finetti’s
Theorem to treat an infinite exchangeable sequence as a conditionally $i.i.d$
sequence of random variables. In data applications, exchangeability implies
that the future is like the past.

##### Distributional Symmetries

It is important to highlight here that exchangeability is one of the four
basic probabilistic symmetries defined by Kallenberg (2006); namely,
stationarity, contractability, exchangeability, and rotatability. These four
distributional symmetries correspond to invariance properties under the four
respective transformations of shifts, sub-sequences, permutations, and
isometries. The four probabilistic symmetries form a hierarchy with
stationarity at the bottom and rotatability, the strongest invariance
property, at the top. Exchangeability in particular allows us, through
conditioning, to transform the joint density into a mixture of densities.

##### Example of a Bayesian model

When $\mathcal{Y}$ is finite, there is usually a natural unique measure $\Pi$
we can obtain. More specifically, if $\mathcal{Y}=\\{1,2,\ldots,K\\}$, then
$\mathcal{PM}(\mathcal{Y})=\\{(p_{1},\ldots,p_{K}):0\leq p_{k}\leq 1,\sum
p_{k}=1\\}$. That is, the space of probability measures corresponds to a
simplex parametrized by a $K-1$ dimensional vector
$\boldsymbol{p}=(p_{1},\ldots,p_{K-1})$. A natural prior $\Pi$ to specify on
$\boldsymbol{p}$ is the Dirichlet distribution. For example, consider the
Bayesian model($\mathcal{P},\Pi$) where the observation model $\mathcal{P}$ is
the _Categorical distribution_ defined on the sample space
$\mathcal{Y}=\\{1,2,\ldots,K\\}$, and the prior $\Pi$ is the _Dirichlet
distribution_ defined on the simplex $\Theta=\\{(p_{1},\ldots,p_{K}):0\leq
p_{k}\leq 1,\sum p_{k}=1\\}$

$\displaystyle\boldsymbol{p}$
$\displaystyle\sim\operatorname{Dirichlet}(\boldsymbol{\alpha})$
$\displaystyle y_{i}\mid\boldsymbol{p}$
$\displaystyle\sim\operatorname{Categorical}(\boldsymbol{p})$
$\displaystyle\textit{where}\quad\boldsymbol{\alpha}$
$\displaystyle=\left(\frac{\alpha}{K},\cdots,\frac{\alpha}{K}\right)$
$\displaystyle\boldsymbol{p}$ $\displaystyle=\left(p_{1},\cdots,p_{K}\right)$

Figure 1: Dirichlet Distributions

Here, the Dirichlet distribution is the conjugate prior of the Categorical
distribution and the concentration hyperparameter vector $\boldsymbol{\alpha}$
represents the number of pseudo-observations, the a priori weights, for each
of the $K$ clusters. Figure 1 shows density and sampling plots corresponding
to Dirichlet distributions for $K=3$ and several choices of
$\boldsymbol{\alpha}$.

The Dirichlet distribution is a conjugate prior that allows us to sample $iid$
from $\mathcal{Y}=\\{1,2,\ldots,K\\}$. To generalize, we can make
$\mathcal{Y}=\operatorname{\mathbb{R}}$ thus obtaining the space of all
measures $\mathcal{M}(\operatorname{\mathbb{R}})$ defined on the measurable
space $(\operatorname{\mathbb{R}},\mathcal{B})$, which is the real line
equipped with the Borel $\sigma$-algebra. It can be shown (Ghosh and
Ramamoorthi, 2006) that for every partition $B_{1},B_{2},\ldots,B_{k}$ of the
real line $\operatorname{\mathbb{R}}$ by Borel sets, there exists a unique
measure $DP_{\alpha}$ on $\mathcal{M}(\operatorname{\mathbb{R}})$ called the
Dirichlet Process with parameter $\alpha$ satisfying

$(P(B_{1}),P(B_{2}),\ldots,P(B_{k}))=DP(\alpha B_{1},\alpha
B_{2},\ldots,\alpha B_{k}))$

That is, when $\mathcal{Y}=\operatorname{\mathbb{R}}$, then
$\mathcal{M}(\operatorname{\mathbb{R}})$, is the space of all probability
measures on $\operatorname{\mathbb{R}}$. If the sample space
$\mathcal{Y}=\operatorname{\mathbb{R}}$ is partitioned into measurable
subsets, then for every partition $(B_{1},B_{2},\ldots,B_{k})$, the prior
probability measure $\Pi$ on $(p(B_{1}),p(B_{2}),\ldots,p(B_{k}))$ is a
Dirichlet process prior.

The Dirichlet process $DP(G_{0},\alpha)$ is the simplest distribution - in
terms of the extent of independence it assumes - that can have this form. That
is probably why it has been referred to as _the normal distribution of
Bayesian nonparametrics_. More specifically, the location of the atoms, the
$\theta^{\prime}s$, are sampled $iid$ but the weights are sampled as
independent proportions. Now, the weights cannot be sampled $iid$ because then
they will not sum up to one. However, the next best thing to $iid$ sampling is
to obtain independent proportions, as is done using the stick-breaking
representation.

##  3 The Probabilistic Approach to Latent Variable Models

We begin by describing the probabilistic modeling approach. Next we explain
what a latent variable is and provide a general formulation of the latent
variable models. We end the section by giving a brief history of latent
variable modeling.

###  3.1 The Probabilistic Modeling Approach

A quick insight into the essential aspect of the probabilistic modeling
perspective can be gained by contrasting it with what _it is not_. The
probabilistic approach does not construct statistical models using algebraic
relations. In particular, reasoning about how observed and unobserved
quantities relate to each other and how their dependencies are induced does
not involve specifying structural relationships using algebraic formulas.
Instead, all observed and unobserved events that are thought to represent the
phenomenon of interest are denoted by random variables and each event is
allocated a finite probability of occurring. That is to say, the generative
probabilistic modeling approach is Bayesian in spirit (Pearl, 2014; Lauritzen,
1996; Jebara, 2012). In the probabilistic approach, we start with the joint
distribution of all observed and unobserved quantities that govern the
distribution of the data.

$p\left(\mathbf{y},\boldsymbol{\theta}\right)$

We construct a probabilistic model in two steps. In the first step, we
factorize the joint density into a product product of conditional and
marginals distributions. Through conditioning, we introduce dependencies that
transform the joint density into a mixture of densities. Implicit in this
mixture is a data-generating process that gives rise to the observations. As
Joe Blitzstein has often said, conditioning is the soul of statistics. Indeed,
it is the most important tool in a statistician’s toolkit for it can make
inference more tractable by allowing us to factorize any joint distribution
into the product of several conditional probability distributions. In the
second step, we specify particular distributions to the factorized densities.

The generative approach requires much more thought to apply in practice since
there are a large number of ways we can decompose the joint distribution into
a product of factors. Since each possible decomposition implicitly encodes a
particular number of modeling assumptions, the researcher is forced to examine
the validity of and the motivation for all modeling assumptions. Furthermore,
each component of the factorization needs to be assigned a probability
distribution that respects the plausibility of what is known about the
putative mechanism that gives rise to the data. For instance, in order to
treat the measurements as random in the first place, we would need to first
quantify the sources of uncertainty in inherent in the randomness using a
statistical model. Randomness can arise from either _experimental error_ ,
_model inadequacy_ , _parametric variability_ induced by random sampling from
a population, or even _parameter uncertainty_ arising from our inability to
exactly know the true values of the parameters (Kennedy and O’Hagan, 2001).

However, before we start making assumptions about the structure of the model
or the distributional form of its component factors, we need to examine two
important questions first; namely, how the data was collected, and what should
its grouping structure be. We start with the first issue.

####  3.1.1 Recording Measurements in Space and Time

In the probabilistic modeling framework, we begin the task of building a
statistical model, not with a particular model in mind, but rather with the
data. Consider once again the data matrix

$\displaystyle\mathbf{Y}=\begin{bmatrix}y_{1,1}&\dots&y_{1,P}\\\
\vdots&\ddots&\vdots\\\ y_{n,1}&\dots&y_{n,P}\\\ \vdots&\ddots&\vdots\\\
y_{N,1}&\dots&y_{N,P}\end{bmatrix}$ (9)

We approach the task of building a statistical model for the data by first
examining the individual measurement events that produce the recorded
observations. Since the event of taking a measurement occurs within a four-
dimensional spatiotemporal frame of reference, we denote a single measurement
not as $y_{np}$, but rather as $y(d_{1},d_{2},d_{3},t)$, a measurement
recorded at spatial location $\mathbf{d}=(d_{1},d_{2},d_{3})$ and time $t$,
where $\mathbf{d}\in D\subset\operatorname{\mathbb{R}}^{3}$ and $t\in
T\subset\operatorname{\mathbb{R}}$. Note that we do not consider here point
process models in which $D$ is random such that $\mathbf{d}$ denotes the
location of random events.

In many applications, spatiotemporal information is assumed to have a
miniminal effect on the observed data and is therefore omitted. For example,
if we wish to measure the heights of a number of individuals at a certain
point in time, it could be safe to assume that the spatial arrangement of the
subjects and the brief time interval that passes between the measurements are
inconsequential. If on the other hand, the height measurements for each
subject are taken multiple times over a period of months, then the time index
would need to be retained. To illustrate, if the height of individual $i$ is
recorded repeatedly over three occasions, we would then have the measurements
$[y(t^{*}_{1}),y(t^{*}_{2}),y(t^{*}_{3})]$. For uniform time intervals, the
notation can be simplified and the vector of repeated measures can be denoted
as $[y_{1},y_{2},y_{3}]$. In other applications, the location of measurements
is an essential aspect of understanding the phenomenon and cannot be ignored.

At any rate, the form in which the data is presented implicitly imposes
certain assumptions on us. For example, if the data come in the form of a data
matrix without an accompanying spatiotemporal index for each measured response
$y_{np}$, then we might incorrectly assume that the measurements were
uniformly sampled both in time or space. When the spatiotemporal indices are
absent such is the case for the data matrix $\mathbf{Y}$ given in Equation 9,
it would be good practice to consider the following few questions before
proceeding to model the data.

  * •

Are the $P$ measurements that make up a given observation $\mathbf{y}_{n}$
collected at the same time? If not, are the time lags between them
inconsequential?

  * •

Do the spatial indices $\mathbf{d}$ vary continuously over the three
dimensional frame of reference $D$ or is the space $D$ discretized?

  * •

Does the time index $t$ vary continuously over its domain $T$ or is the time
domain $T$ discretized?

  * •

In case of discrete time, are we dealing with multivariate time series (i.e.
the spatial coordinates are the same for each time period) or cross-sectional
data (i.e. the spatial coordinates vary over the time periods)?

In what follows, we will be assuming, when applicable, that both the temporal
domain $T$ and spatial domain $D$ are discretized uniformly; more
specifically, that $D$ is divided into a regular $P$ dimensional lattice where
$\mathbf{d}_{p}$ represents the spatial coordinates of the $p$th block and
that $T$ is divided into equal time intervals, where each measurement
represent an average over an time interval and a spatial block. We will also
be assuming that the spatial coordinates are fixed for all the time periods
and that there is no lag between the $P$ measurements. Consequently, the time
and space indices would no longer be needed since the row and column indices
of the data matrix are sufficient to preserve the spatiotemporal structure of
the measurements, in the sense that the complex 3D pattern of spatial
dependencies can still be inferred from the one dimensional response vector.

####  3.1.2 The Grouping Structure of the Data

In statistical studies, the observed data is usually characterized by a
grouping structure that determines what the observational units are. An
observation could be either a univariate random variable
$y_{n}\in\operatorname{\mathbb{R}}$, a multivariate random vector
$\mathbf{y}_{n}=(y_{n1},\ldots,y_{nP})\in\operatorname{\mathbb{R}}^{P}$, or
even a function $f_{n}\in L^{2}\subset\operatorname{\mathbb{R}}^{\infty}$,
where $f_{n}$ in this case is a point in the Lebesgue space $L^{2}$, an
infinite dimensional function space. Although the grouping structure is
usually determined by the experimental design or by how the data was
collected, the statistical model we specify for the data can either ignore,
relax, or even impose a particular grouping structure on the data. Our choice
of what constitutes an observational unit determines whether we approach the
modeling endeavor from a univariate, a multivariate, or a functional data
analysis perspective.

To illustrate, let us begin with the joint probability distribution of all
$K\left(=N\times P\right)$ variables that populate matrix $\mathbf{Y}$

$p\left(y_{11},\ldots,y_{np},\ldots,y_{NP},\boldsymbol{\theta}\right)$

If we map the two-subscript index set to a single -subscript index set, we can
express the joint as such

$p\left(y_{1},\ldots,y_{k},\ldots,y_{K},\boldsymbol{\theta}\right)$

Now if we refrain from making any distributional symmetry assumptions
regarding the $K$ random variables, we then have to treat all $K\ $
measurements as a _single_ high-dimensional sample from a joint distribution
$p_{\boldsymbol{\theta}}$ governed by parameter $\boldsymbol{\theta}$. That
is, since we do not assume that the joint distribution of the variables is
invariant to permutation or any other transformation of the indices, our data
would effectively consist of one single $K$ dimensional observation sampled
from a distribution $p_{\boldsymbol{\theta}}$

$\left(y_{1},\ldots,y_{k},\ldots,y_{K}\mid\boldsymbol{\theta}\right)\sim\quad
p_{\boldsymbol{\theta}}$

Such a model would be plausible if our single observation is a matrix-valued
random variable $\mathbf{Y}\in\operatorname{\mathbb{R}}^{N\times P}$ such that
of an image or a two-dimensional spatial data observation. In principle, the
joint model can account for any complex pattern of dependencies among the $K$
measurements irrespective of their row or column index. Nonetheless, the joint
distribution of a large number of random variables would require a large
number of parameters to enable it to capture the dependencies among all
possible combination of the variables. Furthermore, given the complex form of
dependencies that could be present in the data, statistical inference would be
a hopeless task with only one single sample. Note that in general the
parameter $\boldsymbol{\theta}$ can be a scalar, a multidimensional vector, or
even an infinite dimensional vector function. However, since for any scalar
observation, no more than one parameter can be learned, we are effectively
limited by the number of parameters we can learn given the $K$ measurements.
That is, for the parameter vector

$\boldsymbol{\theta}=\left(\theta_{1},\ldots,\theta_{d},\ldots,\theta_{D}\right)$

the dimensionality of the parameter $\boldsymbol{\theta}$ cannot exceed the
number of data points $K$ (i.e. $D\leqslant K$). An example of such a joint
model is the matrix normal model given by

$\mathrm{vec}(\mathbf{Y})\sim\mathcal{N}_{NP}(\mathrm{vec}(\mathbf{M}),\Sigma\otimes\Omega)$

where $\mathbf{M}$ is $N\times P$; $\Sigma$ is $N\times N$; $\Omega$ is
$P\times P$; and $\otimes$ denotes the Kronecker product. Accordingly, our
parameter vector of interest would consist of

$\boldsymbol{\theta}=[\mathrm{vec}(\mathbf{M)},\mathrm{vech}(\Sigma),\mathrm{vech}(\Omega)]$

with a dimension equal to $\frac{2NP+N(N+1)+P(P+1)}{2}$ \- which is always
bigger than $N\times P$. In general, note that when $P\gg N$, the number of
features is greater than the number of observations, the measurement density
$\alpha=\left(\frac{N}{P}\right)$ is finite. When $\alpha$ does not converge
to infinity as $N\to\infty$, we are faced with a high-dimensional inference
problem for which much of theory of classical statistics does not apply
(Advani and Ganguli, 2016).

Now consider the opposite scenario. If we have no grouping information about
the $K$ measurements then the most reasonable thing to do is to assume
exchangeablility and thus treat all the $K$ measurements as if drawn from a
common distribution $p_{\boldsymbol{\theta}}$

$y_{k}\overset{i.i.d.}{\sim}\quad p_{\boldsymbol{\theta}}$

Compared to the previous situation in which our data consisted of a single
$K$-dimensional vector observation, here we have a total of $K$ one-
dimensional scalar observations. With so many samples, inference should be
much more manageable. Furthermore, we now have much more flexibility in
choosing the dimension of $\boldsymbol{\theta}$. The dimensionality of the
parameter vector $\boldsymbol{\theta}$ controls the flexibility of the model’s
distribution. For instance, we can specify a low information capacity model
such as a univariate Gaussian model with unknown mean. In this case,
$\boldsymbol{\theta}$ would be one dimensional. At the other end of the
spectrum, we can specify a very flexible kernel density estimator to learn the
projected data’s density. Regardless, by collapsing high-dimensional data into
a one-dimension marginal distribution, we would still be unable to learn the
global patterns of dependencies that govern the data distribution in its
original high dimensional space.

In most statistical studies however, we would expect the data matrix to be
characterized by a grouping structure that preserves the distinction between a
feature and an observation, such that the rows of the matrix correspond to the
observations and the columns correspond to the $P$ features. The grouping
structure implicit in a data matrix effectively restricts how the indices can
be shuffled and transformed. Nonetheless, the data matrix can still be
characterized by complicated grouping structures. For example, the $P$
variables can be partitioned into a set of views (i.e. data sets) after which
the observations within a view can be partitioned into a set of clusters or
sub-populations (Mansinghka et al., 2016).

In what follows, we will treat the $N$ rows of the matrix as observations and
the $P$ columns as features such that the joint distribution is given by

$p\left(\mathbf{y}_{1},\mathbf{y}_{2},\ldots,\mathbf{y}_{N},\boldsymbol{\theta}\right)$

We can also partition the $P$-dimensional observation vector into $G$ multiple
groups, where $G\leq P$.

$\mathbf{y}_{n}=\left(\mathbf{y}_{n,1}^{\intercal},\ldots,\mathbf{y}_{n,g}^{\intercal},\ldots,\mathbf{y}_{n,G}^{\intercal}\right)$

Or even partition the $N$ observations into $C$ clusters such that

$\mathbf{y}=\big{(}\\{\mathbf{y}_{n}\\}_{n\in
N_{1}},\ldots,\\{\mathbf{y}_{n}\\}_{n\in N_{C}}\big{)}$

To obtain more complicated grouping structures, we can cluster the
observations after we group the features such that observations for group $g$
would belong to cluster $C_{g}$. The data can then be grouped as such.

$\mathbf{y}=\\{\mathbf{y}_{(n,g)}\\}_{n\in N_{C_{g}}}$

By partitioning the data matrix into groups and assuming partial
exchangeability, we can build a hierarchical model for all the groups. In a
hierarchical model, the measurements in each group are exchangeable within a
corresponding sub-model and the parameters of the sub-models are exchangeable
in their prior distribution (Gelman et al., 2014).

###  3.2 What is a Latent Variable?

Many definitions for a latent variable have been proposed in the literature
(see Bollen, 2002, for a review). Although the definitions that have been put
forward range from the informal to the mathematical, many of them rely on
distinctions that can be conceptually problematic when working in a
probabilistic modeling framework. For example, defining a variable as latent
based on whether it is treated as random rather than fixed does not carry over
to the probabilistic framework where all variables, known or unknown, are
treated as random. Furthermore, definitions that are based on the notion of
unobservability introduce some conceptual issues that complicate any attempt
to give a simple definition based on that concept alone. In a probabilistic
model, a variable is either observable or unobservable, but not all
unobservable variables in a probabilistic model correspond to what we think of
as a latent variable when we are working outside the probabilistic framework.
In particular, definitions based on unobservability break down when
considering how to treat the random noise term $\boldsymbol{\epsilon}_{n}$,
whether it is found in a simple regression model or in an implicit model that
transforms the distribution of the error into a more complicated distribution
$g(\boldsymbol{\epsilon}_{n};\boldsymbol{\theta})$ (Goodfellow et al., 2014).
In both cases, the error term is an unobservable random quantity whose
distribution deterministically induces a distribution on the response
$\mathbf{y}_{n}$ through an injective function $g$ whose inverse is $g^{-1}$.
Nonetheless, it cannot be considered a latent variable since the posterior

$\displaystyle
p(\boldsymbol{\epsilon}_{n}|\mathbf{y}_{n})=\mathbb{I}[\boldsymbol{\epsilon}_{n}=g^{-1}(\mathbf{y}_{n})]$
(10)

is a point mass, a deterministic not a random function (Tran et al., 2017).
Note that for a function to have an inverse, the dimension of
$\boldsymbol{\epsilon}_{n}$ should equal the dimension of $\mathbf{y}_{n}$. It
could be argued that for variable to be considered a latent variable it needs
be an unobserved random variable whose posterior distribution is not
degenerate. For example, in the Tobit model (Tobin, 1958)

$\displaystyle\begin{split}\epsilon_{i}\sim N(0,\sigma^{2})\\\ y_{i}^{*}=\beta
x_{i}+\epsilon_{i}\\\ y_{i}=max(0,y_{i}^{*})\end{split}$ (11)

the error term $\epsilon_{i}$ accordingly would not be considered a latent
variable since the function $g(x)=m+sx$ is invertible given that both
$\epsilon_{i}$ and $y_{i}^{*}$ are scalars (i.e. both of dimension one). In
contrast, $y_{i}^{*}$ can be considered a latent variable since its posterior
is not a point mass and the function $g(x)=max(0,x)$ is not invertible, even
though both $y_{i}^{*}$ and the observed variable $y_{i}$ are scalars. To
complicate things further, one could argue that in the Tobit model, we are in
fact dealing with a partially observable variable model $y_{i}=max(0,\beta
x_{i}+\epsilon_{i})$ rather than with a latent variable model per se.

Requiring a latent variable to be an unobservable random variable with a non-
degenerate posterior distribution can provide us with an adequate working
definition of a latent variable. Still we argue that a latent variable can be
better understood by the function it serves in a model. In particular, we
propose that the inclusion of a latent variable in a statistical model allows
us to accomplish two tasks: First, to capture statistical dependencies in the
data and second, to learn a latent structure that potentially governs the data
generating mechanism for the observations (i.e. learn the posterior
distribution $p(\mathbf{z}|\mathbf{y})$). With this in mind, note that the
notion of a latent variable as widely understood in the majority of modeling
context is most intuitively explained by the local independence definition of
a latent variable (Lord, 1952). We qualify the preceding statement by the word
majority because in the local independence definition does not always hold for
non-probabilistic models such as SEM where the latent variables can regressed
on each other to induce dependencies among the latent variable. Furthermore,
the local independence definition might not be as intuitive in the case of
models, whether estimable or not, in which there are a great number of latent
variables that combine together to give rise to an observed response of a much
lower dimension (Thomson, 1916).

To illustrate the basic idea, consider the following toy example. We observe
$N$ multivariate observations, $N$ samples of three dependent variables
$\\{(y_{n1},y_{n2},y_{n3})\\}_{n=1}^{N}$ (Figure 2(a)). The inclusion of a
scalar latent variable $z_{n}$ for the $n$th observation allows us to
factorize the joint distribution of ($y_{n1},y_{n2},y_{n3}$) into a set of
conditionally independent factors, as shown by the directed graph in Figure
2(b). For example, if we have a set of $P$ dependent variables where the
dependency is restricted to positive linear association among them, then
conditioning on just one single latent variable can be sufficient to explain
the pattern of dependencies in the $P$ dimensional data. The resulting latent
variable model for the $N$ observations can be concisely expressed by the
directed graph in Figure 2(c).

$y_{n1}$$y_{n2}$$y_{n3}$$n=1,\cdots,N$ (a) marginal model

$z_{n}$$y_{n1}$$y_{n2}$$y_{n3}$$n=1,\cdots,N$ (b) conditional model

$z_{n}$$y_{np}$$P$$N$ (c) Finite-data model

Figure 2: (a) Marginal graphical model for the three dependent variables. (b)
Directed graph of the observed variables conditional on the latent variable.
(c) The latent variable model as a directed graphical model for a sample of
$N$ multivariate observations, each a $P$-dimensional vector.

###  3.3 What is a Probabilistic Latent Variable Model?

Let us define probabilistic latent variable model as a statistical model of
observable variables that incorporates unobservable variables. Accordingly,
any parametric or nonparametric hierarchical Bayesian models subsumed under
the decomposition given by Equations 6 and 7 can be considered as special
cases in the general class of probabilistic latent variable models. To
formulate a general modeling framework for the family of probabilistic latent
variable models, we incorporate a D-dimensional local latent variable vector
$\mathbf{z}_{n}=(z_{1},\ldots,z_{D})$ for each observation $\mathbf{z}_{n}$
and a M-dimensional global latent variable vector
$\boldsymbol{\theta}=(\theta_{1},\ldots,\theta_{M})$ shared across all
observations (Blei et al., 2016). The local latent variable, $\mathbf{z}_{n}$,
serves the function of capturing local dependencies within the observation
vector $\mathbf{y}_{n}$ (e.g. clusters or common factors) whereas the global
latent variable $\boldsymbol{\theta}$ captures the dependencies across all, or
even a subset, of the $N$ observable random vectors. Accordingly, the joint
probability distribution for a probabilistic latent variable model can be
expressed as such

$\displaystyle
p(\mathbf{y}_{1},\cdots,\mathbf{y}_{N},\mathbf{z}_{1},\cdots,\mathbf{z}_{N},\boldsymbol{\theta})$
(12)

Now, there are countless many ways we could factor the joint distribution
given by Equation 12. Any chosen factorization encodes a set of implicit
assumptions that make the particular decomposition possible. For example, if
the observations in
$\mathbf{y}=\\{\mathbf{y}_{1},\mathbf{y}_{2},\cdots,\mathbf{y}_{N}\\}$ are
sequential in nature such that their distribution is no longer invariant under
permutation of indices (i.e. exchangeable), but is still invariant under
shifts (i.e. stationary), then a very general class of temporal generative
models (Gemici et al., 2017) can be formulated by the following decomposition.

$\displaystyle p(\mathbf{y}_{\ \leq N},\mathbf{z}_{\ \leq
N},\boldsymbol{\theta})=p(\boldsymbol{\theta})\prod_{n=1}^{N}\
p(\mathbf{y}_{n}\mid f_{x}(\mathbf{z}_{\ \leq n},\mathbf{y}_{\
<n}),\boldsymbol{\theta})p\big{(}\mathbf{z}_{n}\mid f_{z}(\mathbf{z}_{\
<n},\mathbf{y}_{\ <n}),\boldsymbol{\theta}\big{)}$ (13)

The family of generative temporal models give by Equation 13 include hidden
Markov models, non-linear state space models (Tornio et al., ), and the
Variational Recurrent Neural Network (VRNN) (Chung et al., 2015). If we are
not modeling temporal or sequential data, then we can assume that the N
multivariate observations are indeed exchangeable and obtain the following
decomposition instead.

$\displaystyle
p(\mathbf{y}_{1},\cdots,\mathbf{y}_{N},\mathbf{z}_{1},\cdots,\mathbf{z}_{N},\boldsymbol{\theta})=p(\boldsymbol{\theta})\prod_{n=1}^{N}\
p(\mathbf{y}_{n}\mid\mathbf{z}_{n},\boldsymbol{\theta})p(\mathbf{z}_{n}\mid\boldsymbol{\theta})$
(14)

This decomposition subsumes a very broad class of latent variable models such
as hierarchical Bayesian regression models (Gelman et al., 2014), Bayesian
mixture models, additive matrix decomposition models (e.g. factor analysis,
PCA, and CCA), and Bayesian nonparametric models. A graphical representation
of this general latent variable family of models is shown in Figure 3.

$\boldsymbol{\theta}$$z_{n}$$y_{np}$$P$$N$ (a) Unidimensional Latent Variable
Model

$\boldsymbol{\theta}$$\mathbf{z}_{n}$$\mathbf{y}_{n}$$N$ (b) Multi-dimensional
Latent Variable Model

Figure 3: (a) For each sample observation $n$, $z_{n}$ is a local latent
variable, and the $y_{np^{\prime}}$ are $P$ observed variables.
$\boldsymbol{\theta}$ is a vector of global latent variables (parameters) (b)
Directed graph of the observed variables conditional on the latent variable.

Notice that the decomposition given by Equation 14 implies a second
assumption; specifically, that the $N$ local latent variables $\mathbf{z}_{n}$
are exchangeable in their joint distribution. The joint distribution here
refers to the prior $p(\mathbf{z}_{n}\mid\boldsymbol{\theta})$ which is
governed by the global latent variable $\boldsymbol{\theta}$. The role of
global latent variable in the model is to control the extent of information
pooling _across_ the observations. To illustrate, if we set the variance of
the prior distribution to infinity, we obtain a separate model for each group
as a result. That is, a different model is fit to each observation vector,
where $\mathbf{y}_{n}$ can be thought of as an “experiment” with $P$
measurements. On the other extreme, setting the variance to zero would result
in complete pooling of information, effectively giving us a single model for
all $N\times P$ collapsed observations (Gelman et al., 2014). Incidentally, in
the context of mixed effects models, it is interesting to point out that a
fixed effect can be thought of as a special case of a random effect, but one
with an infinite variance prior.

###  3.4 Implicit Vs Explicit Models

Most of the time, we propose models whose likelihoods and priors have an
explicit closed-form probability densities. In some occasions, however, if we
are able to sample from a simulator that is better at modeling the generative
mechanism that gives rise to the data, then we can refrain from specifying an
explicit probability density for either the data distribution or the prior.
Models that do not specify explicit densities are known in the literature as
implicit models, first introduced by Diggle and Gratton (1984) to approximate
an intractable likelihood and later extended to approximating the posterior
distribution by Tavaré et al. (1997) in a method that is now known as
Approximate Bayesian Computation (ABC) Pritchard et al. (1999). In an implicit
model, instead of defining a likelihood, we define a deterministic function
$f(.)$ whose input is a random noise source $\boldsymbol{\epsilon}_{n}$ and
whose output is the observation

$\displaystyle\mathbf{y}_{n}=f(\boldsymbol{\epsilon}_{n}\mid\mathbf{z}_{n},\boldsymbol{\theta});\quad\boldsymbol{\epsilon}_{n}\sim\
s(\textperiodcentered)$ (15)

The function $f(.)$ can for example be implemented as a deep neural network.
Furthermore, the density for the local unobserved variable can also be defined
implicitly as such

$\displaystyle\mathbf{z}_{n}=f(\boldsymbol{\epsilon}_{n}\mid\boldsymbol{\theta});\quad\boldsymbol{\epsilon}_{n}\sim\
s(\textperiodcentered)$ (16)

We can sample from the posterior
$q(\mathbf{z}_{n}\mid\mathbf{y}_{n};\boldsymbol{\nu})$ for observation $n$

$\displaystyle\mathbf{z}_{n}\sim\mathcal{N}_{d}\big{(}\mathbf{z}_{n}|\mathbf{NN}_{\mu}(\mathbf{y}_{n};\boldsymbol{\nu}),\mathbf{NN}_{\sigma}(\mathbf{y}_{n};\boldsymbol{\nu})\big{)}$
(17)

by using the implicit function $g$

$\displaystyle\mathbf{z}_{n}=g(\epsilon_{n},\mathbf{y}_{n};\boldsymbol{\nu})=\mathbf{NN}_{\mu}(\mathbf{y}_{n};\boldsymbol{\nu})+\mathbf{NN}_{\sigma}(\mathbf{y}_{n};\boldsymbol{\nu})\odot\epsilon_{n}$
(18)

where $\odot$ denotes the elementwise product and

$\displaystyle\epsilon_{n}\sim\mathcal{N}_{d}\big{(}0,\mathbf{I}\big{)}$ (19)

If the neural networks $\mathbf{NN}$ are not invertible, as is usually the
case, then the function $g$ is not invertible and the likelihood is therefore
intractable (i.e. implicit). Almost all latent variable models can be
considered implicit models in the sense that their marginal distributions do
not have an explicit form. The most notable exception would be a factor
analysis model for which we have an explicit form for the marginal
distribution after the latent variables get integrated out.

###  3.5 Deep Vs Shallow Models

Also note that the decomposition given by Equation 14 above presumes that the
elements of the local latent variable vector $\mathbf{z}_{n}$ form a single
layer. _Deep models_ in contrast contain several hidden layers that form a
rich hierarchical structure. To obtain a joint factorization that includes
deep models as a general case, we just need to divide the vector of local
unobserved variables $\mathbf{z}_{n}$ into $L$ subsets
$\mathbf{z}_{n}=\\{\mathbf{z}_{n,1},\,\cdots,\mathbf{z}_{n,L}\\}$, such that
each $\mathbf{z}_{n,1}$ corresponds to a separate level in the sampling
hierarchy. Unlike shallow models, deep models with a hierarchy of latent
variable dependencies are better able to capture rich latent structure
(Ranganath et al., 2015).

##  4 Brief History of Latent Variable Models

Although latent variable modeling is a branch of multivariate statistics, many
of the field’s major developments originated and were motivated by research
problems, not within statistics, but rather from the psychological and social
sciences (Izenman, 2008). Whereas the early breakthroughs mainly came from
psychology, educational measurement, and economics, the more recent advances
are increasingly coming from the field of machine learning, especially from
the very active research area of unsupervised learning in which researchers
are trying to accomplish the goal of endowing computers the ability to learn
complex data patterns in an unsupervised automatic manner (Lake et al., 2015;
Rezende et al., 2016).

The quest to understand learning and intelligence, whether in humans or
machines, has always been a motivating force in the field from its early
beginnings. Indeed, the psychologist Charles Spearman, influenced by Francis
Galton’s development of the correlation coefficient, introduced _factor
analysis_ in 1904 as a statistical method to study the organization of mental
abilities in humans (Spearman, 1904). He proposed that from the correlations
of multiple observed variables (e.g. tests of general ability), we can infer a
single latent factor that can explain the common variance shared among all the
observations. He referred to the latent variable common to all the measured
ability tests as _g_ , the general factor of intelligence.

###  4.1 Summary of Key Developments

The factor analytic approach paved the way for, what seemed to many in the
social and psychological sciences, an objective and quantitative framework for
measuring unobservable _hypothetical constructs_ such as “self-esteem”,
“verbal ability”, and “quality of life”. Indeed Spearman’s single factor
model, in which both the one latent variable and the observed variables are
assumed continuous, has laid the foundation for a long research program in
psychometrics that has spanned well over a century. To account for the
inadequacy of a single factor to account for the covariance pattern in testing
data, the psychologist Thurstone extended the single factor model to multiple
factors (Thurstone and L.L., 1947). The factor model was constantly refined
and generalized, culminating in the _structural equation model_ introduced by
(Jöreskog, 1970) in which linear structural relationships between the latent
variables are also modeled. In educational testing, _item response theory_
models (Birnbaum, 1968) were developed to accommodate categorical observed
variables encountered in testing in which test items assume two values,
correct and incorrect. In sociology, clustering analysis motivated the
development of _latent class models_ (i.e. mixture models) in which the latent
variables are assumed to be categorical (Lazarsfeld, 1950). Since then, a
large body of work has been built on top of these foundational early advances
and many authors (Muthén, 2002; Bartholomew et al., 2011; Hwang and Takane,
2014) have attempted to bring together the various models under one roof in a
general framework. The most encompassing as of yet is the Generalized Linear
Latent and Mixed model (GLLMM) introduced by Rabe-Hesketh et al. (2004) as a
unified general framework that brings together not just factor analysis, IRT,
structural equations, and latent class models commonly encountered in
psychometrics and social quantitative methods, but also mainstream statistical
random effects models such as multilevel regression, and longitudinal models.

###  4.2 The Statistical Myth of Measuring Hypothetical Constructs

The enthusiasm regarding the factor analytic approach to uncover underlying
theoretical constructs was not universally shared. In particular, the
interpretation of a latent variable as a hypothetical construct was much
criticized. One of the earliest criticism came from Spearman’s colleague, Karl
Pearson, widely regarded as the founder of modern mathematical statistics. The
two men constantly feuded with each other throughout their careers due to
their opposing scientific philosophies (Hägglund, 2001). Whereas Spearman
strove to discover the fundamental laws of psychology and invent a method by
which he can objectively measure latent traits such as intelligence, Pearson
did not see a place for unobserved and subsequently unmeasurable phenomena in
the scientific enterprise. Incidentally, Pearson (1901) had pioneered a method
closely related to, and often conflated with Factor Analysis (FA) - _principle
component analysis_ (PCA). Pearson’s PCA, which was later formulated as
multivariate technique by Hotelling (1933), was introduced not as an
explanatory model, but as a dimensionality reduction technique in which the
data is geometrically transformed (i.e. projected) into a lower dimensional
subspace.

The criticisms nonetheless did not end with Pearson. The English educational
psychologist Thomson (1916) proposed a sampling model that generated the same
correlation structure as a general factor model but that did not assume a
single factor. In particular, whereas Spearman’s model, and later extensions,
proposed that positive correlations among P variables suggest the existence of
one or more factors fewer than P, Thomson’s sampling model assumed the
opposite: that there are more factors than we have variables and that each
variable is the result of a great number of independent causes, much greater
than P. Even though it seems more plausible to assume that there are a great
number of factors (mental subprocesses) that contribute to a score on any one
ability test or item questionnaire, psychometricians following Spearman’s
footsteps ignored these criticisms and pursued a research program based on the
indefensible premise that factor analysis and its extensions are statistical
methods that can explain correlations in the observed data, rather than just
describe them. It can be argued that the explanatory approach in which the
existence of a hypothetical construct is inferred by performing exploratory
analysis on limited sample self-report data cannot substitute the rigorous
process of developing operational definitions of unobservable constructs on
theoretical grounds or even on experimental measurements of the cognitive
processes underlying the observed responses. Indeed, throughout the history of
psychometrics, there have been repeated calls to guard against the
_reification fallacy_ in which a latent variable (e.g. mental ability) that
captures a common variability in the data is treated as if it were a physical
quantity (e.g. height) that is amenable to direct measurement (Gould, 1996;
Sobel, 1994).

Moreover, given that the definition of a statistical model provided in section
2 does not say much about the nature of the relationship between a statistical
model and its target phenomenon observed in the physical universe, we would
like to point out here that the modeling perspective we adopt is founded on
the premise that a statistical model is a _phenomenological model_ of reality
describing relationships among observed variables in a mathematical formalism
that is not derived from first principles (Hilborn and Mangel, 1997). We
consider a model phenomenological in the Kantian sense (Kant and Guyer, 1998)
in which the word _phenomenon_ denotes an empirical object that is known
through the senses and whose mental representation is constructed by
experience. The _phenomenon_ is contrasted with the _noumenon_ , the “thing-
in-itself” that is not directly knowable by the senses but which can be known
by the mind in the form of a theoretical construct. For example, whereas the
electromagnetic force is a theoretical construct arrived at from first
principles and provides a causal explanation to empirical observations (e.g.
effects of magnetic and electric fields), the empirical observations
themselves constitute experienced phenomenon whose correlative nature can
provide nothing more than evidence for noumenal knowledge.

###  4.3 The First Latent Variable Model

If we restrict latent variable models to those models in which latent
variables denote hypothetical constructs, arguably Spearman’s factor model
could be considered the first formulation of a latent variable model. However,
such a restriction based on a narrowly defined interpretation is not
universally accepted neither within psychometrics, the field Spearmen helped
start, nor within the larger statistics community (Skrondal and Rabe-Hesketh,
2004). If we just consider the statistical form of the model, then it could be
argued that the first latent variable model should be attributed to the
astronomer Airy, who in 1861 proposed a variance components model (i.e. random
effects model) for the repeated measurement of an astronomical object over
several nights (Airy, 1861). Both Airy’s variance components model and
Spearman’s factor model have the same graphical representation depicted in
Figure 2(c). The main difference lies in that the latent variable in Airy’s
model is a random effect that captures unobserved heterogeneity whereas in
Spearman’s model the latent variable denotes a hypothetical construct which we
seek to measure.

In his model, Airy makes $P$ repeated measurements of an astronomical
phenomenon for each of $N$ nights. See Searle et al. (2009) for a more details
and the original formulation. Assuming a balanced design, the model can be
expressed as a probabilistic generative model by the following two-stage
sampling hierarchy.

$\displaystyle\begin{split}\mathbf{y}_{n}\mid z_{n}&\sim\
\mathcal{N}_{P}(\mathbf{1}\mu+\mathbf{1}z_{n},\,\sigma^{2}_{\epsilon}\mathbf{I})\\\
z_{n}&\sim\mathcal{N}(0,\sigma^{2}_{z})\quad\textit{for
}n=1,\cdots,N\end{split}$ (20)

where $\mathbf{y}_{n}\in\mathbb{R}^{P}$ is the vector of repeated measurements
and $z_{n}\in\mathbb{R}$ is a latent variable (i.e. random effect).

###  4.4 Latent Variable Models in Machine Learning

Although the GLLMM framework does indeed treat random coefficient models
(which include hierarchical and longitudinal models) as latent variable
models, the framework leaves out important multivariate data reduction methods
such PCA, Independent Component Analysis (ICA), Canonical Correlation Analysis
(CCA) from the formulation. The argument put forward by Skrondal and Rabe-
Hesketh (2004) contends that PCA and similar procedures are not latent
variable models because they cannot be formulated as statistical models with
latent variables corresponding to lower dimensional space. However, their
assertion does not seem to have taken into consideration results that had been
previously published in the literature. More specifically, in recent years,
research in latent variable modeling has undergone a very active resurgence in
the field of unsupervised machine learning. Indeed, great progress has been
made in developing new latent variable models and recasting existing
clustering and dimensionality reduction methods in a probabilistic graphical
model formulation. For example, building on previous research in neural
networks (MacKay, 1995a; Bishop et al., 1998), Tipping and Bishop (1999)
demonstrated that PCA can be indeed be formulated as a probabilistic latent
variable model. In their paper, they showed that Probabilistic PCA is a
special case of Factor Analysis (FA). More specifically, whereas in FA the
error covariance matrix is diagonal, in Probabilistic PCA, it is
$\sigma\mathbf{I}$ and standard PCA is obtained when $\sigma\rightarrow 0$.
Roweis and Ghahramani (1999) review several common methods showing that not
just PCA and FA, but also ICA, Kalman Filters, and Hidden Markov Models (HMMs)
can be considered as special cases of a general linear-Gaussian latent
variable model. A probabilistic formulation of CCA (Bach and Jordan, 2005)
later followed. Recently, a Group Factor Analysis (GFA) model has been
proposed by Klami et al. (2015) that generalizes Probabilistic CCA to more
than two data sets.

It is these developments that motivate the material in the next two sections.
In the social sciences, latent variable models have been traditionally been
viewed and used as tools to study hypothetical constructs. However, ongoing
research in machine learning, deep learning, and probabilistic graphical
modeling has given a new life and meaning to latent variable modeling. These
recent developments present latent variables as capable of endowing
statistical models not only with the ability to perform essentials data
analysis tasks such as dimensionality reduction, statistical clustering, and
accounting for unobserved heterogeneity, but also with explanatory and
interpretative powers. The inclusion of latent variables allows us to capture
statistical dependencies among observed variables and thus learn the hidden
structure underlying the data generating mechanism. This is especially true
when latent variable models are cast in a generative probabilistic framework,
a perspective that the next section attempts to describe in detail. The
following sections forms the foundation for subsequent sections where many of
the most commonly encountered latent variable models are constructed from the
most basic elements.

##  5 Constructing Probabilistic Latent Variable Models

In this section, we will attempt to show how starting from a joint probability
density, we can formulate a spatiotemporal models, multivariate models,
regression models, and latent variable models. We will also attempt to
describe how regression models are related to latent variable models through
reduced rank regression. We end the section by giving an overview of latent
variable models most commonly encountered in practice.

To begin, let us consider a two-group partition of the observation vectors
that populate the data matrix $\mathbf{Y}$. To simply notation, we can
effectively create a two-group partition by concatenating the data matrix with
a matrix of co-variates such that each row would consist of a P-dimensional
observation vector $\mathbf{y}_{n}=(y_{n1},\cdots,y_{nP})$ and an
M-dimensional covariate vector $\mathbf{x}_{n}=(x_{n1},\cdots,x_{nM})$. The
concatenated data would be an $N$ by $(P+M)$ matrix consisting of a total of
$N\times(P+M)$ measurements

$\displaystyle\begin{bmatrix}(\mathbf{y}_{1}^{\rm T},\mathbf{x}_{1}^{\rm
T})\\\ \vdots\\\ (\mathbf{y}_{n}^{\rm T},\mathbf{x}_{n}^{\rm T})\\\ \vdots\\\
(\mathbf{y}_{N}^{\rm T},\mathbf{x}_{N}^{\rm
T})\end{bmatrix}=\begin{bmatrix}y_{11}&\dots&y_{1P}&x_{11}&\dots&x_{1M}\\\
\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\
y_{i1}&\dots&y_{nP}&x_{n1}&\dots&x_{nM}\\\
\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\
y_{N1}&\dots&y_{NP}&x_{N1}&\dots&x_{NM}\end{bmatrix}$ (21)

with a joint distribution given by

$p\big{(}\left(\mathbf{y}_{1},\mathbf{x}_{1}\right),\cdots,\left(\mathbf{y}_{N},\mathbf{x}_{N}\right),\boldsymbol{\theta}\big{)}$

###  5.1 Spatiotemporal Models

We can treat the data as spatiotemporal in nature by letting the row index
denote time and mapping the column index to a given spatial coordinate. Using
a spatio-temporal mode, we can learn both the temporal and spatial
dependencies of the data. In a spatiotemporal model we do not treat the
observations that populate the data matrix as exchangeable neither in the
temporal nor in the spatial dimension. In the spatial dimension, the patterns
of dependencies tend to vary in complexity depending on the type of data. In
the temporal dimension in contrast, the observations tend to be characterized
by a less complex pattern of sequential dependencies. For example, the
observation $\mathbf{y}_{n}$ can be dependent on the preceding observation
$\mathbf{y}_{n-1}$ only.

If we wanted a very general model, we can formulate a non-separable spatio-
temporal model in which the spatial and temporal components are allowed to
interact. If we assume some form of stationarity in the temporal process such
that the joint probability distribution of a finite sequence of consecutive
observations does not change when they shifted in time, then a family of non-
separable spatio-temporal models can be specified by the following
decomposition.

$\displaystyle p(\mathbf{y}_{\ \leq N};\boldsymbol{\theta})=\prod_{n=1}^{N}\
p(\mathbf{y}_{n}\mid f(\mathbf{y}_{\ <n});\boldsymbol{\theta})$ (22)

where we let the subscript $n$ index time such that the response
$\mathbf{y}_{n}$ can be dependent on preceding responses $\mathbf{y}_{\ <n}$
through a nonlinear function $f$.

Such a model is commonly encountered and is appropriate for modeling dynamical
systems in which the current spatial component of the process evolves as a
function of past spatial locations. Although non-separable models tend to be
realistic for many real world applications, their estimation can be
computationally intractable, especially in the case when the dimensionality of
N or P is high.

Alternatively, when adequate motivation exists, a separable spatio-temporal
model can be used instead. We can formulate a separable model by assuming that
the parameters that govern the distribution of the temporal component are
independent of those that govern the spatial component. To obtain a separable
model, we let
$\boldsymbol{\theta}=(\boldsymbol{\theta}^{S},\boldsymbol{\theta}^{T})$ such
that $\boldsymbol{\theta}^{S}$ governs the distribution of the $P$ elements of
$\mathbf{y}_{n}$ whereas $\boldsymbol{\theta}^{T}$ governs the distribution of
the $N$ observations.

$\displaystyle p(\mathbf{y}_{\ \leq N};\
\boldsymbol{\theta}^{S},\boldsymbol{\theta}^{T})=\prod_{n=1}^{N}\
p(\mathbf{y}_{n}\mid f_{\boldsymbol{\theta}^{T}}(\mathbf{y}_{\ <n}),\
\boldsymbol{\theta}^{S})$ (23)

###  5.2 Multivariate Models

A spatiotemporal model can be viewed either as a temporally varying spatial
model or as spatially varying time series model (Banerjee et al., 2014). In
either case, the data modeled is assumed to be characterized by both spatial
and temporal dependencies. In contrast, in multivariate models, we relax the
stationarity assumption that allows us to model the sequential dependencies in
the temporal dimension. Instead, we assume that the observations are
exchangeable in their joint distribution. Exchangability give us the following
factorization.

$p\big{(}\left(\mathbf{y}_{1},\mathbf{x}_{1}\right),\cdots,\left(\mathbf{y}_{N},\mathbf{x}_{N}\right),\boldsymbol{\theta}\big{)}=\prod_{n=1}^{N}p\big{(}\left(\mathbf{y}_{n},\mathbf{x}_{n}\right)\mid\boldsymbol{\theta}\big{)}p(\boldsymbol{\theta})$

Note that the parameter vector $\boldsymbol{\theta}$ governs the distribution
of the pairs of $\mathbf{y}_{n}$ and $\mathbf{x}_{n}$ jointly. However if we
wish to predict $\mathbf{y}_{n}$ from $\mathbf{x}_{n}$, we need to assume that
the parameters that govern the distribution of $\mathbf{y}_{n}$ are
independent of those that govern the distribution of $\mathbf{x}_{n}$. That
is, given the partition
$\boldsymbol{\theta}=(\boldsymbol{\theta}^{Y},\boldsymbol{\theta}^{X})$, we
have

$p(\boldsymbol{\theta}^{Y},\boldsymbol{\theta}^{X})=p(\boldsymbol{\theta}^{Y})p(\boldsymbol{\theta}^{X})$

and we can factorize the joint distribution as such.

$\displaystyle
p\big{(}\left(\mathbf{y}_{1},\mathbf{x}_{1}\right),\cdots,\left(\mathbf{y}_{N},\mathbf{x}_{N}\right),\boldsymbol{\theta}\big{)}=\prod_{n=1}^{N}p(\mathbf{y}_{n},\mid
f_{\theta}(\mathbf{x}_{n}),\boldsymbol{\theta}^{Y})\prod_{n=1}^{N}p(\mathbf{x}_{n}\mid\boldsymbol{\theta}^{X})p(\boldsymbol{\theta}^{Y})p(\boldsymbol{\theta}^{X})$
(24)

The factorization implies that the distribution of $\mathbf{x}_{n}$ does not
influence the conditional distribution of $\mathbf{y}_{n}$. The assumption
allows us to ignore the distribution of $\mathbf{x}_{n}$ when the focus is on
determining the parameter $\theta^{Y}$ \- such as the case in many regression
or path analysis models for example, where $\mathbf{x}_{n}$ is considered an
exogenous variable (Wright, 1921).

####  5.2.1 Normal Linear Models

Let us split the vector of parameters into two components such that

$\theta^{Y}=[vec(\mathbf{B}),\ vech(\boldsymbol{\Sigma})]$

where

$\mathbf{B}=\begin{bmatrix}\boldsymbol{\beta}_{1}^{\rm T}\\\
\boldsymbol{\beta}_{2}^{\rm T}\\\ \vdots\\\ \boldsymbol{\beta}_{P}^{\rm
T}\end{bmatrix}=\begin{bmatrix}\beta_{11}&\cdots&\beta_{1m}\\\
\beta_{21}&\cdots&\beta_{2m}\\\ \vdots&\ddots&\vdots\\\
\beta_{p1}&\cdots&\beta_{pm}\end{bmatrix}$

and $\boldsymbol{\Sigma}$ is positive definite matrix, and let us further
assume that

$\displaystyle\mu(\mathbf{x}_{n};\mathbf{B})=\mathbf{B}\mathbf{x}_{n}\quad\textit{and}\quad\boldsymbol{\sigma}(\boldsymbol{\Sigma})=\boldsymbol{\Sigma}$

we arrive at following data distribution for a linear regression model

$\displaystyle
p\left(\mathbf{y}_{n}\mid\mathbf{x}_{n},\theta^{Y}\right)\equiv\mathcal{N}_{p}\left(\mathbf{y}_{n}\mid\mathbf{B}\mathbf{x}_{n},\
\boldsymbol{\Sigma}\right)$

Note that since both $\mathbf{y}_{n}$ and $\mathbf{x}_{n}$ are multivariate
vectors, $\mu(\mathbf{x}_{n};\mathbf{B})$ is vector-valued function that
transforms $\mathbf{x}_{n}$ into $\mathbf{y}_{n}$. Also note that we dropped
the index $n$ from both functions. The joint probability distribution can now
be factored as such.

$\displaystyle\prod_{n=1}^{N}\mathcal{N}_{p}(\mathbf{y}_{n}\mid\mathbf{B}\mathbf{x}_{n},\
\boldsymbol{\Sigma})p\left(\mathbf{B},vech(\boldsymbol{\Sigma})\right)\prod_{n=1}^{N}\
p(\mathbf{x}_{n}\mid\theta^{X})p(\theta^{X})$ (25)

Let us simplify things further and let
$\boldsymbol{\Sigma}=\operatorname{Diag}(\sigma^{2}_{1},\cdots,\sigma^{2}_{p})$
so that the elements of the response vector now become uncorrelated. Also, let
us group $\mathbf{B}$ as set of $P$ $M$-dimensional observations
$\left(\boldsymbol{\beta}_{1},\boldsymbol{\beta}_{2},\ldots,\boldsymbol{\beta}_{P}\right)$,
and let us consider the predictor $\mathbf{x}_{n}$ to be fixed, so that we
just focus on the data distribution for now. These assumptions give us the
following decomposition

$\displaystyle\prod_{n=1}^{N}\mathcal{N}_{p}(\mathbf{y}_{n}\mid\mathbf{B}\mathbf{x}_{n},\
\boldsymbol{\Sigma})p\left(\boldsymbol{\beta}_{1},\boldsymbol{\beta}_{2},\ldots,\boldsymbol{\beta}_{P},\sigma^{2}_{1},\cdots,\sigma^{2}_{P}\right)$

By grouping the pairs $\left(\boldsymbol{\beta}_{p},\sigma^{2}_{p}\right)$ and
then assuming conditional exchangeability of the parameter pairs, we can
further decompose the model into

$\displaystyle\prod_{n=1}^{N}\mathcal{N}_{p}(\mathbf{y}_{n}\mid\mathbf{B}\mathbf{x}_{n},\
\boldsymbol{\Sigma})=\prod_{n=1}^{N}\prod_{p=1}^{P}\mathcal{N}(y_{np}\mid\boldsymbol{\beta}^{\rm
T}_{p}\mathbf{x}_{n},\
\sigma^{2}_{p})p\left(\boldsymbol{\beta}_{p},\sigma^{2}_{p}\mid\xi\right)p(\xi)$

If we assume that the priors over the parameters
$\left(\boldsymbol{\beta}_{p},\sigma^{2}_{p}\right)$ to be uninformative, we
further reduce the decomposition to

$\displaystyle\prod_{n=1}^{N}\mathcal{N}_{p}(\mathbf{y}_{n}\mid\mathbf{B}\mathbf{x}_{n},\
\boldsymbol{\Sigma})=\prod_{n=1}^{N}\prod_{p=1}^{P}\mathcal{N}(y_{np}\mid\boldsymbol{\beta}^{\rm
T}_{p}\mathbf{x}_{n},\ \sigma^{2}_{p})$

The above simplification implies that a multivariate regression model can be
broken down into a set of $P$ multiple regression problems. For example, if we
only consider one response variable $y\in\mathbf{y}$, we can then drop the
subscript and obtain the following decomposition for a multiple linear
regression model.

$\displaystyle\prod_{n=1}^{N}\ p(y_{n}\mid\mathbf{x}_{n})\ $
$\displaystyle=\prod_{n=1}^{N}\mathcal{N}(y_{n}\mid\boldsymbol{\beta}^{\rm
T}\mathbf{x}_{n},\ \sigma^{2})$

Note that the decomposition above is not exactly that of a hierarchical
Bayesian regression model. In a hierarchical regression, the $N$ multivariate
observation vector pairs $(\mathbf{y}_{n},\mathbf{x}_{n})$ are partitioned
into $J$ clusters and for each cluster $j$ we introduce a regression
coefficient $\boldsymbol{\beta}_{jp}$. Accordingly, the decomposition can be
expressed as

$\displaystyle\prod_{n=1}^{N_{j}}\prod_{j=1}^{J}\mathcal{N}_{p}(\mathbf{y}_{n}\mid\mathbf{B}\mathbf{x}_{n},\
\boldsymbol{\Sigma})=\prod_{n=1}^{N_{j}}\prod_{j=1}^{J}\prod_{p=1}^{P}\mathcal{N}(y_{njp}\mid\boldsymbol{\beta}^{\rm
T}_{jp}\mathbf{x}_{nj},\
\sigma^{2}_{p})p\left(\boldsymbol{\beta}_{jp},\sigma^{2}_{p}\mid\xi\right)p(\xi)$

###  5.3 Reduced Rank Regression (RRR)

Our starting point is the decomposition given in Equation 25 for the
multivariate regression model - which is reproduced here.

$\displaystyle\prod_{n=1}^{N}\mathcal{N}_{p}(\mathbf{y}_{n}\mid\mathbf{B}\mathbf{x}_{n},\
\boldsymbol{\Sigma})p(\mathbf{B},\boldsymbol{\Sigma})\prod_{n=1}^{N}\
p(\mathbf{x}_{n}\mid\theta^{X})p(\theta^{X})$

and the parameters that govern $\mathbf{x}_{n}$ and $\mathbf{y}_{n}$ to be
independent, and

$\displaystyle\prod_{n=1}^{N}\mathcal{N}_{p}(\mathbf{y}_{n}\mid\mathbf{B}\mathbf{x}_{n},\
\boldsymbol{\Sigma})$

If we allow the rank of the regression coefficient matrix to be deficient such
that $rank(\mathbf{B})=d\ <\ min(m,p)$, then the matrix can be decomposed into
two non-unique matrices as such $\overset{p\times
m}{\mathbf{B}}=\overset{p\times d}{\mathbf{W}}\ \overset{d\times
m}{\mathbf{D}}$. The decomposition results in what is known as Reduced Rank
Regression (RRR) (Izenman, 1975)

$\displaystyle\prod_{n=1}^{N}\mathcal{N}_{p}(\mathbf{y}_{n}\mid\mathbf{W}\mathbf{D}\mathbf{x}_{n},\
\boldsymbol{\Sigma})$

Now if set the predictor vector to be equal to the response vector then the
resulting data distribution of the reduced rank regression can be expressed as

$\displaystyle\prod_{n=1}^{N}\mathcal{N}_{p}(\mathbf{y}_{n}\mid\mathbf{W}\mathbf{D}\mathbf{y}_{n},\
\boldsymbol{\Sigma})$

By letting ${\mathbf{z}_{n}}=\overset{d\times m}{\mathbf{D}}\overset{m\times
1}{\mathbf{y}_{n}}$, we can now think of $\mathbf{y}_{n}$ as a message;
${\mathbf{z}_{n}}$ as an encoded transformation of the message; and
$\mathbf{W}$ as the decoding transformation. The model can then be expressed
as

$\displaystyle\prod_{n=1}^{N}\mathcal{N}_{p}(\mathbf{y}_{n}\mid\mathbf{W}\mathbf{z}_{n},\
\boldsymbol{\Sigma})$

This representation allows us to see that several latent variable models are
special cases of $RRR$ (Izenman, 2008).

###  5.4 Principle Component Analysis (PCA)

To obtain our first and simplest latent variable model, we specify a
distribution to the latent vector since the lower dimensional encoded message
is unobservable. Assigning a standard Gaussian to the latent variables and
restricting the covariance of the errors to be isotropic gives us
probabilistic PCA (Tipping and Bishop, 1999). The model can be expressed as

$\displaystyle\prod_{n=1}^{N}\mathcal{N}_{p}(\mathbf{y}_{n}\mid\mathbf{W}\mathbf{z}_{n},\
\sigma^{2}\mathbf{I})\mathcal{N}_{d}(\mathbf{z}_{n}\mid\boldsymbol{0},\,\mathbf{I})$

or as a generative model

$\displaystyle\begin{split}\mathbf{y}_{n}\mid\mathbf{z}_{n}&\sim\
\mathcal{N}_{p}(\mathbf{W}\mathbf{z}_{n},\,\sigma^{2}\mathbf{I})\\\
\mathbf{z}_{n}\ &\sim\ \mathcal{N}_{d}(\boldsymbol{0},\,\mathbf{I})\
\quad\textit{for }n=1,\cdots,N\end{split}$

Standard PCA is obtained as the maximum likelihood solution if we let
$\sigma^{2}\rightarrow 0$ and constrain $\mathbf{W}$ to be orthogonal.

###  5.5 Factor Analysis (FA)

If we relax the constraint of homogeneous variance (i.e. $\sigma^{2}$ is the
same for all components of the vector $\mathbf{y}_{n}$), the covariance matrix
can be specified as

$\displaystyle\boldsymbol{\Sigma}=\operatorname{Diag}(\sigma^{2})=\begin{bmatrix}\sigma^{2}_{1}&0&\ldots&0\\\
0&\sigma^{2}_{2}&\ldots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&\ldots&0&\sigma^{2}_{p}\end{bmatrix}$

The result is a multidimensional factor analysis model with a diagonal
covariance and $d$ latent factors.

$\displaystyle\prod_{n=1}^{N}\mathcal{N}_{p}(\mathbf{y}_{n}\mid\mathbf{W}\mathbf{z}_{n},\
\operatorname{Diag}(\sigma^{2}))\mathcal{N}_{d}(\mathbf{z}_{n}\mid\boldsymbol{0},\,\mathbf{I})$

or as a generative model

$\displaystyle\begin{split}\mathbf{y}_{n}\mid\mathbf{z}_{n}&\sim\
\mathcal{N}_{p}(\mathbf{W}\mathbf{z}_{n},\,\operatorname{Diag}(\sigma^{2}))\\\
\mathbf{z}_{n}\ &\sim\ \mathcal{N}_{d}(\boldsymbol{0},\,\mathbf{I})\
\quad\textit{for }n=1,\cdots,N\end{split}$ (26)

Note that the standard normal prior over the latent variables resolves
invariance in the scale and location of the latent space, but not in rotation.

###  5.6 Independent Component Analysis (ICA)

Independent Component Analysis (Hyvärinen, 2015) is a widely used model for
signal source separation. In order to model signals as latent variables, their
distribution cannot be normal because the Gaussian is the most random of all
distributions. In fact, according to _Jayne’s principle of maximum entropy_
(Jaynes, 2003), a Gaussian random variable has maximum entropy
(i.e.uncertainty) among all random variables with equal variance. In contrast,
entropy tends to be small for spiked, or super-Gaussian distributions -
distributions that characterize signals! So if we consider the _generalized
Gaussian distribution_ as a prior over the latent variables, we obtain a
highly peaked supergaussian distribution (if the shape parameter $\alpha$ set
to be less than 2) that is appropriate for modeling signals that give rise to
sparse data. Note that the version of the generalized Gaussian distribution we
are considering is that of parametric symmetric distributions. The normal
distribution is a special case obtained when $\alpha=2$ and the Laplace
distribution is another special case when $\alpha=1$ (Gómez et al., 1998). ICA
can be thought of as a linear latent generative model, a type of a Factor
Analysis model but with nongaussian priors over the latent variables (i.e.
sources). The model has the following decomposition.

$\displaystyle\prod_{n=1}^{N}\mathcal{N}_{p}(\mathbf{y}_{n}\mid\mathbf{W}\mathbf{z}_{n},\
\operatorname{Diag}(\sigma^{2}))\prod_{n=1}^{N}\prod_{d=1}^{D}\mathcal{GG}_{d}(z_{d,n}\mid\alpha_{d})\quad\alpha_{d}<2$

or as a generative model

$\displaystyle\begin{split}\mathbf{y}_{n}\mid\mathbf{z}_{n}&\sim\
\mathcal{N}_{p}(\mathbf{W}\mathbf{z}_{n},\,\operatorname{Diag}(\sigma^{2}))\\\
z_{d,n}\ &\sim\ \mathcal{GG}_{d}(\alpha_{d})\quad\alpha_{d}<2\quad d=1:D;\
\quad n=1,\cdots,N\end{split}$

where $\mathcal{GG}$ is the generalized Gaussian distribution.

###  5.7 Canonical Correlation Analysis (CCA)

Hotelling (1936) introduced Canonical Correlation Analysis (CCA) as a
dimensionality reduction method that projects two sets of random variables
(i.e. two data sets) into a shared subspace such that the projections are
maximally correlated. Although almost always commonly thought as data
reduction technique, CCA can indeed be given a latent variable interpretation
as Bach and Jordan (2005) have shown. More specifically, they showed that the
equivalent probabilistic latent variable model can be formulated as

$\displaystyle\begin{split}\mathbf{z}_{n}&\sim\mathcal{N}_{d}\left(\boldsymbol{0},\mathbf{I}\right)\quad
n=1,\cdots,N\\\ \begin{bmatrix}{\mathbf{y}_{n}}^{(1)}\mid\mathbf{z}_{n}\\\
{\mathbf{y}_{n}}^{(2)}\mid\mathbf{z}_{n}\end{bmatrix}&\sim\mathcal{N}_{(p1+p2)}\left(\begin{bmatrix}\mathbf{W}^{(1)}\\\
\mathbf{W}^{(2)}\end{bmatrix}\mathbf{z}_{n},\begin{bmatrix}\boldsymbol{\Sigma}^{(1)}&\boldsymbol{0}\\\
\boldsymbol{0}&\boldsymbol{\Sigma}^{(2)}\end{bmatrix}\right)\end{split}$

where ${\mathbf{y}_{n}}^{(1)}\in\mathbb{R}^{p_{1}}$ and
${\mathbf{y}_{n}}^{(2)}\in\mathbb{R}^{p_{2}}$, each denote a random vector of
observations for Data set 1 and Data set 2, respectively. Note that by
replacing the block diagonal covariance matrix of the data with a diagonal
matrix, the probabilistic CCA model reduces to a standard factor analysis
model!

They also show that since any covariance matrix of rank $D$ can be written as
the following decomposition of the marginal covariance in the marginal joint
distribution of the data

$\displaystyle\begin{bmatrix}{\mathbf{y}_{n}}^{(1)}\\\
{\mathbf{y}_{n}}^{(2)}\end{bmatrix}$
$\displaystyle\sim\mathcal{N}\left(\boldsymbol{0},\begin{bmatrix}\mathbf{W}^{(1)}{\mathbf{W}^{(1)}}^{\intercal}+\boldsymbol{\Sigma}^{(1)}&\mathbf{W}^{(1)}{\mathbf{W}^{(2)}}^{\intercal}\\\
\mathbf{W}^{(2)}{\mathbf{W}^{(1)}}^{\intercal}&\mathbf{W}^{(2)}{\mathbf{W}^{(2)}}^{\intercal}+\boldsymbol{\Sigma}^{(2)}\end{bmatrix}\right)$

it follows that the MLE solution for the joint covariance matrix is such that
the cross-covariance matrix is of rank $D$ and consists of the canonical
correlation directions that would be obtained in regular CCA.

###  5.8 Inter-Battery Factor Analysis (IBFA)

If we split the vector of latent variables $\mathbf{z}_{n}$ into three subsets
$\\{\mathbf{z}_{n,0},\mathbf{z}_{n,1},\mathbf{z}_{n,1}\\}$ such that the
latent vector $\mathbf{z}_{n,0}$ captures the common structure shared across
both groups of variables, ${\mathbf{y}_{n}}^{(1)}$ and
${\mathbf{y}_{n}}^{(2)}$, while $\mathbf{z}_{n,1}$ and $\mathbf{z}_{n,1}$ each
captures the latent structure specific to each group, then CCA model can be
modified to have an equivalent specification to the Inter-Battery Factor
Analysis (IBFA) model (Tucker, 1958; Browne, 1979).

$\displaystyle\begin{split}\begin{bmatrix}\mathbf{z}_{n,0}\\\
\mathbf{z}_{n,1}\\\
\mathbf{z}_{n,2}\end{bmatrix}&\sim\mathcal{N}_{d}\left(\boldsymbol{0},\mathbf{I}\right)\quad
n=1,\cdots,N\\\
\begin{bmatrix}\mathbf{x}^{(1)}\mid\mathbf{z}_{n,0},\mathbf{z}_{n,1}\\\
\mathbf{x}^{(2)}\mid\mathbf{z}_{n,0},\mathbf{z}_{n,2}\end{bmatrix}&\sim\mathcal{N}\left(\begin{bmatrix}\mathbf{W}^{(1,0)}&\mathbf{W}^{(1,1)}&\boldsymbol{0}\\\
\mathbf{W}^{(2,0)}&\boldsymbol{0}&\mathbf{W}^{(2,2)}\end{bmatrix}\begin{bmatrix}\mathbf{z}_{n,0}\\\
\mathbf{z}_{n,1}\\\
\mathbf{z}_{n,2}\end{bmatrix},\begin{bmatrix}\boldsymbol{\Sigma}^{(1)}&\boldsymbol{0}\\\
\boldsymbol{0}&\boldsymbol{\Sigma}^{(2)}\end{bmatrix}\right)\end{split}$

Note that the above generative formulation makes it apparent that IBFA can be
viewed as a special case of CCA, one obtained when we restrict the projection
matrix $\mathbf{W}$ to have a the structure given above; namely

$\displaystyle\mathbf{W}=\begin{bmatrix}\mathbf{W}^{(1)}\\\
\mathbf{W}^{(2)}\end{bmatrix}=\begin{bmatrix}\mathbf{W}^{(1,0)}&\mathbf{W}^{(1,1)}&\boldsymbol{0}\\\
\mathbf{W}^{(2,0)}&\boldsymbol{0}&\mathbf{W}^{(2,2)}\end{bmatrix}$

###  5.9 Multi-Battery Factor Analysis (MBFA)

The Inter-Battery Factor Analysis model can be readily extended to multiple
groups of variables. So if we let
${\mathbf{y}_{n}}^{(g)}\in\mathbb{R}^{p_{g}}$ for $g=1,\cdots,G$ denote the
groups of variables (i.e. data sets), the vector of observations can be
stacked as $P=\sum_{g}^{G}p_{g}$ dimensional vector and the Multi-Battery
Factor Analysis (MBFA) model (Browne, 1980) can be expressed as follows.

$\displaystyle\begin{split}\mathbf{z}_{n}&\sim\mathcal{N}_{d}\left(\boldsymbol{0},\mathbf{I}\right)\quad
n=1,\cdots,N\\\ \begin{bmatrix}{\mathbf{y}_{n}}^{(1)}\mid\mathbf{z}_{n}\\\
{\mathbf{y}_{n}}^{(2)}\mid\mathbf{z}_{n}\\\ \vdots\\\
{\mathbf{y}_{n}}^{(G)}\mid\mathbf{z}_{n}\\\
\end{bmatrix}&\sim\mathcal{N}_{P}\left(\begin{bmatrix}\mathbf{W}^{(1)}\\\
\mathbf{W}^{(2)}\\\ \vdots\\\
\mathbf{W}^{(G)}\end{bmatrix}\mathbf{z}_{n},\begin{bmatrix}\boldsymbol{\Sigma}^{(1)}&0&\ldots&0\\\
0&\boldsymbol{\Sigma}^{(2)}&\ldots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&\ldots&0&\boldsymbol{\Sigma}^{(G)}\end{bmatrix}\right)\\\ \end{split}$

Again note that by replacing the block diagonal covariance matrix of the data
with a diagonal matrix, MBFA reduces to a standard factor analysis model.

###  5.10 Group Factor Analysis (GFA)

Klami et al. (2015) proposed Group Factor Analysis (GFA) as a generalization
of multi-battery factor analysis (MBFA). The generative model does still have
the same general structure of the preceding latent variable models

$\displaystyle\begin{split}\mathbf{z}_{n}&\sim\mathcal{N}_{d}\left(\boldsymbol{0},\mathbf{I}\right)\quad
n=1,\cdots,N\\\ \begin{bmatrix}{\mathbf{y}_{n}}^{(1)}\mid\mathbf{z}_{n}\\\
{\mathbf{y}_{n}}^{(2)}\mid\mathbf{z}_{n}\\\ \vdots\\\
{\mathbf{y}_{n}}^{(G)}\mid\mathbf{z}_{n}\\\
\end{bmatrix}&\sim\mathcal{N}_{P}\left(\begin{bmatrix}\mathbf{W}^{(1)}\\\
\mathbf{W}^{(2)}\\\ \vdots\\\
\mathbf{W}^{(G)}\end{bmatrix}\mathbf{z}_{n},\begin{bmatrix}\sigma^{2}_{1}\mathbf{I}&0&\ldots&0\\\
0&\sigma^{2}_{2}\mathbf{I}&\ldots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&\ldots&0&\sigma^{2}_{G}\mathbf{I}\end{bmatrix}\right)\end{split}$

where

$\displaystyle\mathbf{W}=\begin{bmatrix}\mathbf{W}^{(1)}\\\
\mathbf{W}^{(2)}\\\ \vdots\\\
\mathbf{W}^{(G)}\end{bmatrix}=\begin{bmatrix}\mathbf{w}^{(1)}_{1}&\mathbf{w}^{(1)}_{2}&\dots&\mathbf{w}^{(1)}_{D}\\\
\mathbf{w}^{(2)}_{1}&\mathbf{w}^{(2)}_{2}&\dots&\mathbf{w}^{(2)}_{D}\\\
\vdots&\vdots&\vdots&\vdots\\\
\mathbf{w}^{(G)}_{1}&\mathbf{w}^{(G)}_{2}&\dots&\mathbf{w}^{(G)}_{D}\end{bmatrix}$

There are two differences however. First, GFA restricts the group specific
error covariance to isotropic structure. That is,
$\boldsymbol{\Sigma}^{(g)}=\sigma^{2}_{g}\mathbf{I}$ for $g=1,\cdots,G$.
Second, in order to model the dependency structure between any subsets of
groups and resolve indeterminacy issues, GFA imposes a structural sparsity
prior over the columns of the projection matrix $\mathbf{W}$. The Automatic
Relevance Determination (ARD) prior (MacKay, 1995b) over the $G\times D$ group
specific columns of the projection matrix $\mathbf{W}$ can be specified as
follows.

$\displaystyle\begin{split}\mathbf{w}^{(g)}_{d}&\sim\mathcal{N}\left(\boldsymbol{0},{(\alpha^{(g)}_{d})}^{-1}\mathbf{I}\right)\quad
g=1,\cdots,G;\quad d=1,\cdots,D\\\
log(\boldsymbol{\alpha})&=UV^{\intercal}+\mu_{u}\textbf{1}^{\intercal}+\textbf{1}\mu_{v}^{\intercal}\\\
p(U)&=\prod_{g}^{G}{\prod_{r}^{R}{\mathcal{N}(u_{g,r}|0,(\lambda=0.1)^{-1})}}\\\
p(V)&=\prod_{d}^{D}{\prod_{r}^{R}{\mathcal{N}(v_{d,r}|0,(\lambda=0.1)^{-1})}}\end{split}$

for $\boldsymbol{\alpha}\in\mathbb{R}^{G\times D}$, $U\in\mathbb{R}^{G\times
R}$, and $V\in\mathbb{R}^{D\times R}$ and $R\ll min(G,D)$.

###  5.11 Structural Equation Models (SEM)

The preceding models can be extended by allowing the latent variables to
regress on each other. The resulting general model is known as a Structural
Equation Model (SEM) and was first introduced by Jöreskog (1970). SEM models
are constructed using algebraic equations, not by formulating a probabilistic
statistical model as is done in Bayesian Graphical Models (Pearl, 2014) for
example.

A popular general formulation of SEMs is given by the LISREL model. The LISREL
model consists of two parts. The first part is a structural component in which
dependencies among the latent variables are created by regressing them on each
other. The vector of latent variables we introduce into the model is divided
into two groups of independent and dependent variable
$\mathbf{z}_{n}=(\mathbf{z}_{n,1}^{\intercal},\mathbf{z}_{n,2}^{\intercal})$
of dimensions $D_{1}$ and $D_{2}$, respectively, such that $D_{1}+D_{2}=D$.
The second part consists of a measurement model that is similar in form to the
factor analysis model (i.e.
$\mathbf{y}_{n}=\mathbf{W}\mathbf{z}_{n}+\boldsymbol{\epsilon}_{n}$) but for
which we also partition the $P$-dimensional observation vector (and
correspondingly the error vector) into two groups
$\mathbf{y}_{n}=(\mathbf{y}_{n,1}^{\intercal},\mathbf{y}_{n,2}^{\intercal})$
and
$\boldsymbol{\epsilon}_{n}=(\boldsymbol{\epsilon}_{n,1}^{\intercal},\boldsymbol{\epsilon}_{n,2}^{\intercal})$
such that $\mathbf{y}_{n,1}$ depends only on $\mathbf{z}_{n,2}$ and
$\mathbf{y}_{n,2}$ depends only on $\mathbf{z}_{n,1}$. The dimension of two
groups is $P_{1}$ and $P_{2}$, such that $P_{1}+P_{2}=P$. The general
algebraic form of the LISREL model is as follows.

$\displaystyle\begin{split}\begin{bmatrix}\mathbf{y}_{n,1}\\\
\mathbf{y}_{n,2}\end{bmatrix}&=\begin{bmatrix}\mathbf{W}_{1}&\boldsymbol{0}\\\
\boldsymbol{0}&\mathbf{W}_{2}\end{bmatrix}\begin{bmatrix}\mathbf{z}_{n,1}\\\
\mathbf{z}_{n,2}\end{bmatrix}+\begin{bmatrix}\boldsymbol{\epsilon}_{n,1}\\\
\boldsymbol{\epsilon}_{n,2}\end{bmatrix}\\\
\mathbf{B}\mathbf{z}_{n,2}&=\mathbf{C}\mathbf{z}_{n,1}+\boldsymbol{\xi}_{n}\end{split}$

where it is assumed that B is non-singular and

$\displaystyle\begin{split}&\operatorname{\mathbb{E}}(\mathbf{y}_{n})=\mathbf{0};\operatorname{\mathbb{E}}(\mathbf{z}_{n})=\mathbf{0}\;\
\operatorname{\mathbb{E}}(\boldsymbol{\epsilon}_{n})=\mathbf{0};\
\operatorname{\mathbb{E}}(\boldsymbol{\xi}_{n})=\mathbf{0}\\\
&\operatorname*{\mathrm{Cov}}(\boldsymbol{\epsilon}_{n})=\Psi;\
\operatorname*{\mathrm{Cov}}(\mathbf{z}_{n,1})=\Phi_{z_{1}};\
\operatorname*{\mathrm{Cov}}(\boldsymbol{\xi}_{n})=\Phi_{\boldsymbol{\xi}}\end{split}$

and the covariance matrices are diagonal. The regressions imply the following
covariance matrix $\boldsymbol{\Sigma}_{\mathbf{y}\mathbf{y}}$ for the
observations.

$\displaystyle\begin{bmatrix}\mathbf{W}_{1}\Phi_{z_{1}}\mathbf{W}_{1}^{\intercal}+\Psi_{1}&\mathbf{W}_{1}\Phi_{z_{1}}\mathbf{C}^{\intercal}{(\mathbf{B}^{\intercal})}^{-1}\mathbf{W}_{2}^{\intercal}\\\
\mathbf{W}_{2}\Phi_{z_{1}}\mathbf{C}^{\intercal}{(\mathbf{B}^{\intercal})}^{-1}\mathbf{W}_{1}^{\intercal}&\mathbf{W}_{2}(\mathbf{B}^{-1}\mathbf{C}^{\intercal}\Phi_{z_{1}}\mathbf{C}{(\mathbf{B}^{\intercal})}^{-1}+\mathbf{B}^{-1}\Phi_{\boldsymbol{\xi}}{(\mathbf{B}^{\intercal})}^{-1})\mathbf{W}_{2}^{\intercal}+\Psi_{2}\end{bmatrix}$

Notice that how the implied covariance structure of the data is determined by
the coefficient matrices ( i.e. $\mathbf{B}$, $\mathbf{C}$, and $\mathbf{W}$)
and the covariance matrices (i.e. $\Psi$, $\Phi_{z_{1}}$ and
$\Phi_{\boldsymbol{\xi}}$). More importantly, notice also that an SEM model
implicitly assumes that the relationships between the observations are linear
since covariances are measures of linear association. Another aspect of SEM
models that is worth pointing out is that although we can assess whether a
particular implied decomposition is consistent with the data, we cannot
determine whether it is the only decomposition that is equally consistent with
the data. That is, the same data can be generated given different latent
structural relations (Bartholomew et al., 2011).

##### Generalized Structured Component Analysis (GSCA)

A generalization of the SEM framework has recently been proposed by Hwang and
Takane (2004). The goal of the generalization is to unify two distinct
approaches to SEM; namely the factor-based and component-based (i.e. partial
least squares). The Generalized Structured Component Analysis (GSCA) model
they propose has the following structural formulation.

$\displaystyle\begin{split}\begin{bmatrix}\mathbf{I}\\\
\mathbf{W}\end{bmatrix}\mathbf{y}_{n}=\begin{bmatrix}\mathbf{C}\\\
\mathbf{B}\end{bmatrix}\mathbf{W}\mathbf{y}_{n}+\begin{bmatrix}\boldsymbol{\epsilon}_{n}\\\
\boldsymbol{\xi}_{n}\end{bmatrix}\end{split}$ (27)

where the latent variables are obtained as component scores using what the
authors refer to as weighted relation model
$\mathbf{z}_{n}=\mathbf{W}\mathbf{y}_{n}$. The assumption regarding the means
and covariance of the variables are identical to the LISERL model. It is
interesting to note that the weighted relation model equation is no different
conceptually than the encoding transformation that Izenman (2008) used to
construct latent variable models from the reduced rank regression framework (
Section 5.3).

We argue that by relying on complicated algebraic equations to induce linear
dependencies in the latent structure, the SEM framework risks losing
statistical rigor especially when modeling complex phenomena. Real data, most
often than not, is noisy and characterized by complex nonlinear dependencies.
We think that we have a better and more statistically rigorous alternative in
the probabilistic generative perspective that approaches the modeling endeavor
by first proposing a probabilistic model for all the observed variables and
any latent variables that we wish to incorporate in the model. As part of the
probabilistic model, we can also specify functional relations between the
variables in order to capture any nonlinear dependencies in the data or to
construct a rich latent space with complex structure. For example, a recent
method, the Inverse Autoregressive Flow (IAF)(Kingma et al., 2016), has been
proposed in the generative modeling literature that enables us to create
models with rich latent structure. The key idea behind the IAF is that any
multivariate normal variable with diagonal covariance can be transformed into
a multivariate normal with linear dependencies (i.e. full covariance matrix)
by expressing it as an autoregressive model

$\displaystyle
z_{d}=\mu_{d}(\mathbf{z}_{<d})+\sigma_{d}(\mathbf{z}_{<d})\ast\epsilon$

where $\mu_{d}$ and $\sigma_{d}$ are the conditional mean vector and
covariance matrix of $z_{d}$ on the other elements of the latent vector
$\mathbf{z}_{n}=(z_{1},z_{2},\cdots,z_{d},\cdots,z_{D})$. Accordingly, a
general SEM model can be conceptually viewed as the following generative
model.

$\displaystyle\begin{split}\mathbf{y}_{n}\mid\mathbf{z}_{n}&\sim\mathcal{N}_{p}\Big{(}\mathbf{W}\mathbf{z}_{n},\operatorname{Diag}(\sigma^{2})\Big{)}\\\
\mathbf{z}_{n}&\sim\mathcal{N}_{r}(\boldsymbol{0},\boldsymbol{\Sigma}_{z})\end{split}$

By including structural latent regressions, a specific SEM model would then
correspond to a particular structure of the covariance matrix
$\boldsymbol{\Sigma}_{z}$. In the next section, we briefly describe Deep
Latent Gaussian Models (DLGMs) as an alternative general framework to SEM, one
that is more statistically disciplined and potentially much more powerful,
especially when large amounts of data is available.

###  5.12 Deep Latent Gaussian Models (DLGM)

A Deep Latent Gaussian Model (DLGM) (Rezende et al., 2014) is a directed
graphical model with Gaussian latent variables that are stacked in $L$
stochastic layers and that are transformed at each layer by a differentiable
nonlinear function parameterized by a neural network
$\mathbf{NN}(\mathbf{z};\boldsymbol{\theta})$ consisting of $K^{(l)}$
deterministic hidden layers. Whereas the latent variables in the first DLGMs
such as the Variational Autoencoder (VAE) (Kingma and Welling, 2013) are
assumed to be independent, recent extensions (Kingma et al., 2016; Maaløe et
al., 2016) have introduced a deep latent hierarchy that endows the latent
space with complex rich structure that can capture complex patterns of
dependencies between the latent variables. A general DLGM model has the
following generative formulation.

$\displaystyle\begin{split}\mathbf{z}_{n}^{(L)}\ &\sim\
\mathcal{N}_{d_{(L)}}\left(\boldsymbol{0},\,\mathbf{I}\right)\quad
n=1,\cdots,N\\\ \mathbf{z}_{n}^{(l)}\ &\sim\
\mathcal{N}_{d_{(l)}}\Big{(}\mathbf{z}_{n}^{(l)}\mid\mathbf{NN}^{(l)}(\mathbf{z}_{n}^{(l+1)};\boldsymbol{\theta}^{(l)}),\,\Sigma^{(l)}\Big{)}\quad
l=1,\cdots,L-1\\\ \mathbf{y}_{n}\mid\mathbf{z}_{n}&\sim\
\mathbf{Expon}\Big{(}\mathbf{y}_{n}\mid\mathbf{NN}^{(0)}(\mathbf{z}_{n}^{(1)};\boldsymbol{\theta}^{(0)})\Big{)}\end{split}$

where

$\displaystyle\begin{split}\mathbf{NN}(\mathbf{z};\theta)&=h_{K}\circ
h_{K-1}\circ\ldots\circ h_{0}(\mathbf{z})\\\ \quad
h_{k}(\mathbf{y})&=\sigma_{k}(\mathbf{W}^{(k)}\mathbf{y}+\mathbf{b}^{(k)})\\\
\boldsymbol{\theta}&=\\{(\mathbf{W}^{(k)},\mathbf{b}^{(k)})\\}_{k=0}^{K}\end{split}$

Such that the observation model $\mathbf{Expon}$ refers to the exponential
family of models and the neural network
$\mathbf{NN}(\mathbf{z};\boldsymbol{\theta})$ consists of the composition of
$K$ non-linear transformation functions, each of which is represented by a
hidden layer $h_{k}(\mathbf{y})$. The function $\sigma_{k}$ is an element-wise
activation function such as the sigmoid function or the rectified linear unit
(ReLU). The constant $d_{(l)}$ is the number of latent variables at the $l$
layer so that the total number of latent variables equals
$d=\sum_{l=1}^{L}d_{(l)}$

A Deep Latent Gaussian Model is considered a deep model in two distinct
senses. First, the latent variables that make up the stochastic layers can be
stacked in a deep hierarchy. Second, the neural networks that make up the
deterministic nonlinear transformations can be made up of several hidden
layers. In contrast, a nonlinear regression model is a shallow nonlinear model
in the sense that it uses only one non-linear transformation while a factor
analysis model is a shallow linear model in the sense that the latent
variables $\mathbf{z}_{n}$ are linearly projected by a matrix $\mathbf{W}$ to
a higher dimensional vector.

## References

  * Advani and Ganguli (2016) Madhu Advani and Surya Ganguli.  Statistical Mechanics of Optimal Convex Inference in High Dimensions.  _Physical Review X_ , 6(3):031034, 2016.
  * Airy (1861) George Biddell Airy.  _On the algebraical and numerical theory of errors of observations and the combination of observations_.  Macmillan&Company, 1861.
  * Athreya and Lahiri (2006) Krishna B Athreya and Soumendra N Lahiri.  _Measure theory and probability theory_.  Springer Science & Business Media, 2006.
  * Bach and Jordan (2005) Francis R Bach and Michael I Jordan.  A probabilistic interpretation of canonical correlation analysis.  techreport 688, Department of Statistics, University of California, Berkeley, Berkeley, 2005.
  * Banerjee et al. (2014) Sudipto Banerjee, Bradley P Carlin, and Alan E Gelfand.  _Hierarchical modeling and analysis for spatial data_.  Crc Press, 2014.
  * Bartholomew et al. (2011) David J Bartholomew, Martin Knott, and Irini Moustaki.  _Latent variable models and factor analysis: A unified approach_ , volume 904.  John Wiley & Sons, 2011.
  * Birnbaum (1968) A. Birnbaum.  Some latent trait models and their use in inferring an examinee’s ability., 1968.
  * Bishop et al. (1998) Christopher M Bishop, Markus Svensén, and Christopher KI Williams.  Gtm: the generative topographic mapping.  _Neural computation_ , 10(1):215–234, 1998.
  * Blei et al. (2016) David M Blei, Alp Kucukelbir, and Jon D McAuliffe.  Variational inference: A review for statisticians.  _Journal of the American Statistical Association (to appear). arXiv preprint arXiv:1601.00670_ , 2016.
  * Bollen (2002) Kenneth A Bollen.  Latent variables in psychology and the social sciences.  _Annual review of psychology_ , 53(1):605–634, 2002.
  * Browne (1979) Michael W Browne.  The maximum-likelihood solution in inter-battery factor analysis.  _British Journal of Mathematical and Statistical Psychology_ , 32(1):75–86, 1979.
  * Browne (1980) Michael W Browne.  Factor analysis of multiple batteries by maximum likelihood.  _British Journal of Mathematical and Statistical Psychology_ , 33(2):184–199, 1980.
  * Chung et al. (2015) Junyoung Chung, Kyle Kastner, Laurent Dinh, Kratarth Goel, Aaron C Courville, and Yoshua Bengio.  A recurrent latent variable model for sequential data.  In _Advances in neural information processing systems_ , pages 2980–2988, 2015.
  * Diggle and Gratton (1984) Peter J Diggle and Richard J Gratton.  Monte carlo methods of inference for implicit statistical models.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , pages 193–227, 1984.
  * Durrett (2010) Rick Durrett.  _Probability: theory and examples_.  Cambridge university press, 2010.
  * Gelman et al. (2014) Andrew Gelman, John B Carlin, Hal S Stern, and Donald B Rubin.  _Bayesian data analysis_ , volume 2.  Chapman & Hall/CRC Boca Raton, FL, USA, 2014.
  * Gemici et al. (2017) Mevlana Gemici, Chia-Chun Hung, Adam Santoro, Greg Wayne, Shakir Mohamed, Danilo J Rezende, David Amos, and Timothy Lillicrap.  Generative temporal models with memory.  _arXiv preprint arXiv:1702.04649_ , 2017.
  * Ghosh and Ramamoorthi (2006) J.K. Ghosh and R.V. Ramamoorthi.  _Bayesian Nonparametrics_.  Springer Science & Business Media, 2006.
  * Gómez et al. (1998) E Gómez, MA Gomez-Viilegas, and JM Marin.  A multivariate generalization of the power exponential family of distributions.  _Communications in Statistics-Theory and Methods_ , 27(3):589–600, 1998.
  * Goodfellow et al. (2014) Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  In _Advances in neural information processing systems_ , pages 2672–2680, 2014.
  * Gould (1996) Stephen Jay Gould.  _The mismeasure of man_.  WW Norton & Company, 1996.
  * Hägglund (2001) G Hägglund.  Milestones in the history of factor analysis.  In Robert Cudeck, K. G. Jöreskog, and Dag Sörbom, editors, _Structural equation modeling, present and future: A festschrift in honor of Karl Jöreskog_ , pages 11–38. Scientific Software International, Lincolnwood, IL, 2001.
  * Hilborn and Mangel (1997) Ray Hilborn and Marc Mangel.  _The ecological detective: confronting models with data_ , volume 28.  Princeton University Press, 1997.
  * Hotelling (1933) Harold Hotelling.  Analysis of a complex of statistical variables into principal components.  _Journal of Educational Psychology_ , 24(6):417–441, 1933.
  * Hotelling (1936) Harold Hotelling.  Relations Between Two Sets of Variates.  _Biometrika_ , 28(3/4):321, dec 1936.
  * Hwang and Takane (2004) Heungsun Hwang and Yoshio Takane.  Generalized structured component analysis.  _Psychometrika_ , 69(1):81–99, 2004.
  * Hwang and Takane (2014) Heungsun Hwang and Yoshio Takane.  Generalized structured component analysis: A component-based approach to structural equation modeling, 2014.
  * Hyvärinen (2015) Aapo Hyvärinen.  A unified probabilistic model for independent and principal component analysis.  _Advances in Independent Component Analysis and Learning Machines_ , pages 1–9, 2015.
  * Izenman (1975) Alan Julian Izenman.  Reduced-rank regression for the multivariate linear model.  _Journal of multivariate analysis_ , 5(2):248–264, 1975.
  * Izenman (2008) Alan Julian Izenman.  _Modern multivariate statistical techniques_.  Springer, 2008.
  * Jaynes (2003) Edwin T Jaynes.  _Probability theory: The logic of science_.  Cambridge university press, 2003.
  * Jebara (2012) Tony Jebara.  _Machine learning: discriminative and generative_ , volume 755.  Springer Science & Business Media, 2012.
  * Jöreskog (1970) K. G. Jöreskog.  A general method for analysis of covariance structures.  _Biometrika_ , 57(2):239–251, 1970.
  * Kallenberg (2006) Olav Kallenberg.  _Probabilistic symmetries and invariance principles_.  Springer Science & Business Media, 2006.
  * Kant and Guyer (1998) Immanuel Kant and Paul Guyer.  _Critique of pure reason_.  Cambridge University Press, 1998.
  * Kennedy and O’Hagan (2001) Marc C. Kennedy and Anthony O’Hagan.  Bayesian calibration of computer models.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 63, 2001.
  * Kingma and Welling (2013) Diederik P Kingma and Max Welling.  Auto-Encoding Variational Bayes.  _ICLR_ , (Ml):1–14, dec 2013.
  * Kingma et al. (2016) Diederik P Kingma, Tim Salimans, and Max Welling.  Improving variational inference with inverse autoregressive flow.  _arXiv preprint arXiv:1606.04934v2_ , 2016.
  * Klami et al. (2015) Arto Klami, Seppo Virtanen, Eemeli Leppäaho, and Samuel Kaski.  Group factor analysis.  _IEEE transactions on neural networks and learning systems_ , 26(9):2136–2147, 2015.
  * Lake et al. (2015) Brenden M Lake, Ruslan Salakhutdinov, and Joshua B Tenenbaum.  Human-level concept learning through probabilistic program induction.  _Science_ , 350(6266):1332–1338, 2015.
  * Lauritzen (1996) Steffen L Lauritzen.  _Graphical models_ , volume 17.  Clarendon Press, 1996.
  * Lazarsfeld (1950) Paul F Lazarsfeld.  The logical and mathematical foundation of latent structure analysis.  _Measurement and prediction_ , 4:362–412, 1950.
  * Lord (1952) Frederic M Lord.  The relation of test score to the trait underlying the test.  _ETS Research Report Series_ , 1952(2):517–549, 1952.
  * Maaløe et al. (2016) Lars Maaløe, Casper Kaae Sønderby, Søren Kaae Sønderby, and Ole Winther.  Auxiliary deep generative models.  _arXiv preprint arXiv:1602.05473_ , 2016.
  * MacKay (1995a) David JC MacKay.  Bayesian neural networks and density networks.  _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , 354(1):73–80, 1995a.
  * MacKay (1995b) David JC MacKay.  Probable networks and plausible predictions - a review of practical bayesian methods for supervised neural networks.  _Network: Computation in Neural Systems_ , 6(3):469–505, 1995b.
  * Mansinghka et al. (2016) Vikash Mansinghka, Patrick Shafto, Eric Jonas, Cap Petschulat, Max Gasner, and Joshua B Tenenbaum.  Crosscat: A fully bayesian nonparametric method for analyzing heterogeneous, high dimensional data.  _Journal of Machine Learning Research_ , 17:1–49, 2016\.
  * Muthén (2002) Bengt O Muthén.  Beyond sem: General latent variable modeling.  _Behaviormetrika_ , 29(1):81–117, 2002.
  * Pearl (2014) Judea Pearl.  _Probabilistic reasoning in intelligent systems: networks of plausible inference_.  Morgan Kaufmann, 2014.
  * Pearson (1901) Karl Pearson.  LIII. On lines and planes of closest fit to systems of points in space.  _Philosophical Magazine Series 6_ , 2(11):559–572, nov 1901.
  * Pritchard et al. (1999) Jonathan K Pritchard, Mark T Seielstad, Anna Perez-Lezaun, and Marcus W Feldman.  Population growth of human y chromosomes: a study of y chromosome microsatellites.  _Molecular biology and evolution_ , 16(12):1791–1798, 1999.
  * Rabe-Hesketh et al. (2004) Sophia Rabe-Hesketh, Anders Skrondal, and Andrew Pickles.  Generalized multilevel structural equation modeling.  _Psychometrika_ , 69(2):167–190, 2004.
  * Ranganath et al. (2015) Rajesh Ranganath, Linpeng Tang, Laurent Charlin, and David M Blei.  Deep exponential families.  In _AISTATS_ , 2015.
  * Rezende et al. (2014) Danilo Jimenez Rezende, Shakir Mohamed, and Daan Wierstra.  Stochastic backpropagation and approximate inference in deep generative models.  _arXiv preprint arXiv:1401.4082_ , 2014.
  * Rezende et al. (2016) Danilo Jimenez Rezende, Shakir Mohamed, Ivo Danihelka, Karol Gregor, and Daan Wierstra.  One-shot generalization in deep generative models.  _arXiv preprint arXiv:1603.05106_ , 2016.
  * Roweis and Ghahramani (1999) Sam Roweis and Zoubin Ghahramani.  A Unifying Review of Linear Gaussian Models.  _Neural Computation_ , 1999.
  * Searle et al. (2009) Shayle R Searle, George Casella, and Charles E McCulloch.  _Variance components_ , volume 391.  John Wiley & Sons, 2009.
  * Skrondal and Rabe-Hesketh (2004) Anders Skrondal and Sophia Rabe-Hesketh.  _Generalized Latent Variable Modeling_ , volume 20041561 of _C &H/CRC Monographs on Statistics & Applied Probability_.  Chapman and Hall/CRC, may 2004.
  * Sobel (1994) Michael E Sobel.  Causal inference in latent variable models.  1994\.
  * Spearman (1904) C Spearman.  ”General Intelligence,” Objectively Determined and Measured.  _The American Journal of Psychology_ , 15(2):201, apr 1904.
  * Tavaré et al. (1997) Simon Tavaré, David J Balding, Robert C Griffiths, and Peter Donnelly.  Inferring coalescence times from dna sequence data.  _Genetics_ , 145(2):505–518, 1997.
  * Thomson (1916) Godfrey H. Thomson.  A HIERARCHY WITHOUT A GENERAL FACTOR.  _British Journal of Psychology, 1904-1920_ , 8(3):271–281, sep 1916.
  * Thurstone and L.L. (1947) Thurstone and L.L.  Multiple Factor Analysis., 1947.
  * Tipping and Bishop (1999) Michael E. Tipping and Christopher M. Bishop.  Probabilistic Principal Component Analysis.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 61(3):611–622, aug 1999.
  * Tobin (1958) James Tobin.  Estimation of relationships for limited dependent variables.  _Econometrica: journal of the Econometric Society_ , pages 24–36, 1958.
  * (66) Matti Tornio, Antti Honkela, and Juha Karhunen.  Time series prediction with variational bayesian nonlinear state-space models.
  * Tran et al. (2017) Dustin Tran, Rajesh Ranganath, and David M Blei.  Deep and hierarchical implicit models.  _arXiv preprint arXiv:1702.08896_ , 2017.
  * Tucker (1958) Ledyard R Tucker.  An inter-battery method of factor analysis.  _Psychometrika_ , 23(2):111–136, 1958.
  * Wright (1921) Sewall Wright.  Correlation and causation.  _Journal of agricultural research_ , 20(7):557–585, 1921.