Generating t-distributed random numbers: lesser numbers near zero.

According to a prescription given on Wikipedia, I tried generating Student's t-distributed random numbers with three degrees of freedom. I also generated these numbers using numpy's in-built random number generator for t-distribution. However, I feel that the numbers that I generated using the Wikipedia prescription don't align well with the analytical form of the distribution. The following figure shows five sample histograms for each case: manual (Wikipedia) generation (upper row) and python's in-built generation (lower row). I feel that around $x = 0$, there is some issue. Am I right in my observation? If yes, what is wrong with the way I am doing it? My python code is also given below..

Each histogram contains $10000$ numbers and the green curve represents the analytical distribution

import numpy as np
import matplotlib.pyplot as plt
plt.style.use("seaborn")
import seaborn as sns
from scipy.stats import t

"""Generate t distributed values"""
def f(x, mu):
n = len(x)
return np.sqrt(n) * (x.mean()-mu)/ x.std()

mu = 0
df = 3

for i in range(5):
plt.subplot(2,5,i+1)
t_vals = [f(np.random.normal(loc = mu, size = df + 1), mu) for i in range(10000)]
sns.distplot(t_vals, kde = False, norm_hist = True)
x = np.linspace(-5, 5, 100)
plt.plot(x, t.pdf(x, df))
plt.xlim([-5, 5])
plt.xlabel(r"$x$")
if i == 0:
plt.ylabel(r"$p(x)$")
if i == 2:
plt.title("Manually generated")

for i in range(5):
plt.subplot(2,5,i+6)
t_vals = np.random.standard_t(df, size = 10000)
sns.distplot(t_vals, kde = False, norm_hist = True)
x = np.linspace(-5, 5, 100)
plt.plot(x, t.pdf(x, df))
plt.xlim([-5, 5])
plt.xlabel(r"$x$")
if i == 0:
plt.ylabel(r"$p(x)$")
if i == 2:
plt.title("Generated using python")

plt.tight_layout()
plt.savefig("t_dists.pdf", bboxinches = "tight")

• There does not appear to be any appreciable difference between the theoretical density and the histograms you produce. Use (much) narrower bins to check.
– whuber
Mar 22, 2018 at 17:57
• You also seem to have picked wider bins for the histograms in the top row. Mar 22, 2018 at 19:57
• In both cases, I am using Freedman-Diaconis rule and so there is essentially no problem about binning. This is automated in seaborn. Mar 24, 2018 at 4:37
• @Peaceful, there is a problem with binning (they are not sufficiently narrow). If you wish to stick to that rule then generate more samples. Mar 29, 2018 at 8:51

Short version: the problem stands with NumPy x.std() which does not divide by the right degrees of freedom.

Repeating the experiment in R shows no discrepancy: either by comparing the histogram with the theoretical Student's $t$ density with three degrees of freedom

or the uniformity of the transform of the sample by the theoretical Student's $t$ cdf with three degrees of freedom

or the corresponding QQ-plot:

The sample of size 10⁵ was produced as follows in R:

X=matrix(rnorm(4*1e5),ncol=4)
Z=sqrt(4)*apply(X,1,mean)/apply(X,1,sd)


A Kolmogorov-Smirnov test also produces an acceptance of the null:

> ks.test(Z,"pt",df=3)

One-sample Kolmogorov-Smirnov test

data:  Z
D = 0.0039382, p-value = 0.08992
alternative hypothesis: two-sided


for one sample and

>  ks.test(Z,"pt",df=3)

One-sample Kolmogorov-Smirnov test

data:  Z
D = 0.0019529, p-value = 0.8402
alternative hypothesis: two-sided


for the next.

However..., the reason is much more mundane: it just happens that NumPy does not define the standard variance in the standard (Gosset's) way! Indeed it uses instead the root of $$\frac{1}{n}\sum_{i=1}^n (x_i-\bar{x})^2$$ which leads to a $t$ distribution inflated by $$\sqrt\frac{n}{n-1}$$ and hence to the observed discrepancy:

> ks.test(sqrt(4/3)*Z,"pt",df=3)

One-sample Kolmogorov-Smirnov test

data:  Z
D = 0.030732, p-value < 2.2e-16
alternative hypothesis: two-sided


While I have no personal objection to using $n$ instead of $n-1$ in the denominator, this definition clashes with Gosset's one and hence with the definition of the Student's $t$ distribution.

• I appreciate the efforts. Unfortunately, it does not provide answer to my question. I am showing five different realizations just to show that there seems to be a mismatch consistently. However, I sometimes do get good histograms, but less often. Your R code is right but does that point out error in my Python code (if any)? Binning is using Freedman-Diaconis rule in both the cases. Mar 24, 2018 at 4:41
• Exactly. And for the same reason, your answer doesn't help me if you say that you are getting a correct answer using other implementation. Mar 24, 2018 at 5:09