Implicit Statistical Validation¶
You collect data in real-time and build models to analyze this data. Your decisions depend on this data. But most statistical theory relies on data honoring parametric assumptions; assumptions that the data (or more correctly the residuals) are normally distributed. Second that this data is stationary and not drifting in its characteristics. Because if it is drifting, the assumptions for when the model was built and validated may not stand true at the time of scoring.
Instead of assuming the de facto validity and stability of the parameters, it is critical that you build-in safety measures to actively monitor, alert, and trigger any violations from these assumptions. Below we build a few of these implicit tests. We will demonstrate conducting three tests --
- First that the observations are unbiased and sampled from a normal distribution.
- Second that the attributes/features of observation are discriminant enough for inclusion into the model.
- Third, assess the distribution to best describe the observations if they do not follow the normal distribution.
Normality Validation¶
Central limit theorem is a pivotal statistical theorem that confers a level of assurance that most random unbiased samples tend to aggregate (sum) into a normal curve (with a specific mean and standard deviation) and that as the number of samples increase, the distribution approximates normality ever more closely.</p>
Given this observation, if we believe our samples (data measures) are truly from an unbiased sample, it is safe to assume our data is normal. There are two tests we can do easily if the normality assumption holds true. They are --
- Is the data "normal" -- a nice Gaussian distribution?
- Is the data "stationary" over time -- meaning has it not drifted over time?
Let us accomplish a test for both.
Testing the "bell" shape¶
The test for normality, in simple terms, is a test to see if the curve looks like a nice symmetric bell with a defined peak and a curvature. We will create synthetic datasets, one normal and another non-normal, to test.
from scipy import stats
# Create the normal series of 1000 numbers, around the mean (peak) of 10 and a standard deviation of 5 (spread)
n_series = np.random.normal(10, 5, 1000)
# Create 100 numbers randomly
non_n_series = stats.uniform.rvs(size=1000) * 20
Let us plot the two synthetic distributions to visually gauge the curvature and peak of the bells.
%matplotlib inline
from scipy import stats, integrate
# Plot the two series
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 3))
sns.distplot(n_series, ax=ax1)
# Visually, you should already see central limit theorem at 100 values (the kernel desnity)
# almost looking like a normal distribution
sns.distplot(non_n_series, ax=ax2);
Is the data "normal"?¶
Spot a bell curve?
from scipy import stats
# Test if the X is normal; centered around zero and nicely distributed.
def test_norm(x):
# Conduct Kolmogorov-Smirnov test for normality
(kstat, pvalue) = stats.mstats.normaltest(x)
return pvalue > 0.05
# Plot two side-by-side plots -- one normal and one non-normal
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 3))
sns.distplot(n_series, ax=ax1)
# Put the results of our investigation into normality as plot title
ax1.set_title('''The known normal series is tested to be {0}normal'''.format(
"" if test_norm(n_series) else "not "))
sns.distplot(non_n_series, ax=ax2)
# Put the results of our investigation into non-normality as plot title
ax2.set_title('''The known non-normal series is tested to be {0}normal'''.
format("" if test_norm(non_n_series) else "not "))
plt.show();
Is the data "stationary"?¶
We also want to spot if a series that is presumably stationary, meaning the peak and the shape that is relatively static, still remains so over time. Because we do not have the luxury of time (in this notebook), let us create two more more synthetic normal series.
- One normal series has the same standard deviation (of 5) but different mean as our prior good distribution -- meaning the expanse/curvature of the distribution has not changed too substantively.
- The other series has the same mean (of 10) but different standard deviation to show the curvature has changed but peak remains at the same mean.
# Make a copy of our known prior series
known_n_series = n_series
# Create a new series with same mean of 10 but wider standard deviation
n_series_with_diff_stddev = np.random.normal(10, 8, 1000)
# Create another new series with same standard deviation but mean higher than previous
n_series_with_diff_mean = np.random.normal(15, 5, 1000)
# Plot the series to visually show the series; share X axis to we can compare mean and variance
f, (ax1, ax2, ax3) = plt.subplots(3, 1, sharex=True, figsize=(12, 6))
# Original series is green
sns.distplot(known_n_series, color="g", hist=False, ax=ax1)
ax1.set_title('Known normal series')
# Wider series is yellow
sns.distplot(n_series_with_diff_stddev, color="y", hist=False, ax=ax2)
ax2.set_title('A new series with wider std deviation but same mean')
# Shifted series is blue
sns.distplot(n_series_with_diff_mean, color="b", hist=False, ax=ax3)
ax3.set_title('A new series with same std deviation but different mean');
Has the series shifted its "center"/mean?¶
# Do two series have same mean?
def test_mean_drift(x, y):
(tstat, pvalue) = stats.ttest_ind(x - x.mean(), y - x.mean())
print "The series have" + (" different "
if pvalue <= 0.05 else " same ") + "mean"
return pvalue <= 0.05
#Test drift when series has not moved; can have different variance
test_mean_drift(known_n_series, n_series_with_diff_stddev)
#Test drift when series has moved 5 unit
test_mean_drift(known_n_series, n_series_with_diff_mean);
Has the series' variance/expanse changed?¶
# Do two series have same variance?
def test_var_drift(x, y):
#Use levene test to test variance; also called homoscedasticity
(tstat, pvalue) = stats.levene(x - x.mean(), y - y.mean())
print "The series have" + (" different " if pvalue <= 0.05 else " same "
) + "standard deviation"
return pvalue <= 0.05
# Test variance drift when series has not moved
test_var_drift(known_n_series, n_series_with_diff_mean)
# Test variance drift when series has moved
test_var_drift(known_n_series, n_series_with_diff_stddev);
Feature Evaluation¶
Given two normal series, it is relatively easy to spot the stationarity or significant differences. However, when you have multiple normal series, the permutorial combinations for significant difference tests can be substantively computational and error prone unless variance (scatter distance) across all the series are normalized.
Consider the following example of a wine dataset. There are three kinds/classes of wines each with different (or we believe) levels of malic acid, ash, magnesium, phenols, color hue etc. In order to determine if malic acid is a substantive feature that distinguishes the class of wine, we need to compare 3 distributions, wine 1 with wine 2, wine 2 with wine 3, and wine 1 with wine 3. And we need to repeat the exercise for every feature like acidity, hue etc before tossing it into the model.
Underneath, we show we may use pairwise tests to check for significant differences between series. This test isn't so much the stationarity test although similar theory is at play. It is used to identify significant attributes and differences relevant for the model (model that signifies the class of wine in our example).
Based on a real wine dataset (http://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.names), analyze if there is substantive difference in alcohol and malic acid levels of three kinds of wine.
# Download the wine dataset
# Initialize the data model that one would expect from the dataset
col_names = [
'Class', 'Alcohol', 'Malic acid', 'Ash', 'Alcalinity of ash', 'Magnesium',
'Total phenols', 'Flavanoids', 'Nonflavanoid phenols', 'Proanthocyanins',
'Color intensity', 'Hue', 'OD280_OD315', 'Proline'
]
# Fetch the data
wines = pd.read_csv(
'http://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data',
names=col_names,
index_col=False)
# Change the class of the wine to a categorical type from numeric type
wines['WineClass'] = wines.Class.apply(lambda x: 'Wine ' + str(x))
# Display a preview
pd_display(wines, 'Wines dataset from UCI machine learning website')
Group Differences¶
Can we visually, along each feature attribute, see if the means and standard deviations are different across the three wine groups?
# Find all "continuous/numerical features" that can predict the class of wine
all_measures = wines[[x for x in wines.select_dtypes(include=['number']).columns.values if 'Class' not in x]]
# What are we plotting means against
class_feature = 'WineClass'
# Plot the violin plots to spot feature differences
fig, axes = plt.subplots(nrows=len(all_measures.columns)/4 + 1, ncols=4, sharex=True, figsize=(12,12))
for i, col in enumerate(all_measures.columns):
sns.violinplot(data=wines, y=col, x='WineClass', ax=axes[i/4, i%4], linewidth=1)
sns.despine()
plt.tight_layout()
Do the wines differ?¶
Just visually looking at the malic acid and ash plots, the means and standard deviations of wine 1 and wine 2 for malic acid; and wine 1 and wine 3 for ash do not seem too far off from each other. The plots and these tests, besides telling us what individual class differences, can be useful in spotting the useful features to build a wine classification automaton if we were to keep the model parsimonious (aka feature selection).
Let us also compute (besides visual plotting before) if the means and standard deviations are different across the three wine groups along other features. If the features are distinctly different along different wine classes, that means the feature is discriminatory enough to be included into the automaton system.
# Are the feature/per-feature measurement "clumps" similar between wines
# We will conduct Tukey test to find pairwise differences. Tukey test is similar to ttest
# But since we have multiple groups (wine 1, 2, and 3 leading to 1-2, 2-3, 1-3 comparisons)
# it is better to use tukey test (or bonferroni test) since stddev is normalized across the groups
from statsmodels.stats.multicomp import pairwise_tukeyhsd
# Along each feature, compare the wines and see if the mean/variance is significantly different
for feature in all_measures:
compare_result = pairwise_tukeyhsd(wines[feature], wines[class_feature])
display(
HTML(
'<span class="{2} text-center simpletable"><caption class="caption">When comparing <u><i>{0}</i></u>, we conclude {4} groups are {1}</caption><div class="row center text-center col-lg-12">{3}</div></span>'.
format(feature, 'different' if any(compare_result.reject) else
'same', 'text-primary' if any(compare_result.reject) else
'text-muted', compare_result._results_table.as_html(), np.count_nonzero(compare_result.reject))))
It so turns out all features of the wine are distinct enough between the various wine classes.
Distribution Checks¶
Data is not always "clean". Perhaps there is a stochastic process that drifts into an anomalous state; or perhaps data is just spotty and the NaN practices lead to wrongful values in the data. Either case, missing data and wrongful data can cause our normality assumptions to be easily violated. Different parametric/distribution assumptions must be applied if we determine that the normality is invalid. Below, we spot missing values, outlier values, and best distribution to describe incoming data.
We will use the air quality dataset at https://archive.ics.uci.edu/ml/datasets/Air+Quality. The datasets are chemical readings at various time intervals.
import requests, zipfile, StringIO
# Let us read the air quality dataset from UCI -- https://archive.ics.uci.edu/ml/datasets/Air+Quality
uri = 'https://archive.ics.uci.edu/ml/machine-learning-databases/00360/AirQualityUCI.zip'
# Read the data from the web
air_quality = pd.read_excel(zipfile.ZipFile(StringIO.StringIO(requests.get(uri).content)).open('AirQualityUCI.xlsx'))
# Display a preview of the data
pd_display(air_quality, 'Air quality readings of ~ 1 year between 2004-2005')
Time Series¶
The data is temporal in nature with continuous readings every hour. But the Pandas DataFrame treats the day and time of recording distinctly. In order to plot continuous maps (and study the trends), it is important to turn this data into a "time-series".
# Create time index
air_quality_ts = air_quality.copy()
# Convert date and time into string forms
air_quality_ts['Date'] = air_quality_ts.Date.astype('str')
air_quality_ts['Time'] = air_quality_ts.Time.astype('str')
# Create a Datetime column
air_quality_ts['datetime'] = air_quality_ts.apply(lambda x : pd.to_datetime(x.Date + ' ' + x.Time), axis=1)
# Change the index into a time series
air_quality_ts.set_index('datetime', inplace=True)
# Drop Date and Time since we already have the index
air_quality_ts.drop(['Date', 'Time'], axis=1, inplace=True)
# Convert everything to numeric datatype
air_quality_ts = air_quality_ts.astype('float')
# Display a preview
pd_display(air_quality_ts, 'Air quality data rendered as time series.')
Visual "outliers"¶
Let us spot the trends visually to spot any outlier data. As is evident in the tabular display above, there are some readings at -200.0 leading us to suspect sensor readings. Plotting may reveal if the distribution is indeed "normal".
%matplotlib inline
# Plot the time series
sns.set(color_codes=True)
# Create chart area
fig, ax = plt.subplots(2, 1, figsize=(12, 8))
# Plot the CO(GT) column -- Hourly averaged concentration CO
series = air_quality_ts['CO(GT)']
# Plot both trend line
series.plot(ax=ax[0], linewidth=0.25)
# And the kernel density distribution
sns.kdeplot(series, shade=True, ax=ax[1], linewidth=0.5);
Plot Observations¶
The plot seems to suggest that the distribution is bipolar -- the normal range mostly lies between [0, 10] and abnormal range is indicated with a -200.
Separating the outliers¶
Let us use simple three sigma measure to separate the outliers. If you see the https://en.wikipedia.org/wiki/68%E2%80%9395%E2%80%9399.7_rule, it suggests that most of all normal values lie within two standard deviations of the mean. Anything outside these "normal" distributions is perhaps an outlier. Let us spot where the outliers exist.
from scipy import stats
# Convert our original values into a dataframe
co_df = series.to_frame()
# Compute a column where the (absolute) deviation of value is beyond two standard deviations away
co_df['is_outlier'] = (np.abs(stats.zscore(co_df['CO(GT)'])) > 2)
# Show abnormal/outlier values
pd_display(co_df[co_df.is_outlier], "Points where the sensor recorded an outlier value")
If we exclude the outlier values and treat the observations, are the observations normal?
# Also plot "normal" values for visually cleaner dataset
good_values = co_df[~co_df.is_outlier]['CO(GT)']
ax = sns.distplot(good_values, rug=True, hist=True);
# Check for normality of the good values to make sure we may now do parametric tests
ax.set_title("The dataset is {0} normal".format("NOT" if stats.mstats.normaltest(good_values)[1] < 0.05 else ""));
Findings¶
It seems to suggest from (preview of the dataframe) that the sensor is misbehaving at 4AM in the morning usually. Something to consider?
Moreover, it looks the data is not really "normally" distributed -- the distinct bell shape is not apparent. It seems to be skewed. What distribution is this sensor exhibiting?
# Create a function to run ks test against all known distributions
from math import log
def run_dist_test(x, cdf):
#fit our data set against every probability distribution
dist = getattr(stats, cdf)
#estimate the parameters for known samples -- aka mean, variance, skew etc
params = dist.fit(good_values)
# apply the Kolmogorov Smirnov test
(D, p) = stats.kstest(x, cdf, args=params)
# if the distribution is significant return the D values else pass null
return log(p) if not np.isnan(p) and p <= 0.05 and p > 0 else None
# Get all known distributions from the web; for demo we will only check 5 distributions
cdfs = pd.read_html('https://docs.scipy.org/doc/scipy-0.16.1/reference/stats.html')[1][0].to_frame('cdf').head()
# Evaluate the distribution
cdfs['log_pvalue'] = cdfs['cdf'].apply(lambda cdf: run_dist_test(good_values, cdf))
# Display the preview of the cdf evaluations
pd_display(cdfs.dropna().sort_values('log_pvalue'), 'Log likelihood of the values following a known distribution -- smaller is better')
Alpha Distribution¶
From the tests of all known distributions, the observed sensor values seem to follow the "alpha" distribution the best. See more about alpha distribution at (http://www.itl.nist.gov/div898/software/dataplot/refman2/auxillar/alppdf.htm). Here is how the theoretical distribution and the observed distributions look side-by-side.
# Plot the observed values
sns.distplot(good_values, rug=False, hist=True);
# Plot the theoretical values
Image('https://docs.scipy.org/doc/scipy-0.16.1/reference/_images/scipy-stats-alpha-1.png')
Conclusions¶
Do not just leave the implicit assumptions of your model go unchecked. The data could be drifting, missing, over selecting (features) or changing its inherent shape over time. Build data-awareness in the models.