⇦ Back

Student’s t-Distributions are a generalisation of the standard normal distribution. Their exact shapes depend on the sample sizes - more correctly, the number of degrees of freedom (which is usually equal to the sample size minus one) - of the samples they model. The larger the number of degrees of freedom the more similar the t-distribution will become to the standard normal distribution, and at an infinite number of degrees of freedom they are the same. Student’s t-distributions are useful in three main areas of statistics:

1 Python Packages

The code on this page uses the numpy, matplotlib and scipy packages. These can be installed from the terminal with:

$ python3.11 -m pip install numpy
$ python3.11 -m pip install matplotlib
$ python3.11 -m pip install scipy

where python3.11 corresponds to the version of Python you have installed and are using.

Import these packages into Python as follows:

import numpy as np
from matplotlib import pyplot as plt
from scipy import stats as st

Generally, NumPy only lets you draw random numbers from a distribution while SciPy let’s you work with the PDF, CDF, PPF, etc, as exact mathematical formulas.

2 Probability Density Functions (PDFs)

The PDFs for t-distributions of varying degrees of freedom (dof) can be generated with SciPy:

# We want the distributions to appear on the plot from -3 to 3
x = np.linspace(-3, 3, 1000)
# The t-distribution for 1 degree of freedom
dof = 1
pdf_1 = st.t.pdf(x, dof)
# The t-distribution for 3 degrees of freedom
dof = 3
pdf_3 = st.t.pdf(x, dof)
# The t-distribution for 8 degrees of freedom
dof = 8
pdf_8 = st.t.pdf(x, dof)
# The t-distribution for 30 degrees of freedom
dof = 30
pdf_30 = st.t.pdf(x, dof)

For comparison, we will also look at the standard normal distribution. The expectation is that for larger sample sizes the t-distribution will be similar to this.

# We want the distributions to appear on the plot from -3 to 3
x = np.linspace(-3, 3, 1000)
# The standard normal distribution
pdf_norm = st.norm.pdf(x)

We can sample random numbers from a t-distribution using NumPy. Note that if you use a dof value less than 3 you will get values that don’t closely align with the PDFs because the mean and standard deviation of a t-distribution are defined differently for dof < 3.

# Set the seed so that we get the same random numbers each time this code runs
np.random.seed(20230811)
# Sample 1000 random numbers from the t-distribution with dof = 30
random_numbers = np.random.standard_t(30, 1000)

# As dof ≥ 3, the mean and standard deviation should be approximately 0 and 1
print(f'This should be close to zero: {random_numbers.mean()}')
print(f'This should be close to one: {random_numbers.std()}')
## This should be close to zero: 0.003504072170736343
## This should be close to one: 0.9948240325422479

Let’s take a look:

# Create plot
ax = plt.axes()
# Random numbers
label = '1000 Random Numbers'
ax.hist(random_numbers, density=True, label=label, color='gray', alpha=0.3)
# Probability density functions (PDFs)
ax.plot(x, pdf_norm, 'k', label='Standard Normal Distribution')
ax.plot(x, pdf_30, label='$t$-distribution, dof = 30')
ax.plot(x, pdf_8, label='$t$-distribution, dof = 8')
ax.plot(x, pdf_3, label='$t$-distribution, dof = 3')
ax.plot(x, pdf_1, label='$t$-distribution, dof = 1')
# Mean
ax.axvline(0, 0, st.norm.pdf(0) / 0.45, c='k', ls='--', label='Mean')
# Axes
ax.set_ylim([0, 0.45])
ax.set_xlim([-3, 3])
# Format
title = """Probability Density Functions (PDFs) of $t$-Distributions
and the Standard Normal Distribution"""
ax.set_title(title)
ax.set_ylabel('Relative Likelihood')
ax.set_xlabel('Distance of the Value of a Random Variable from the Mean')
ax.legend(frameon=False, fontsize='x-small')
plt.show()

Firstly, as the number of degrees of freedom increases from 1 to 3 to 8 to 30 we can indeed see that the t-distribution becomes more similar to the standard normal distribution (the solid black line). When dof = 30 they are almost indistinguishable (the blue and black lines are very close), hence there is a rule-of-thumb in statistics that says that for n > 30 you can use the normal distribution to model your data while for n < 30 you should use a t-distribution.

Secondly, our 1000 random numbers drawn from the dof = 30 t-distribution do indeed appear to follow the relevant distribution (the blue line).

3 Cumulative Distribution Functions (CDFs)

Much like the PDFs, SciPy lets us generate CDFs:

# We want the distributions to appear on the plot from -3 to 3
x = np.linspace(-3, 3, 1000)
# The t-distribution for 1 degree of freedom
dof = 1
cdf_1 = st.t.cdf(x, dof)
# The t-distribution for 3 degrees of freedom
dof = 3
cdf_3 = st.t.cdf(x, dof)
# The t-distribution for 8 degrees of freedom
dof = 8
cdf_8 = st.t.cdf(x, dof)
# The t-distribution for 30 degrees of freedom
dof = 30
cdf_30 = st.t.cdf(x, dof)

Let’s again compare these to the standard normal distribution:

# We want the distributions to appear on the plot from -3 to 3
x = np.linspace(-3, 3, 1000)
# The standard normal distribution
cdf_norm = st.norm.cdf(x)

We can also take the area under our random numbers:

# Cumulative sum of random numbers
random_numbers.sort()
auc = [np.trapz([0.1] * i, random_numbers[:i]) for i in range(300)]

Here’s the plot:

# Create plot
ax = plt.axes()
title = """Cumulative Distribution Functions (CDFs) of $t$-Distributions
and the Standard Normal Distribution"""
ax.set_title(title)
ax.set_ylabel('Cumulative Relative Likelihood')
ax.set_xlabel('Distance of the Value of a Random Variable from the Mean')
# Plot the cumulative sum of the random numbers
label = '1000 Random Numbers'
ax.hist(
    random_numbers, bins=40, density=True, cumulative=True,
    histtype='step', label=label, color='grey'
)
# Plot the CDFs
ax.plot(x, cdf_norm, 'k', label='Standard Normal Distribution')
ax.plot(x, cdf_30, label='$t$-distribution, dof = 30')
ax.plot(x, cdf_8, label='$t$-distribution, dof = 8')
ax.plot(x, cdf_3, label='$t$-distribution, dof = 3')
ax.plot(x, cdf_1, label='$t$-distribution, dof = 1')
# Plot the percent point function (PPF)
label = 'Percent point function (PPF)'
quantile = 0.4
ax.axvline(st.t.ppf(quantile, 30), 0, quantile, c='k', ls='--', label=label)
ax.axhline(quantile, 0, (st.t.ppf(quantile, 30) - -3) / 6, c='k', ls='--')
# Axes
ax.set_ylim([0, 1])
ax.set_xlim([-3, 3])
# Finish
ax.legend(frameon=False, fontsize='x-small')
plt.show()

You’ll notice that this plot contains a ‘percent point function’. This is the inverse of the CDF - for a given relative likelihood it will return the distance from the mean (in standard deviations) below which a random variable taken from the distribution in question will be with this likelihood. In the plot above, a relative likelihood of 0.4 (ie a quantile of 0.4 out of 1) corresponds to a distance of -0.26 standard deviations from the mean.

3.1 Aside: Limits of Agreement

On this Wikipedia page about the limits of agreement when calculating inter-rater reliability it states that they can be approximated using 1.96 standard deviations for sample sizes greater than 60 while 2 standard deviations should be used for sample sizes less than or equal to 60. The reason that \(n=60\) is the cut-off for using 2 standard deviations is because this is almost the last t-distribution where a relative likelihood of 0.025 (which is used for 95% limits of agreement) corresponds to a value with a magnitude greater than 2:

print(st.t.ppf(0.025, 59))
print(st.t.ppf(0.025, 60))
print(st.t.ppf(0.025, 61))
## -2.0009953770482105
## -2.000297821058262
## -1.9996235841149783

As you can see, \(n=61\) (ie \(dof=60\)) is the last t-distribution where 0.025 corresponds to more than 2 standard deviations of distance from the mean.

3.2 Aside: Prediction Intervals and Confidence Intervals

For normally distributed variables, the 95% prediction interval (95% PI) can be calculated as:

\[ 95\%\ PI = mean \pm t_{0.975,\ n-1} \times \sqrt{\frac{n+1}{n}}\times SD \]

where \(t_{0.975,\ n−1}\) is the 97.5% quantile of a Student’s t-distribution with \(n-1\) degrees of freedom. Here’s a partial example:

# Sample size
n = 20
# Critical t-statistic
t_crit = st.t.ppf(q=0.975, df=n - 1)
print(t_crit)
## 2.093024054408263

Another option is to use the interval() function with a confidence of 0.95:

# Sample size
n = 20
# Critical t-statistics as an interval
t_crit_interval = st.t.interval(confidence=0.95, df=n - 1)
print(t_crit_interval)
## (-2.093024054408263, 2.093024054408263)

Similar calculations can be done when estimating 95% confidence intervals.

4 More Information

For more info on Student’s t-distributions and how to work with them in Python, see the Wikipedia page, SciPy documentation and NumPy documentation.

⇦ Back