7.1. Expectations and joint distributions#

To motivate the considerations in this section, suppose that we have two random variables $X$ and $Y$ defined on a probability space $S$ and a real-valued function $g:{\mathbb{R}}^2\to \mathbb{R}$. We may then “transform” the random variables to obtain the random variable

where $(X,Y)$ denotes the $2$-dimensional random vector with $X$ and $Y$ as its components. Notice that the notation $g(X,Y)$ strictly interpreted is an abuse of notation; as the arrows above indicate, this notation actually stands for the composite $g\circ (X,Y)$, where we view $(X,Y)$ as a function in accordance with Definition 6.1.

In any case, provided that this “transformed” random variable $Z=g(X,Y)$ is continuous (say) with density $f_{g(X,Y)}(z)$, we may compute its expectation (according to the definition!) as

However, this formula is quite inconvenient to use in practice, due to the necessity of the density $f_{g(X,Y)}(z)$. We wonder if there is an alternate method to compute this expectation, one that uses the joint density $f(x,y)$. Indeed, there is! It is outlined in the following bivariate generalization of the LotUS from Theorem 3.1.

Though the argument in the discrete case is very similar to the one given for the univariate version of Theorem 3.1, we will not give it here. See if you can work it out on your own. You can imagine that the univariate and bivariate LotUS’s are special cases of a general multivariate LotUS that computes expecations of random variables of the form $g(X_1,\dots ,X_n)$, where each $X_i$ is a random variable and $n\ge 1$. I will leave you to imagine what the statement of this multivariate LotUS looks like. For those who might be interested, in the most general case, all of these LotUS’s are consequences of the general change-of-variables formula for Lebesgue integrals.

Our first application of the bivariate LotUS is to show that the expectation operator is multiplicative on independent random variables:

The proof is a simple computation using the LotUS; it appears in your homework for this chapter.

Our second application of the bivariate LotUS is to tie up a loose end from Section 3.9 and prove the full-strength version of “linearity of expectation.”

Let’s finish off the section by working through an example problem:

7.2. Expectations and conditional distributions#

In the previous section, we learned how to use joint distributions in computations of expectations. In this section, we define new types of expectations using conditional distributions.

In both cases, notice that conditional expected values depend on the particular choice of observed value $x$. This means that these conditional expectations yield functions:

where $A$ is the subset of $\mathbb{R}$ consisting of all $x$ for which either the conditional mass function $p(y|x)$ or density function $f(y|x)$ exists. There is a bit of notational awkwardness here, since I rather arbitrarily chose the name $h$ for this function. It would be much preferable to have notation that is somehow built from $E(Y\mid X=x)$, but there isn’t a commonly accepted one.

Let’s take a look at a practice problem before proceeding.

The input $x$ to the conditional expectation function

is an observation of a random variable, and thus this function is also “random.” To make this precise, let’s suppose for simplicity that $h$ is defined for all $x\in \mathbb{R}$, so that $h:\mathbb{R}\to \mathbb{R}$. Then, if $S$ is the probability space on which $X$ and $Y$ are defined, we may form the function

obtained as the composite $E(Y\mid X)\stackrel{\text{def}}{=}h\circ X$. Note that $E(Y\mid X)$ is a real-valued function defined on a probability space and thus (by definition!) it is a random variable. This is a bit confusing the first time you see it, so make sure to pause and ponder it for bit! To help, let’s take a look at a simple example problem where we unwind everything:

Now, since $Z=E(Y\mid X)$ is a random variable, it has an expected value. Provided that it is continuous (for example) with density function $f_{E(Y\mid X)}(z)$, we would compute its expectation (by definition!) as the integral

However, this formula is of little practical value, since we would need to compute the density $f_{E(Y\mid X)}(z)$ which can be quite difficult. So, we wonder: Is there an alternate way to compute this expectation? In fact there is! The key is to not compute the expectation according to the definition (7.3), but rather to use the LotUS:

However, if you push through the computations beginning with (7.4), you will find that the “iterated” expectation reduces to the expectation of $Y$. This is the content of:

Before we begin the proof, observe that there is a problem with the interpretation of the integral (7.4) since it extends over all $x$-values, whereas the function $E(Y\mid X=x)$ is only defined for those $x$-values for which the conditional density $f(y|x)$ exists (at least according to the conventions we use in this class). We ran into a similar problem when we discussed the Law of Total Probability in Theorem 6.5. However, we will fix this problem by declaring $E(Y\mid X=x)=0$ for those $x$-values outside the domain of the conditional expectation.

7.3. Density transformations#

It is very often the case that one needs to compute the density of a transformed continuous random variable. Actually, we saw such a situation already in Theorem 4.7 where we computed the distribution of an affine transformation of a normal random variable. We proved that theorem by explicitly computing the density of the transformed random variable. The main result in this section gives us a direct formula:

We will not prove the theorem, as it will end up being a special case of the generalized transformation theorem given in Theorem 7.6 below. But we should also say that there are simpler proofs of the theorem that do not use the heavy machinery that the proof of Theorem 7.6 relies upon; see, for example, Section 4.7 in [DBC21].

We are assuming that the support $T$ contains the range of the random variable $X$ so that the random variable $Y$ is defined; indeed, remember that $Y=r(X)$ is an abuse of notation that stands for the composite function $r\circ X$. (See the margin note directly to the right.)

Observe that the equality (7.5) is “one half” of what it means for $s$ to be the inverse of $r$. But in fact, since the domain of $s$ is precisely the range of $r$, as you may easily check we automatically have the other equality:

for all $y\in U$. Hence, $r$ is invertible with inverse function $s$.

Note that we are explicitly assuming that the inverse $s$ is continuously differentiable. But in other versions of the Density Transformation Theorem (and its generalization Theorem 7.6 below), certain conditions are placed on $r$ that guarantee differentiability of $s$ via the Inverse Function Theorem. We have elected to give the version of the theorem presented here because in practice (or, at least in many textbook problems) it is often quicker to simply check differentiability of $s$ directly, rather than verify the preconditions of the Inverse Function Theorem hold.

Let’s do some examples:

We now begin moving toward the generalization of The Density Transformation Theorem to higher dimensions. The statement of this generalization uses the following object:

Note that the gradient matrix is the transpose of the Jacobian matrix of $\mathbf{s}$ at $\mathbf{y}$, denoted $(\partial \mathbf{s}/\partial \mathbf{y})(\mathbf{y})$:

The Jacobian matrix is the matrix representation (with respect to the standard bases) of the derivative of $\mathbf{s}$, provided that the latter is differentiable. We will not talk about differentiability and derivatives in higher dimensions; for that, see Chapter 2 in [Spi65]. All we will say is that existence of the partial derivatives in the Jacobian matrix is necessary for differentiability, but not sufficient. If, however, the partial derivatives are all continuous, then $s$ is differentiable. (See Theorem 2-8 in [Spi65].)

The gradient matrix is used to define the affine tangent space approximation of $\mathbf{s}$ near a point of differentiability ${\mathbf{y}}^{\star }$ via the equation

Notice the similarity to the affine tangent line approximation studied in single-variable calculus.

With gradient matrices in hand, we now state the generalization of Theorem 7.5:

The following proof uses mathematics that are likely unfamiliar. Look through it if you like, but also feel free to skip it as well. I’ve included it really only out of principle, because I could not find a satisfactory proof in the standard references on my bookshelf.

7.4. Moment generating functions#

We begin with:

Take care to notice that I am not claiming that all of these moments exist for all random variables. Notice also that the first moment of $X$ is precisely its expectation, while its second central moment is its variance. The “higher moments” are more difficult to interpret. The situation with them is analogous to the “higher derivatives” of a function $y=f(x)$. I have good intuition for what the first two derivatives $f^{\prime }(x)$ and $f^{″}(x)$ measure, but I have much less intuition for what the thirty-first derivative $f^{(31)}(x)$ measures!

Actually, this analogy with derivatives can be carried further, so let’s leave the world of probability theory for a moment and return to calculus. Indeed, as you well know if a function $y=f(x)$ has derivatives of all orders at $x=0$, then we can form its Taylor series centered at $x=0$:

You also learned that this series may, or may not, converge to the original function $y=f(x)$ on an open interval about $x=0$. For example, the Taylor series of $y=e^x$ actually converges to $y=e^x$ everywhere:

On the other hand, there exist functions for which the Taylor series (7.6) exists and converges everywhere, but does not converge to the function on any open interval around $x=0$.

What this Uniqueness Theorem tells us is that complete knowledge of all the values $f^{(k)}(0)$ determines the function $y=f(x)$ uniquely, at least locally near $x=0$. Therefore, even though we don’t have good intuition for what all the higher derivatives $f^{(k)}(0)$ mean, they are still incredibly important and useful objects. Just think about it: If the Taylor series of $f(x)$ converges to $f(x)$ on an open interval containing $x=0$, then uncountably many functional values $f(x)$ are determined by a countable list of numbers $f(0),f^{\prime }(0),f^{″}(0),\dots$. There is an incredible amount of information encoded in these derivatives.

Now, hold this lesson in your mind for just a bit as we return to probability theory. We are about to see something very similar occur in the context of moments of random variables. The gadget that is going to play the role for random variables analogous to Taylor series is defined in:

The reason that the function $\psi (t)$ is said to generate the moments is encapsulated in:

Thus, the moments $E(X^k)$ may be extracted from the moment generating function $\psi (t)$ simply by taking derivatives and evaluating at $t=0$. In this sense, the function $\psi (t)$ generates the moments.

Now, the true power of moment generating functions comes from the following extremely important and useful theorem, which may be seen as an analog of the Taylor Series Uniqueness Theorem stated above as Theorem 7.7. It essentially says that: If you know all the moments, then you know the distribution.

To illustrate how this theorem may be used, we will establish the fundamental fact that sums of independent normal random variables are normal; this is stated officially in the second part of the following theorem.

Before turning to the worksheet to finish this section, it is worth extracting one of the steps used in the previous proof and putting it in its own theorem:

Now, let’s do an example:

7.5. Dependent random variables, covariance, and correlation#

In this section, we begin our study of dependent random variables, which are just random variables that are not independent. This study will continue through Section 9.3 in the next chapter, and then culminate in Chapter 11 where we learn how to use “networks” of dependent random variables to model real-world data.

While the definition of dependence of random variables $X$ and $Y$ is simply that they are not independent, it is helpful to conceptualize dependence as a flow of “information” or “influence” between them. If they are independent, then this flow vanishes and there is no transfer of “information”:

The most straightforward method to guarantee a transfer of “information” between the variables is to link them deterministically via a function $g:\mathbb{R}\to \mathbb{R}$, in the sense that $Y=g(X)$. This means that an observed value $x$ (of $X$) uniquely determines an observed value $y=g(x)$ (of $Y$). For random variables linked in this way, we prove what should be an intuitively obvious result:

What does a pair of functionally dependent random variables look like? For an example, let’s suppose that

Then, let’s simulate a draw of 1000 samples from $X$, toss them into

to obtain the associated $y$-values, and then produce a scatter plot:

Show code cell source Hide code cell source

import numpy as np
import scipy as sp
import seaborn as sns
import pandas as pd
import matplotlib.pyplot as plt
import scipy as sp
import math
import matplotlib_inline.backend_inline
import warnings
plt.style.use('../aux-files/custom_style_light.mplstyle')
matplotlib_inline.backend_inline.set_matplotlib_formats('svg')
warnings.filterwarnings('ignore')
blue = '#486AFB'
magenta = '#FD46FC'

def h(x):
    return x * (x - 1) * (x - 2)

np.random.seed(42)
x = sp.stats.norm.rvs(loc=1, scale=0.5, size=1000)
y = h(x)

sns.scatterplot(x=x, y=y)
plt.xlabel('$x$')
plt.ylabel('$y=g(x)$')
plt.ylim(-1.5, 1.5)
plt.gcf().set_size_inches(w=5, h=3)
plt.tight_layout()

The plot looks exactly like we would expect: A bunch of points lying on the graph of the function $y=g(x)$. However, very often with real-world data, an exact dependence $Y=g(X)$ does not truly hold. Instead, the functional relationship is “noisy”, resulting in scatter plots that look like this:

Show code cell source Hide code cell source

epsilon = sp.stats.norm.rvs(scale=0.15, size=1000)
grid = np.linspace(-0.5, 3)

_, ax = plt.subplots(ncols=2, figsize=(7, 3), sharey=True)

sns.scatterplot(x=x, y=y + epsilon, ax=ax[0])
ax[0].set_ylim(-1.5, 1.5)
ax[0].set_xlabel('$x$')
ax[0].set_ylabel('$y=g(x) + $noise')

sns.scatterplot(x=x, y=y + epsilon, alpha=0.2, ax=ax[1])
ax[1].plot(grid, h(grid), color='#FD46FC')
ax[1].set_xlabel('$x$')

plt.tight_layout()

The “noisy” functional relationship is drawn in the left-hand plot, while on the right-hand plot I have superimposed the graph of the function $y=g(x)$ for reference. Instead of lying directly on the graph of $y=g(x)$, the data is clustered along the graph.

The goal in this chapter is to study “noisy” linear dependencies between random variables; relationships that look like these:

Show code cell source Hide code cell source

grid = np.linspace(-2.5, 2.5)
epsilon = sp.stats.norm.rvs(scale=0.3, size=500)
m = [1, 0, -1]
x = sp.stats.norm.rvs(size=500)
_, ax = plt.subplots(ncols=3, figsize=(10, 3), sharey=True, sharex=True)

for i, m in enumerate(m):
    y = m * x + epsilon
    sns.scatterplot(x=x, y=y, ax=ax[i], alpha=0.3)
    ax[i].plot(grid, m * grid, color='#FD46FC')
    ax[i].set_xlim(-3, 3)
    ax[i].set_ylim(-3, 3)
    ax[i].set_xlabel('$x$')
    ax[i].set_ylabel('$y$')

plt.tight_layout()

Our goal in this section is to uncover ways to quantify or measure the strength of these types of “noisy” linear dependencies between random variables. We will discover that there are two such measures, called covariance and correlation. An alternate measure of more general dependence, called mutual information, will be studied in the next chapter in Section 9.3.

The definition of covariance is based on the following pair of basic observations:

Let $(x,y)$ be an observation of a two-dimensional random vector $(X,Y)$

1. If the observed values of $(X,Y)$ cluster along a line of positive slope, then $x$ and $y$ tend to be large (small) at the same time.

2. If the observed values of $(X,Y)$ cluster along a line of negative slope, then a large (small) value $x$ tends to be paired with a small (large) value $y$.

The visualizations that go along with these observations are:

In order to extract something useful from these observations, it is convenient to center the dataset by subtracting off the means:

Notice that

and similarly $E(Y-{\mu }_Y)=0$, so that when we carry out these replacements, we get random variables with mean $0$. Here’s an example:

Show code cell source Hide code cell source

_, axes = plt.subplots(ncols=2, figsize=(8, 3), sharex=True, sharey=True)
np.random.seed(42)
m = 0.5

x = sp.stats.norm.rvs(loc=2, scale=2, size=500)
y = m * x + epsilon + 2
mean = (x.mean(), y.mean())
sns.scatterplot(x=x, y=y, ax=axes[0], alpha=0.2)
axes[0].scatter(mean[0], mean[1], color='magenta', zorder=3, s=100)
axes[0].set_xlim(-6, 6)
axes[0].set_ylim(-6, 6)
axes[0].set_title('raw data')
axes[0].set_xlabel('$x$')
axes[0].set_ylabel('$y$')

sns.scatterplot(x=x - mean[0], y=y - mean[1], ax=axes[1], alpha=0.2)
axes[1].scatter(0, 0, color=magenta, zorder=3, s=100)
axes[1].set_xlim(-6, 6)
axes[1].set_ylim(-6, 6)
axes[1].set_title('centered data')
axes[1].set_xlabel('$x$')
axes[1].set_ylabel('$y$')
plt.tight_layout()

You can see that the dataset has not changed its shape—it has only shifted so that its mean (represented by the magenta dot) is at the origin $(0,0)$.

The reason that we center the data is that it allows us to conveniently rephrase our observations above in terms of signs:

Let $(x,y)$ be an observation of a centered two-dimensional random vector $(X,Y)$

1. If the observed values of $(X,Y)$ cluster along a line of positive slope, then $x$ and $y$ tend to have the same sign, i.e., $xy>0$.

2. If the observed values of $(X,Y)$ cluster along a line of negative slope, then $x$ and $y$ tend to have opposite signs, i.e., $xy<0$.

The new visualization is:

Essentially, the next definition takes the average value of the product $xy$, as $(x,y)$ ranges over observed pairs of values of a centered random vector $(X,Y)$. If this average value is positive, it suggests a (noisy) linear dependence with positive slope; if it is negative, it suggests a (noisy) linear dependence with negative slope. If the random variables are not centered, then we subtract off their means before computing the product and taking its average value.

Notice that the covariance of a random variable $X$ with itself is exactly its variance:

Before we look at an example, it will be convenient to state and prove the following generalization of Theorem 3.5:

Armed with this formula, let’s do an example problem:

A pair of very useful properties of covariance are listed in the following:

I suggest that you attempt to prove this theorem on your own. A special case appears in your homework for this chapter.

Bilinearity of covariance allows us to generalize Theorem 3.6 on the variance of an affine transformation of a random variable:

In particular, we see that if $\sigma (X_i,X_j)=0$ for all $i\ne j$ (i.e., if the random variables are pairwise uncorrelated; see below), the formula simplifies to:

We now turn toward the other measure of linear dependence, called correlation. To motivate this latter measure, we note that while the signs of covariances are significant, their precise numerical values may be less so, if we are attempting to use them to measure the strength of a linear dependence. One reason for this is that covariances are sensitive to the scales on which the variables are measured. For an example, consider the following two plots:

Show code cell source Hide code cell source

_, axes = plt.subplots(ncols=2, figsize=(8, 3))
np.random.seed(42)
m = 0.5

x = sp.stats.norm.rvs(loc=0, scale=2, size=500)
y = m * x + epsilon
mean = (x.mean(), y.mean())
sns.scatterplot(x=x, y=y, ax=axes[0], alpha=0.5)
axes[0].set_xlim(-5, 5)
axes[0].set_ylim(-5, 5)
axes[0].set_title('raw data')
axes[0].set_xlabel('$x$')
axes[0].set_ylabel('$y$')

sns.scatterplot(x=10 * x, y=10 * y, ax=axes[1], alpha=0.5)
axes[1].set_xlim(-50, 50)
axes[1].set_ylim(-50, 50)
axes[1].set_title('scaled data')
axes[1].set_xlabel('$x$')
axes[1].set_ylabel('$y$')
plt.tight_layout()

The only difference between the two datasets is the axis scales, but the linear relationship between the $x$’s and $y$’s remains the same. In fact, the scales on the axes differ by a factor of $10$, so by Theorem 7.14, we have

Thus, if we use the numerical value of the covariance to indicate the strength of a linear relationship, then we should conclude that the data on right-hand side is one hundred times more “linearly correlated” than the data on the left. But this is nonsense!

The remedy is to define a “normalized” measure of linear dependence:

The correlation ${\rho }_{XY}$ is also sometimes called the Pearson (product-moment) correlation coefficient, to distinguish it from other types of correlation measures. Its key properties are given in the following:

Remember, covariance and correlation were cooked up to measure linear dependencies between random variables. We wonder, then, what is the correlation between two random variables that are perfectly linearly dependent? Answer:

We give a name to two random variables whose correlation is zero:

Let’s take a look at an example before continuing:

You should think of independence as a strong form of non-correlation, and hence correlation is also a strong form of dependence. This is the content of the first part of the following result:

We have only considered a pair $(X,Y)$ of random variables; but what if we have a $d$-dimensional random vector $\mathbf{X}=(X_1,\dots ,X_d)$? Then we may compute the pairwise covariances $\sigma (X_i,X_j)$ and $\rho (X_i,X_j)$, and we can’t resist the urge to put them into matrices: