Show code cell source Hide code cell source

import matplotlib.pyplot as plt
%matplotlib inline
import matplotlib_inline
matplotlib_inline.backend_inline.set_matplotlib_formats('svg')
import seaborn as sns
sns.set_context("paper")
sns.set_style("ticks");

Estimating Predictive Quantiles

Recall the definition of the \(q\)-quantile of a normal random variable we gave in Lecture 4. Recall also that you can say \(100q\)-th percentile instead of \(q\)-quantile. The exact definition holds for any random variable. Let’s repeat it here for convenience. Take \(X\) as a random variable and \(Y=g(X)\) as a function of \(X\). The \(q\)-quantile of \(Y\) is defined to be the number \(\mu_q\) for which:

\[ F(\mu_q) = p(Y \le \mu_q) = q, \]

where \(F(y)\) was defined to be the CDF of \(Y\). For example, the 0.50-quantile (also known as the median) is the value \(\mu_{0.50}\) for which:

\[ F(\mu_{0.50}) = p(Y \le \mu_{0.50}) = 0.5. \]

So, to find the quantiles, we need to 1) know the CDF of \(Y\) and 2) solve a root-finding problem with \(\mu_q\) as the unknown. We have already seen how one can estimate the CDF from samples, so we would only have to worry about the root-finding problem. It is not difficult to do, but since it is already implemented in numpy we will not bother with it. So, here is how you can find the empirical quantiles of \(Y\) for a specific example where \(g(x)\) is given by Example 3.4 of [Robert and Casella, 2004]:

\[ g(x) = \left(\cos(50x) + \sin(20x)\right)^2. \]

and \(X\sim U([0,1])\).

As usual, define the function and take some samples:

import numpy as np

# define the function here:
g = lambda x: (np.cos(50 * x) + np.sin(20 * x)) ** 2

# Maximum number of samples to take
max_n = 10000 
# Generate samples from X
x_samples = np.random.rand(max_n)
# Get the corresponding Y's
y_samples = g(x_samples)

Now let’s, find the 0.5-quantile:

mu_50 = np.quantile(y_samples, 0.5)
print(f"mu_50 = {mu_50:.2f}")

mu_50 = 0.60

# Let's find the 0.025-quantile
mu_025 = np.quantile(y_samples, 0.025)
print(f"mu_025 = {mu_025:.2f}")
# and the 0.975-quantile
mu_975 = np.quantile(y_samples, 0.975)
print(f"mu_975 = {mu_975:.2f}")

mu_025 = 0.00
mu_975 = 3.68

Let’s now mark these quantiles on the histogram of Y:

fig, ax = plt.subplots()
ax.hist(
    y_samples,
    density=True,
    alpha=0.25,
    bins=100
)
ax.plot(
    [mu_50],
    [0],
    "ro",
    markersize=5, 
    label="0.50-quantile (median)"
)
ax.plot(
    [mu_025],
    [0],
    "bo",
    markersize=5,
    label="0.025-quantile"
)
ax.plot(
    [mu_975],
    [0],
    "go",
    markersize=5,
    label="0.975-quantile"
)
ax.set_xlabel(r"$y$")
ax.set_ylabel(r"$p(y)$")
plt.legend(loc="best", frameon=False)
sns.despine(trim=False);

Very often, the predictive intervals are summarized using box plots:

fig, ax = plt.subplots()
# Here the input is percentiles not quantiles though
ax.boxplot(y_samples, whis=[2.5, 97.5], labels=["Y"]);
ax.set_ylabel(r"$y$")
sns.despine(trim=False);

In the plot above, the \(y\)-axis indicates possible values, the median is shown as an orange line, the box encapsulates 50% of the probability around the median, and the whiskers are extreme quantiles (here selected to be the 0.025 and 0.975 quantiles). Finally, the plot also shows the samples that fall outside the extreme quantiles.

An alternative to the box plot is the violin plot, which shows the distribution of the data (essentially a smoothed histogram). The violin plot is shown below.

# Write code that makes the violin plot using y_samples
fig, ax = plt.subplots()
ax.violinplot(y_samples, quantiles=[0.025, 0.5, 0.975])
ax.set_ylabel(r"$y$")
sns.despine(trim=False);

An alternative to

Questions

• How much probability do you have on the left of \(\mu_{0.50}\), i.e., what is \(p(Y \le \mu_{0.50})\)?

• How much probability do you have on the right of \(\mu_{0.50}\), i.e., what is \(p(Y \le \mu_{0.50})\)?

• How much probability do you have on the left of \(\mu_{0.025}\)?

• How much probability do you have on the right of \(\mu_{0.975}\)?

• How much probability do you have between \(\mu_{2.5}\) and \(\mu_{0.975}\)?

• The predictive quantiles are a very nice way to summarize the probability density of a random variable with a few numbers. For example, you can think of \(\mu_{0.50}\) as a central value of \(Y\). Often, we call the interval \([\mu_{0.025}, \mu_{0.975}]\) the 95% predictive interval. You can interpret this interval as \(Y \in [\mu_{0.025}, \mu_{0.975}]\) with 95% probability. Find a \(99\)% predictive interval for the \(Y\) of the example above.