Multivariate Normal Distribution#
from IPython.display import display, HTML
display(HTML(“””
“””))
display(HTML(“””
“””))
Probability defined over vectors and the MVN#
Up until this point we have discussed assigning probabilities to a single random variable, or to a sequence of random variables such that every random variable in the sequence is independent of one another. We have not yet discussed (the very real) case that two or more random variables are related to one another and how to express these relationships.
Given a set of random variables \(X_{1}, X_{2}, \cdots, X_{n}\), a succint approach for defining probabilities to these random variables is by collecting them into a vector
and assigning probabilties on the space \(\mathbb{R}^{n}\). In other words, we will begin to explore mathematical statements like
, or “What is the probability that the first random variable \(X_{1}\) equals the value 1 and the first random variable \(X_{2}\) equals the value 2 and the third random variable \(X_{3}\) equals the value 3 “
Before we begin a serious discussion of probability assignments like this, we need mathematical tools from linear algebra.
Review of vector and matrix algebra#
Vectors#
Vectors are the fundemental mathemtical object that we will work with in this unit. Vectors are defined as a list of items where each item is an element from the real line. For example
is the vector \(v_{1}\) that contains the value one in the first position, value negative one in the second position and the value one-half in the third position.
Vectors can interact with elements from \(\mathbb{R}\) called scalars as well as interact with one another.
\(\alpha x\) and \(\alpha + x\)#
Given a scalar \(\alpha\), vector \(x \in \mathbb{R}^{n}\), then the new vector \(\alpha x\) is defined as
and the new vector \(\alpha + x\) is defined as
\(x+y\) and \(xy\)#
Given vectors \(x\) and \(y\) both in \(\mathbb{R}^{n}\) then the new vector \(x+y\) is defined as
and the new vector \(xy\) is defined as
An approach to defining a set of vectors#
On occasion it is important to describe a set of vectors, called a vector space. The canonical approach to describe a set of vectors is to generate all vectors in a vector space \(V\) as multiplicative and additive combinations of some core set of vectors. In otherwords, we can generate a set of vectors by picking a small number of vectors \(y_{1}, y_{2}, y_{3}\) and building a vector space \(V\) as
For example, we may decide to generate the vector space \(V\) using these two vectors:
Then in the space we created contains vectors like
or the vector
or
One can imagine defining then different vector spaces, and so then it is natural to ask about functions from one vector space to a second.
In otherwords, we can define functions over vectors like
However, an important class of functions are linear functions. A function is linear if
As it turns out, linear functions have a unique representation. Consider a vector space \(V\) generated by the two vectors \(x_{1}\) and \(x_{2}\). In addition, consider a function \(g\) that is linear and maps the vector space \(V\) to the new space \(Q\).
Then for a vector \(x\) in the space \(V\), we can write
If we wanted to apply the function \(g\) to the vector \(x\) then
However, remember that \(g\) is a linear function. This means that we can simplify the above as
But, because \(g\) maps vectors from \(V\) to \(Q\) these new vectors \(g\left(x_{1}\right)\) and \(g\left(x_{2}\right)\) are members of \(Q\). If we suppose that every vector in \(Q\) can be generated from (say) the three vectors \(q_{1}, q_{2}, q_{3}\) then we can rewrite \(g\left(x_{1}\right)\) in terms of these three \(q\)s as well as \(g\left(x_{2}\right)\).
In other words,
or \(g(x)\) can be represented as the vector
Matrices#
A matrix is defined as a stack or ordered collection of vectors. Matrices are written as a rectangular grid of values such as
We can refer to individual elements inside the matrix using this notation \(A_{ij}\). This refers to the element located in row \(i\) and column \(j\). For example, \(A_{24}\) is the value \(a_{24}\) above. We can refer to an entire row of \(A\) usign this notation \(A_{(i,:)}\) and an entire column of \(A\) using this notation \(A_{(:,j)}\).
Ax#
The matrix \(A\) with \(m\) rows and \(n\) columns, referred to as its dimension, maps a vector of length \(n\) to a new vector of length \(m\). From our discussion of linear functions above, we see that, for example,
, or we multiply the first column in \(A\) by the first element in \(x\) plus the second column in \(A\) with the second \(element\) in \(x\) and so on.
Another approach towards matrix multiplication uses the inner product. Given a vector \(x\) and \(y\), the inner product is defined as
It is important to note some properties of the inner product. The inner product is a linear operator or, for constant \(\alpha\),
These properties come from the fact that the inner product is a summation over the individual elements in \(x\) and \(y\). With the inner product, we can write matrix multiplication using the rows of \(A\) as
AB#
Recall that the matrix \(A\) was a represenation of a linear function that input a vector \(x\) and output a vector \(y\). Another common operation is function composition, and function composition can be represented as the product of two matrices.
Consider a space \(V\) that is generated from vectors \(x_{1} and x_{2}\) that are each length 3 vectors. In addition, consider two functions, \(g\) and \(f\), that map a vector \(x\) from \(V\) onto a new space. We’ll assume the new space created by \(g(x)\) for all \(x\) is called \(G\) and the corresponsing space for \(f\) is called \(F\).
Then function composition results in
When we look at \( f[ g(x_{1}) ]\) we notice that \(g(x_{1})\) can be represented by its own set of vectors, say, \(y_{1}\) and \(y_{2}\). Then
and we can further assume \(f(x)\) also has its own set of vectors that can be used to represent it. Lets say these three vectors \(z_{1}, z_{2}, z_{3}\). Then
Again, our goal is to represent this vector using the three \(z\)s and so we want to group by \(z_{1}\),\(z_{2}\), and \(z_{3}\).
Lets use an inner product to rewrite the above,
where \begin{align} \beta_{:,1} = \begin{bmatrix} \beta^{1}{1}\ \beta^{1}{2} \end{bmatrix} ; ; ; G = \begin{bmatrix} \gamma^{1}{1} & \gamma^{2}{1} \ \gamma^{1}{2} & \gamma^{2}{2}\ \gamma^{1}{3} & \gamma^{2}{3}\ \end{bmatrix} \end{align}
In the same way, we would simplify \(\alpha_{2} f[ g(x_{2}) ]\) and find that
and so then
Finally, we define the product between the matrix \(G\) and the matrix \(B\) as
where \(B\) is a matrix that stacks, column-wise, the vectors \(\beta_{:,1}\) and \(\beta_{:,2}\)
Using the inner product, we can also write the matrix product between \(G\) and \(B\) as
\(AB = A A^{-1} = I\)#
Just like in univariate algebra there exists (or not) a function inverse, \(f^{-1}\), there can also exist a matrix inverse. A matrix \(B\) is an inverse to a matrix \(A\) is
We denote the matrix inverse to \(A\) as the matrix \(A^{-1}\).
Expectation algebra (Part 1)#
The expected value, and so the covariance and the variance, have many convienant properties due to the fact that the expected value is linear. i call these identities Expectation algebra.
First, we should make clear that the expected value is linear
Given a random variable \(X\) and constant \(\alpha \in \mathbb{R}\),
and so for X,Y r.v.s and \(\alpha,\beta \in \mathbb{R}\)
If we consider the vectors
then we can write
where \(\mathbb{E}(Z)\) applies the expectation to each random variable in the vector, or
Because the Covariance is also an expected value, it also has important properties.
Given a random variable \(X\), \(\text{Cov}(X,X) = \mathbb{E}[ (X-\mathbb{E}(X)) (X-\mathbb{E}(X)) ] = \mathbb{E}[(X-\mathbb{E}(X))^{2}] = \mathbb{V}(X)\).
The covariance is linear in its first and second position separately. In other words,
We can put these properties above all together to show that
We will further discuss below Expectation algebra that concerns the covariance when it is needed.
Normal + Normal = Normal#
Given two random variables, \(X \sim \mathcal{N}(\mu_{x}, \sigma_{x})\) and \(Y \sim \mathcal{N}(\mu_{y}, \sigma_{y})\), there is a highly unique property of the Normal: the sum of two (or more) Normal densities is itself another Normal distribution.
To be clear, the sum of two variables with the same density rarely, if ever, is distributed the same as the original two random variables. This property is so unique that sometimes the Normal density is defined as that density such that when summed is itself.
To be clear, the new normal density that is produced from a sum has a different expected value and variance than the original two random variables,
where \(\sigma_{xy} = \text{Cov}(X,Y)\).
This result is difficult to prove by definition because
However, we can note that a Normal density is completely defined from its expected value variance. That is, recall that
Lets use expedctation algebra for the random variable \(X+Y\) to find the expected value and variance.
First, the expected value:
and then the variance
This expected value and variance completely define \(X+Y\), and we can further generalize this ideas as
If we form the vector of constants \(\lambda = [\alpha, \beta]\) and the vector of random variables \(X = \begin{bmatrix}X_{1}\\X_{2}\end{bmatrix}\) then we could rewrite the above as
where \(\mu = [\mathbb{E}(X_{1}),\mathbb{E}(X_{2})]'\).
We can also take the same steps for the variance to find
The goal here is to write the above using the vector \(\lambda = [\alpha, \beta]\). The trick is to store the four covariances into a \(2 \times 2\) matrix
\(\Sigma\) contains all four covariances (remember the variance is a special type of covariance). This matrix \(\Sigma\) is called the covariance matrix.
Then we can write
These are the tools that we need to describe the Multivariate Normal as well as any linear transformation of a vector of Normally distributed random variables.
An inner product view point for the MVN#
Define \(\lambda \in \mathbb{R}^{D}\) be a \(1 \times D\) length vector and define \(X\) as a \(D \times 1\) length vector of random variables. Then \(X\) has a multivariate normal distribution if
This is a powerful viewpoint, and one that we can take even further by considering a “stack” of inner products.
Consider the matrix \(A_{m \times d}\) that is a stack of \(m\) vectors, \(\lambda_{1}, \lambda_{2}, \cdots, \lambda_{m}\). To be clear, we can write
Also consider the vector of random variables \(X = [X_{1},X_{2},\cdots,X_{d}]^{'}\). Then
Though \(\lamba x\) creates a single values, \(Ax\) creates a vector and so we need to discuss then a covariance matrix. Recall from above that the Coviariance matrix is a matrix where the entry in row \(i\) and column \(j\) is \(\text{Cov}(X_{i},X_{j})\).
In our example, the \((i,j)\) entry is then \(\text{Cov}(A_{i} X,A_{j}X)\) and we know that this equals \(A_{i} \text{Cov}(X,X)A_{j}^{'}\). Then the covariance matrix is
where we define \(\text{Cov}(X,X)\) as
Expectation algebra (Part 2)#
The part two of expedctation linear algebra concerns the Covariance matrix.
Given a vector of random variables, the covariance matrix defined as above is equal to
Then the familar covariance matrix for \(Ax\) falls out directly from this definition.
Marginal#
With the inner-product view of the MVN, determining the marginal density is trivial. Given a vector \(X \sim \text{MVN}(\mu, \Sigma)\) where \(\mathbb{E}(X) = \mu\) and the covariance matrix for \(X\) is \(\text{Cov}(X,X) = \Sigma\), the marginal density for \(X_{j}\) is \(\lambda X\) where the vector \(\lambda\) is a vector with the value one in the \(j^{\text{th}}\) position and the value zero otherwise. Then,
That is, for a multivariate normal. distribution the marginal density of \(X_{j}\), by integrating out all other random variables, is just a normal density with the jth mean and j,jth entry of the covariance matrix.
Let
Then \(X_{2} \sim \mathcal{N}( 2, 2 )\)
Conditional#
Though the marginal distribution for the MVN is easy to extarct, the conditional distribution of one or more random variables given one or more others is much more difficult to derive. First, we will write down the solution and then below derive this solution.
Given a vector of random variables \([X_{1}, X_{2}, \cdots X_{n}]^{T}\) that we partition into \(X = [X_{1}, X_{2}, \cdots X_{m}]\) and \(Y = [X_{m+1}, \cdots, X_{n}]\), the conditional distribution of \(X|Y\) is again a multivariate distribution
Derivation#
To derive the above result, we start with the probability density function for \([X,Y]^{T}\). Importantly, a distribution that is conditional on \(Y\) is one such that
In other words, the density of this conditional distribution only depends on \(X\). Our goal is to remove any expression that include \(y\) and recognize this density function as a familar distribution.
For the MVN, this means that
First, we removed the normalizing constant and defined \(\Lambda = \Sigma^{-1}\). The matrix \(\Lambda\) is often named the Precision matrix.
Our next steps will be to expand the expression in the exponential.
Lets continue to expand and then remove any items that are a constant with respect to X (ie do not include the variable X).
Now we are left to decide if the above function, up to a proportionality constant, looks like a density that we know. Intuitively, this almost looks like a MVN. Lets explore the MVN a bit to see if we can match the above function and traditional MVN.
MVN up to a constant#
The probability density over a vector \(X\) for an \(\text{MVN}(\mu, \Lambda)\) up to a constant looks like this
We see that, for a traditional multivariate normal density, the precision matrix is uniquely identified as the matrix that sits between \(X^{T}\) and \(X^{T}\). In other words, this term uniquely identifies the precision matrix
We also see that the mean vector, \(\mu\), is identified here
Lets apply this idea to find the precision matrix and mean vector for our conditional density. First, we see that for our conditional density we have a term \(X^{T} \Lambda_{XX}X\). That is, the precision matrix for our conditional density, which we will call \(\Lambda_{X|Y}\) is \(\Lambda_{X|Y} = \Lambda_{XX}\).
Now that we have the precision matrix, we can match up the remaining that accompanies \(2 X^{T}\) to find that
That is, the conditional distribution of \(X|Y\) does indeed have a multivariate normal density like
Though we found the conditional proability distribution of \(X\) given \(Y\), this is in terms of the precision matrix \(\Lambda\). Ideally, we would also want to be able to express this distribution in terms of \(\Sigma\) as well.
Our approach for defining the above in terms of \(\Sigma\) is to recognize that \(\Sigma^{-1} = \Lambda\). In other words, we need to find the inverse of \(\Lambda\) to define \(\Sigma\)
This identitiy implies that, for example,
and so then
and
import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets
from scipy.stats import multivariate_normal, norm
def plot_mvn_marginal_conditional(mu1=0.0, mu2=0.0, rho=0.0, x2_star=0.0):
sigma1 = 1.0
sigma2 = 1.0
mu = np.array([mu1, mu2])
Sigma = np.array([
[sigma1**2, rho * sigma1 * sigma2],
[rho * sigma1 * sigma2, sigma2**2]
])
x = np.linspace(-5, 5, 300)
y = np.linspace(-5, 5, 300)
X, Y = np.meshgrid(x, y)
pos = np.dstack((X, Y))
Z = multivariate_normal(mean=mu, cov=Sigma).pdf(pos)
# Marginals
fx = norm(loc=mu1, scale=sigma1).pdf(x)
fy = norm(loc=mu2, scale=sigma2).pdf(y)
# Conditional: X1 | X2 = x2_star
cond_mean = mu1 + rho * sigma1 / sigma2 * (x2_star - mu2)
cond_var = sigma1**2 * (1 - rho**2)
cond_sd = np.sqrt(cond_var)
f_cond = norm(loc=cond_mean, scale=cond_sd).pdf(x)
fig = plt.figure(figsize=(8, 7))
ax_joint = fig.add_axes([0.15, 0.15, 0.55, 0.55])
ax_x = fig.add_axes([0.15, 0.72, 0.55, 0.18], sharex=ax_joint)
ax_y = fig.add_axes([0.72, 0.15, 0.18, 0.55], sharey=ax_joint)
# Joint density contours
ax_joint.contour(X, Y, Z, levels=3, colors="black")
ax_joint.scatter(mu1, mu2, s=60, color="black")
ax_joint.axhline(0, linewidth=0.5)
ax_joint.axvline(0, linewidth=0.5)
# Show the conditioning value X2 = x2_star
ax_joint.axhline(x2_star, linestyle="--", linewidth=1.5)
ax_joint.set_xlim(-5, 5)
ax_joint.set_ylim(-5, 5)
ax_joint.set_aspect("equal")
ax_joint.set_xlabel(r"$X_1$")
ax_joint.set_ylabel(r"$X_2$")
# Marginal density of X1
ax_x.plot(x, fx, label=r"$f_{X_1}(x_1)$")
ax_x.plot(x, f_cond, linestyle="--", label=rf"$f_{{X_1 \mid X_2}}(x_1 \mid {x2_star:.2f})$")
ax_x.axvline(cond_mean, linestyle="--", linewidth=1)
ax_x.legend(fontsize=8)
ax_x.tick_params(labelbottom=False)
# Marginal density of X2
ax_y.plot(fy, y, label=r"$f_{X_2}(x_2)$")
ax_y.axhline(x2_star, linestyle="--", linewidth=1.5)
ax_y.tick_params(labelleft=False)
fig.suptitle(
rf"$X_1 \mid X_2={x2_star:.2f} \sim N({cond_mean:.2f}, {cond_var:.2f})$",
y=0.98
)
plt.show()
widgets.interact(
plot_mvn_marginal_conditional,
mu1=widgets.FloatSlider(value=0, min=-3, max=3, step=0.1, description=r"$\mu_1$"),
mu2=widgets.FloatSlider(value=0, min=-3, max=3, step=0.1, description=r"$\mu_2$"),
rho=widgets.FloatSlider(value=0, min=-0.95, max=0.95, step=0.05, description=r"$\rho$"),
x2_star=widgets.FloatSlider(value=0, min=-3, max=3, step=0.1, description=r"$x_2^*$")
);