Exploring non-parametric statistics

import numpy as np
import scipy
from matplotlib import pyplot as plt
import seaborn as sns

As a follow-on blog post to our recent podcast about non-parametric systems, we created a visual reference guide to illustrate the discussion we had. In the following examples, we are using Python 3.9 and the excellent scientific computing libraries of numpy, scipy, matplotlib, and seaborn.

Distributions

As discussed on the podcast, which we recommend you listen to prior to reading this post, probability distributions are, simply, mathematical functions that represent the probability of a value occurring. In the podcast we discussed several standard distributions, and how they are often used. Below are illustrations of the common distributions we mentioned:

Discrete Uniform:is the simplest distribution. It is used when all outcomes are equally likely. For example, if you roll a fair die the probability of getting any number is 1/6. Thus the probability of rolling a 1, 2, 3, 4, 5, or 6 is uniform.

This can be formalized as: $f(x|n) = \frac{1}{n}$ where $n$ is the number of possible outcomes.

# Show an example of a uniform distribution
num_samples = 1000
data = np.random.choice([1,2,3,4,5,6], size = num_samples)

sns.countplot(x=data, palette='mako')
plt.axhline(num_samples/6.0, color='red', linestyle='dashed', linewidth=2, label='Expected Value')
plt.title('Dice Rolls')
plt.xlabel('Dice Value')
plt.ylabel('Count')
plt.ylim(0, num_samples)


plt.legend(loc="upper left")
plt.savefig('dice_rolls.png', dpi=300)
plt.show()

Normal

Also called the Gaussian distribution. This is the most common distribution we discuss in statistics. It is used when the data is symmetric and shaped like a "bell". For example, the heights of people in a population are normally distributed as well as their SAT scores.

This distribution is defined by two parameters: the mean and the standard deviation. The mean $(μ)$, typically 0 is the center of the distribution, and the standard deviation $(σ)$, typically 1, is a measure of the spread of the distribution.

The Probability Distribution Function (PDF) of the normal distribution is more formally defined as:

$f(x|\mu,\sigma) = \frac{1}{\sigma\sqrt{2\pi}}e^{\left( \frac{-(x-\mu)^2}{2\sigma^2}\right)}
$

Below is an example of 2022 of the SAT score data with the red dotted line being a 'perfect' normal distribution.

# Show an example of a normal distribution using SAT scores
num_samples = 1000
mean = 1050 # mean SAT score in 2022
std = 216 # std of SAT scores in 2022
data = np.random.normal(mean, std, size = num_samples)

# show the distribution
ax = sns.displot(data, bins=40, color='skyblue', stat='density', label='SAT Scores')
# show the expected values
x_pdf = np.linspace(min(data), max(data), 100)
y_pdf = scipy.stats.norm.pdf(x_pdf, mean, std)
plt.plot(x_pdf, y_pdf, color='red', linestyle='dashed', label='Expected Dist')
plt.legend()
plt.title('SAT Scores Distribution')
plt.xlabel('Score')

plt.savefig('SAT_scores.png', dpi=300)

‍

Source: nces.ed.gov/programs/digest/d22/tables/dt22_226.40.asp

In the plot below, we overlay both the expected distribution of SAT (red dotted line) vs actual distribution (green dotted line) on top of the underlying data: a very close fit!

# Show an example of a normal distribution using SAT scores
# set numpy seed
np.random.seed(10)
num_samples = 1000
mean = 1050 # mean SAT score in 2022
std = 216 # std of SAT scores in 2022
data = np.random.normal(mean, std, size = num_samples)

# show the distribution
ax = sns.displot(data, bins=40, color='skyblue', stat='density', label='Observed SAT Scores', alpha=0.2,linewidth=0.1)
# show the expected values
x_pdf = np.linspace(min(data), max(data), 100)
y_pdf = scipy.stats.norm.pdf(x_pdf, mean, std)
plt.plot(x_pdf, y_pdf, color='red', linestyle='dashed', label='Expected Distribution')
# show the interpreted values
y_pdf = scipy.stats.norm.pdf(x_pdf, np.mean(data), np.std(data))
plt.plot(x_pdf, y_pdf, color='green', linestyle='dashed', label='Observed Distribution')

plt.legend(loc="lower left")
plt.title("Expected vs Observed SAT Scores")
plt.savefig('SAT_scores_dists.png', dpi=300)

Assumptions Baked into Parametric Statistics

Switching gears, below, we will attempt to model the powerlaw distribution using a normal distribution. The powerlaw distribution is used to model the distribution of wealth, the frequency of words in a language, the size of cities, and many other phenomena. It can be formalized as $f(x|\alpha,x_{\min}) = \frac{\alpha-1}{x_{\min}} \left(\frac{x}{x_{\min}}\right)^{-\alpha}$ where $\alpha$ is the scaling factor and $x_{\min}$ is the minimum value.

# Using normal distribution to model a powerlaw distribution
from scipy.stats import powerlaw
samples = 100000
a = 0.659
# powerlaw distribution
data = powerlaw.rvs(a, loc=0, size=samples)
ax = sns.displot(data, color='skyblue', stat='density', alpha=0.4,linewidth=0.1)
# show the expected distribution
x_pdf = np.linspace(powerlaw.ppf(0.01, a), powerlaw.ppf(0.99, a), 100)
y_pdf = scipy.stats.powerlaw.pdf(x_pdf, a=0.659)
plt.plot(x_pdf, y_pdf, color='blue', linestyle='dashed', label='Powerlaw Distribution')
# using normal distribution
y_pdf = scipy.stats.norm.pdf(x_pdf, np.mean(data), np.std(data))
plt.plot(x_pdf, y_pdf, color='green', linestyle='dashed', label='Fitted Normal Distribution')


plt.legend(loc="upper right");
plt.savefig('normal_for_powerlaw.png', dpi=300)

In the plot above you can see that unlike the SAT distributions above, these two lines are no where near each other. Different distributions used for different use cases.

‍

Tests Using Parametric Statistics

In the following section we will briefly touch on concepts of statistical testing. In the example below we calculate and illustrate the probability that an SAT score of 1450 occurred by chance.

By using the probability distribution we can calculate the probability of observing a given value. For example, if we assume that the data is normally distributed we can calculate the probability of obtaining an SAT score over 1400. This is shown with the shaded area in the figure where green is the probability of the outcome.

With a little bit of math, we can use the estimated parameters to calculate the probability of observing a given value or a more extreme value. This is called a p-value. The p-value is calculated using the z-score, which is the number of standard deviations $(σ)$ away from the mean the observed value is. -

$$z = {x- \mu \over \sigma}$$

If the p-value is less than a threshold (usually $\alpha = 0.05$) we say that the result is statistically significant. This means that the probability of observing the data we did is less than 5% assuming the null hypothesis is true. The null hypothesis is the hypothesis that there is no difference between the groups we are comparing.

# Probability that an SAT score of 1450 occured by chance
from scipy.stats import norm
# create a normal distribution
mean = 1050
std = 216
x_pdf = np.linspace(400, 1600, 100)
y_pdf = scipy.stats.norm.pdf(x_pdf, mean, std)
plt.plot(x_pdf, y_pdf, color='red', linestyle='dashed', label='Normal Distribution')
# create a point
point = 1400
# show the point
plt.scatter(point, 0, color='black')
# show the probability
plt.fill_between(x_pdf, 0, y_pdf, where=(x_pdf >= point), color='green', alpha=0.5, label='Probability')
# numerical probability
prob = scipy.stats.norm.cdf(point, mean, std)
print('Probability of getting a score of 1400 or higher is {:.2f}%'.format(100 - prob*100))
# z-score of the point
z = (point - mean)/std
print('The z-score of the point is {:.2f} σ'.format(z))
# the p_value of the point
p_value = scipy.stats.norm.sf(z)
print('The p-value of the point is {:.2f}'.format(p_value))
# show the probability
plt.legend(loc="upper left");
plt.savefig('point_probability.png', dpi=300)

Probability of getting a score of 1400 or higher is 5.26%
The z-score of the point is 1.62 σ
The p-value of the point is 0.05

‍

Non-Parametric Statistics

Tests Using Non-Parametric Statistics

Kolmogorov–Smirnov test

Let's revisit our example comparing the powerlaw distribution with the gaussian distribution. This time we will find the mean $μ$ and std $σ$ of the powerlaw and then generate data that fits those parameters. We will then compare the powerlaw data with the Gaussian data using the non-parametric Kolmogorov–Smirnov test.

from scipy.stats import kstest

samples = 100000

a = 0.659
# powerlaw distribution
data_powerlaw = powerlaw.rvs(a, loc=0, size=samples)
# normal distribution
data_normal = np.random.normal(np.mean(data_powerlaw), np.std(data_powerlaw), samples)

# Plot the data
sns.histplot(data_powerlaw, bins=100, color='skyblue', stat='density', alpha = 0.3)
sns.histplot(data_normal, bins=100, color='red', stat='density', alpha = 0.3)
#sns.kdeplot(data, color = 'blue', alpha = 0, linestyle='dashed', label='Powerlaw Distribution')

# Plot the normal distribution
x_pdf = np.linspace(min(data), max(data), 100)
y_pdf = scipy.stats.norm.pdf(x_pdf, np.mean(data), np.std(data))
plt.plot(x_pdf, y_pdf, color='red', linestyle='dashed', label='Normal Distribution')

# Plot the powerlaw distribution
x_pdf = np.linspace(powerlaw.ppf(0.01, a), powerlaw.ppf(0.99, a), 100)
y_pdf = scipy.stats.powerlaw.pdf(x_pdf, a=0.659)
plt.plot(x_pdf, y_pdf, color='blue', linestyle='dashed', label='Expected Distribution')

plt.xlim(0,1)
plt.legend(loc="upper right");
plt.title("Normal vs. Powerlaw Distribution")

# Perform the KS test
ks_statistic, p_value = kstest(data_powerlaw, data_normal)

# Print the KS statistic and p-value
print("KS Statistic:", ks_statistic, "\np-value:", p_value)

plt.savefig('KS_test.png', dpi=300)

KS Statistic: 0.09362
p-value: 0.0

The Kolmogorov–Smirnov Test checks if a sample is part of a probability distribution, with no assumptions made about the distribution, in the same manner we tested with a z-score above, which assumes normal distribution and that the population standard deviation is known, as well as a sample size, $n>30$. The K-S statistic calculates the distance between the two empirical distribution functions to determine if there is a statistically significant difference between the two. The Kolmogorov–Smirnov Test has a null hypothesis that the data follows the same distribution, with the alternative hypothesis being that the data does not follow the same distribution. The test statistic is

$$D = \max _{1 \leq  i \leq  N}(F(Y_i)-\frac{i-1}{N},\frac{i}{N}- F(Y_i))$$

• $H_0$: the distributions are the same (Null Hypothesis)
• $H_a$: the distributions are different

In the above case we find that the p-value is $\approx0.0$, which means that we can reject the null hypothesis that the data is normally distributed. Since the K-S test is essentially assumption free, it can be used in instances were Z-score would normally be used. The only real negative of using non-parametric test statistics, such as the K-S test is if your parametric assumptions are met, as in the Z-score example above, the parametric tests will produce higher statistical power than a non-parametric test (power essentially meaning finding statistical significance if there is any to be found).

Conclusion

In this post we've provided a visual aid extension to our non-parametric podcast discussion. If this content was helpful, you would like us to go deeper in an area, or if you have questions or requests, please submit feedback on our site.

Until next time,

Andrew & Sid, The AI Fundamentalists.
‍

‍