Chapter 5: Applications to Machine Learning

Introduction

Across the four preceding chapters, we learned the fundamentals of descriptive statistics and probability (Chapter 1), probability distributions (Chapter 2), statistical estimation and hypothesis testing (Chapter 3), and Bayesian statistics (Chapter 4). As the final chapter of this series, this chapter examines how this statistical knowledge is actually used inside concrete machine learning algorithms, through hands-on implementation.

Many machine learning algorithms are built directly on top of statistical concepts. Linear regression is itself a statistical optimization technique called the least squares method, and logistic regression is trained within the statistical estimation framework known as maximum likelihood estimation. The Naive Bayes classifier is a direct application of Bayes' theorem — which we learned in Chapter 1 — to classification problems, and the Gaussian process extends the ideas of probability distributions and Bayesian statistics to quantify uncertainty in regression problems. Furthermore, model evaluation for comparing trained models' performance and A/B testing for validating the effect of interventions both directly apply the statistical testing ideas we learned in Chapter 3.

Through this chapter, I hope you'll come to feel that "learning statistics" and "mastering machine learning" are one continuous path.

The table below summarizes the algorithms covered in this chapter along with their statistical foundations.

1. The Statistical Interpretation of Linear Regression and Least Squares

1.1 What Is the Least Squares Method

The least squares method (Ordinary Least Squares: OLS) determines a model's parameters by minimizing the sum of squared residuals — the differences between observed and predicted values. A simple regression model is expressed as follows.

$$y_i = \beta_0 + \beta_1 x_i + \varepsilon_i$$

Here, $\beta_0$ is the intercept, $\beta_1$ is the slope, and $\varepsilon_i$ is the error term. The least squares method finds the $\beta_0, \beta_1$ that minimize the following residual sum of squares (RSS).

$$\text{RSS}(\beta_0, \beta_1) = \sum_{i=1}^{n}(y_i - \beta_0 - \beta_1 x_i)^2$$

This optimization problem can be solved analytically, and scikit-learn's LinearRegression internally uses this analytical solution (or an efficient numerical computation) to estimate the parameters.

1.2 Statistical Assumptions and Uncertainty in the Coefficients

To treat linear regression not merely as an "optimization" but as "statistical inference," we place several assumptions on the error term $\varepsilon_i$.

• Linearity: the relationship between $y$ and $x$ can be expressed with linear parameters
• Independence: the errors $\varepsilon_i$ are mutually independent
• Homoscedasticity: the variance of the error is constant regardless of the value of $x$
• Normality: the errors follow a normal distribution (needed for constructing hypothesis tests and confidence intervals)

Under these assumptions, the least squares estimator is known to be unbiased and to have the minimum variance among all linear unbiased estimators (the Gauss–Markov theorem). If we further assume normality, we can use the standard error $\text{SE}(\hat{\beta}_1)$ of the coefficient $\hat{\beta}_1$ to construct a t-test or confidence interval, in the same framework we learned in Chapter 3.

The test statistic for the null hypothesis "the slope $\beta_1$ is 0" is given by:

$$t = \frac{\hat{\beta}_1}{\text{SE}(\hat{\beta}_1)}$$

This $t$ follows a t-distribution with $n-p-1$ degrees of freedom (the sample size $n$ minus the number of parameters $p+1$, including the intercept). If the p-value is small, we can conclude that the explanatory variable $x$ contributes to predicting the response variable $y$ in a statistically significant way.

1.3 Python Implementation

Let's implement a simple regression model that predicts an exam score from study time, including the statistical interpretation.

import numpy as np
from scipy import stats
from sklearn.linear_model import LinearRegression

# Generate experimental data: study hours and exam scores (points)
np.random.seed(42)
n = 30
study_hours = np.random.uniform(1, 10, n)
true_intercept, true_slope = 50, 4.5
noise = np.random.normal(0, 6, n)
scores = true_intercept + true_slope * study_hours + noise

X = study_hours.reshape(-1, 1)
y = scores

# Fit with scikit-learn using the least squares method
model = LinearRegression()
model.fit(X, y)

intercept = model.intercept_
slope = model.coef_[0]
y_pred = model.predict(X)

# Coefficient of determination R^2
r_squared = model.score(X, y)

# Residual analysis
residuals = y - y_pred
n_samples, n_features = X.shape
dof = n_samples - n_features - 1  # Degrees of freedom (2 parameters estimated: intercept and slope)
residual_std_error = np.sqrt(np.sum(residuals**2) / dof)

# Standard error of the coefficients (standard formula using the design matrix)
X_design = np.column_stack([np.ones(n_samples), X])
cov_matrix = residual_std_error**2 * np.linalg.inv(X_design.T @ X_design)
se_intercept = np.sqrt(cov_matrix[0, 0])
se_slope = np.sqrt(cov_matrix[1, 1])

# 95% confidence interval and t-test for the slope (null hypothesis: slope = 0)
t_critical = stats.t.ppf(0.975, dof)
ci_slope = (slope - t_critical * se_slope, slope + t_critical * se_slope)
t_stat_slope = slope / se_slope
p_value_slope = 2 * (1 - stats.t.cdf(abs(t_stat_slope), dof))

print(f"Intercept: {intercept:.4f}")
print(f"Slope: {slope:.4f}")
print(f"Coefficient of determination R^2: {r_squared:.4f}")
print(f"Residual standard error: {residual_std_error:.4f}")
print(f"\nStandard error of slope: {se_slope:.4f}")
print(f"95% confidence interval of slope: [{ci_slope[0]:.4f}, {ci_slope[1]:.4f}]")
print(f"t-statistic of slope: {t_stat_slope:.4f}")
print(f"p-value of slope: {p_value_slope:.2e}")

Execution Result:

Intercept: 53.5025
Slope: 3.5713
Coefficient of determination R^2: 0.7671
Residual standard error: 5.0899

Standard error of slope: 0.3719
95% confidence interval of slope: [2.8095, 4.3330]
t-statistic of slope: 9.6036
p-value of slope: 2.33e-10

2. Logistic Regression and Maximum Likelihood Estimation

2.1 Why Linear Regression Falls Short

If we directly apply linear regression to a binary classification problem where the response variable is, say, "pass/fail," the predicted values can fall below 0 or exceed 1, making them impossible to interpret as probabilities. Logistic regression solves this by passing the linear combination through a sigmoid function, which confines the output to the range $[0,1]$.

$$p_i = P(y_i=1|x_i) = \frac{1}{1+e^{-(\beta_0+\beta_1 x_i)}}$$

2.2 Parameter Estimation via Maximum Likelihood Estimation

The parameters of logistic regression are found not by the least squares method but by maximum likelihood estimation (MLE). Assuming the response variable $y_i$ follows a Bernoulli distribution, the likelihood function can be written as follows.

$$L(\beta) = \prod_{i=1}^{n} p_i^{y_i}(1-p_i)^{1-y_i}$$

Taking the logarithm gives us the more tractable log-likelihood.

$$\ell(\beta) = \sum_{i=1}^{n}\left[y_i \log p_i + (1-y_i)\log(1-p_i)\right]$$

There is no analytical solution to the problem of finding the $\beta$ that maximizes this log-likelihood, so scikit-learn's LogisticRegression uses gradient-based numerical optimization (the L-BFGS method by default) to estimate $\beta$ iteratively. The negative of this log-likelihood is called log loss (or cross-entropy loss), and it serves as the loss function for many classification models.

2.3 Interpreting the Coefficients: Log-Odds and Odds Ratios

The coefficient $\beta_1$ in logistic regression can be interpreted as the change in log-odds for a one-unit increase in the explanatory variable $x$. Taking $\exp(\beta_1)$ gives the odds ratio, which can be interpreted intuitively as "how many times larger the odds of success become when $x$ increases by one unit."

2.4 Python Implementation

import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import log_loss

# Generate experimental data: a binary classification problem predicting pass/fail from study hours
np.random.seed(42)
n = 100
study_hours = np.random.uniform(0, 10, n)
# True logistic model: P(pass=1) = sigmoid(-4 + 1.0 * study_hours)
true_logit = -4 + 1.0 * study_hours
true_prob = 1 / (1 + np.exp(-true_logit))
pass_fail = (np.random.uniform(0, 1, n) < true_prob).astype(int)

X = study_hours.reshape(-1, 1)
y = pass_fail

# Logistic regression via maximum likelihood estimation with scikit-learn
# C=1e6 effectively disables regularization, bringing this closer to pure MLE
model = LogisticRegression(C=1e6)
model.fit(X, y)

intercept = model.intercept_[0]
coef = model.coef_[0][0]

# Predicted probabilities and log-likelihood
proba = model.predict_proba(X)[:, 1]
# log_loss is the average negative log-likelihood, so multiplying by n and flipping
# the sign recovers the log-likelihood
log_likelihood = -log_loss(y, proba, normalize=True) * n

accuracy = model.score(X, y)

print(f"Intercept: {intercept:.4f}")
print(f"Coefficient: {coef:.4f}")
print(f"Log-Likelihood: {log_likelihood:.4f}")
print(f"Training accuracy: {accuracy:.4f}")

# Interpreting the coefficient: how many times the odds increase per extra hour of study
odds_ratio = np.exp(coef)
print(f"\nOdds ratio exp(coefficient): {odds_ratio:.4f}")
print(f"-> Each additional hour of study multiplies the odds of passing by about {odds_ratio:.2f}")

# Predict the probability of passing for specific study hours
test_hours = np.array([[2], [4], [6], [8]])
test_proba = model.predict_proba(test_hours)[:, 1]
print("\nPredicted pass probability by study hours:")
for h, p in zip(test_hours.ravel(), test_proba):
    print(f"  {h:.0f} hours: {p:.4f}")

Execution Result:

Intercept: -5.2407
Coefficient: 1.2589
Log-Likelihood: -24.4132
Training accuracy: 0.8800

Odds ratio exp(coefficient): 3.5214
-> Each additional hour of study multiplies the odds of passing by about 3.52

Predicted pass probability by study hours:
  2 hours: 0.0616
  4 hours: 0.4489
  6 hours: 0.9099
  8 hours: 0.9921

3. Implementing a Naive Bayes Classifier

3.1 Applying Bayes' Theorem to Classification

Let's recall Bayes' theorem from Chapter 1. Given features $x_1,\ldots,x_n$, the posterior probability of class $y$ can be written as:

$$P(y|x_1,\ldots,x_n) = \frac{P(y)\,P(x_1,\ldots,x_n|y)}{P(x_1,\ldots,x_n)}$$

The problem here is that as the number of features grows, the amount of data needed to estimate the joint distribution $P(x_1,\ldots,x_n|y)$ grows exponentially. The Naive Bayes classifier sidesteps this problem by making a simplifying (naive) assumption: "given the class $y$, each feature is conditionally independent of the others."

$$P(x_1,\ldots,x_n|y) \approx \prod_{i=1}^{n}P(x_i|y)$$

Under this assumption, classification reduces to comparing the class-wise posterior probabilities — or, more precisely, the product of the prior probability and the likelihood, since the denominator $P(x_1,\ldots,x_n)$ is the same for every class and can be ignored — and selecting the largest one.

$$\hat{y} = \arg\max_{y}\; P(y)\prod_{i=1}^{n}P(x_i|y)$$

3.2 Gaussian Naive Bayes

When the features are continuous, Gaussian Naive Bayes assumes that the distribution of each feature within each class follows a normal distribution. It estimates the mean $\mu_{y,i}$ and variance $\sigma_{y,i}^2$ of each feature for each class $y$ from the training data, and substitutes them into the normal probability density function to compute the likelihood $P(x_i|y)$.

3.3 Python Implementation

We'll implement this using the Iris dataset (measurements of iris flowers, a 3-class classification problem), which is real-world data.

import numpy as np
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB
from sklearn.metrics import accuracy_score

# Use the Iris dataset (iris flower measurements, 3-class classification)
iris = load_iris()
X, y = iris.data, iris.target
class_names = iris.target_names

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=42, stratify=y
)

# Train a Gaussian Naive Bayes classifier
# Assumes each feature follows a normal distribution within each class, and that features are independent
gnb = GaussianNB()
gnb.fit(X_train, y_train)

y_pred = gnb.predict(X_test)
accuracy = accuracy_score(y_test, y_pred)

print(f"Test accuracy: {accuracy:.4f}")

# Learned class-wise prior probabilities (the proportion of each class in the training data)
print("\nPrior probability P(class) for each class:")
for name, prior in zip(class_names, gnb.class_prior_):
    print(f"  {name}: {prior:.4f}")

# Class-wise mean and variance (Gaussian parameters) for the first feature (sepal length)
print("\nMean and variance by class for feature 0 (sepal length):")
for i, name in enumerate(class_names):
    mean = gnb.theta_[i, 0]
    var = gnb.var_[i, 0]
    print(f"  {name}: mean={mean:.4f}, variance={var:.4f}")

# Check the posterior probability of belonging to each class for one test sample
sample = X_test[0].reshape(1, -1)
posterior = gnb.predict_proba(sample)[0]
true_label = class_names[y_test[0]]

print(f"\nFeatures of one test sample: {X_test[0]}")
print(f"True class: {true_label}")
print("Posterior probability P(class|features) for each class:")
for name, p in zip(class_names, posterior):
    print(f"  {name}: {p:.4f}")

Execution Result:

Test accuracy: 0.9111

Prior probability P(class) for each class:
  setosa: 0.3333
  versicolor: 0.3333
  virginica: 0.3333

Mean and variance by class for feature 0 (sepal length):
  setosa: mean=4.9886, variance=0.1033
  versicolor: mean=5.9486, variance=0.2408
  virginica: mean=6.6829, variance=0.4248

Features of one test sample: [7.3 2.9 6.3 1.8]
True class: virginica
Posterior probability P(class|features) for each class:
  setosa: 0.0000
  versicolor: 0.0000
  virginica: 1.0000

4. Quantifying Prediction Uncertainty with Gaussian Processes

4.1 When a Point Estimate Alone Is Not Enough

The linear and logistic regression models we've covered so far are essentially point estimate models — they predict a single "most likely" value. In practice, however, there are many situations where information about uncertainty — "how confident can we be in this prediction?" — matters a great deal. For example, in regions with little training data, or for inputs that lie far outside the distribution of the training data, the reliability of a prediction should naturally be lower.

A Gaussian process (GP) is a regression method that defines a probability distribution over functions themselves, allowing it to naturally output not only the mean of a prediction but also its uncertainty (variance).

4.2 The Idea Behind Gaussian Processes

A Gaussian process is defined as a stochastic process with the property that, for any finite set of input points $x_1,\ldots,x_m$, the function values $f(x_1),\ldots,f(x_m)$ always follow a multivariate normal distribution.

$$f(x) \sim \mathcal{GP}(m(x), k(x, x'))$$

Here $m(x)$ is the mean function (often taken to be 0), and $k(x,x')$ is a covariance function called the kernel function, which expresses how "similar" two points $x$ and $x'$ are. In this chapter we use the most widely used kernel, the RBF kernel (Radial Basis Function kernel). The kernel encodes a smoothness assumption: nearby inputs are strongly correlated, while distant inputs are weakly correlated.

Conditioning on the observed data $D$, the predictive distribution at a new input point $x_*$ is also normal, with its mean giving the point estimate and its standard deviation representing the prediction's uncertainty. A key feature of Gaussian processes is that the standard deviation is small (high confidence) near the training data and grows larger (lower confidence) further away from it.

4.3 Python Implementation

We'll assume a scenario with sensor data at only a few observation points, and check how the prediction uncertainty changes near versus far from those observations.

import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel

# Simulate a situation with only a few observation points (imagine sensor inspection data)
np.random.seed(42)
X_train = np.array([[1.0], [2.0], [3.5], [6.0], [7.5]])
def true_function(x):
    return np.sin(x) * 2 + 0.3 * x

y_train = true_function(X_train).ravel() + np.random.normal(0, 0.15, X_train.shape[0])

# Kernel: RBF kernel (captures smoothness) + WhiteKernel (captures observation noise)
# We constrain the search range moderately so hyperparameters stay in a reasonable
# range even with so few data points
kernel = RBF(length_scale=1.5, length_scale_bounds=(0.5, 5.0)) \
    + WhiteKernel(noise_level=0.05, noise_level_bounds=(1e-3, 0.5))

gpr = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10, random_state=42)
gpr.fit(X_train, y_train)

print(f"Learned kernel parameters: {gpr.kernel_}")
print(f"Log marginal likelihood: {gpr.log_marginal_likelihood(gpr.kernel_.theta):.4f}")

# Compare prediction uncertainty (standard deviation) near observations versus far from them
X_test = np.array([[1.5], [4.5], [7.0], [12.0]])
y_mean, y_std = gpr.predict(X_test, return_std=True)

print("\nMean and standard deviation (uncertainty) for each prediction point:")
for x, m, s in zip(X_test.ravel(), y_mean, y_std):
    print(f"  x={x:.1f}: mean prediction={m:.4f}, std dev={s:.4f}, "
          f"95% prediction interval=[{m - 1.96*s:.4f}, {m + 1.96*s:.4f}]")

print("\nPoints near observed data (x=1.5, 4.5) have small standard deviations, while")
print("a point extrapolated far outside the observed range (x=12.0) shows a larger one.")

Execution Result:

Learned kernel parameters: RBF(length_scale=1.53) + WhiteKernel(noise_level=0.5)
Log marginal likelihood: -13.1404

Mean and standard deviation (uncertainty) for each prediction point:
  x=1.5: mean prediction=1.8016, std dev=0.8472, 95% prediction interval=[0.1412, 3.4621]
  x=4.5: mean prediction=0.3935, std dev=0.9498, 95% prediction interval=[-1.4681, 2.2551]
  x=7.0: mean prediction=2.5160, std dev=0.8721, 95% prediction interval=[0.8067, 4.2253]
  x=12.0: mean prediction=0.0369, std dev=1.2247, 95% prediction interval=[-2.3635, 2.4373]

Points near observed data (x=1.5, 4.5) have small standard deviations, while
a point extrapolated far outside the observed range (x=12.0) shows a larger one.

5. Model Evaluation and Statistical Testing

5.1 Why a Single Evaluation Is Not Enough

If we split the data into training and test sets only once to evaluate a model, the evaluation result can be swayed by chance depending on how the split happened to fall (e.g., an unusually large share of easy samples ended up in the test set). Cross validation addresses this by splitting the data into $k$ folds and repeating the process of using one fold for testing and the rest for training $k$ times. This lets us assess a model's performance including the spread (standard deviation) of the evaluation results.

5.2 Statistically Testing the Performance Difference Between Models

When comparing the cross-validation scores of two models, simply looking at which mean is larger doesn't tell us whether the difference reflects genuine, statistically meaningful signal or just random variation. This is where the hypothesis testing ideas from Chapter 3 come in handy. Because the scores of two models evaluated on the same folds are paired, we use a paired t-test to test the null hypothesis that the mean score difference is 0.

5.3 Python Implementation

We'll compare the performance of ordinary linear regression against regularized Ridge regression using 10-fold cross validation, and test whether the difference is statistically significant.

import numpy as np
from scipy import stats
from sklearn.datasets import make_regression
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.model_selection import KFold, cross_val_score

# Synthetic regression data (with some correlation among features)
X, y = make_regression(
    n_samples=200, n_features=15, n_informative=8,
    noise=15.0, random_state=42
)

# Compare ordinary linear regression against regularized Ridge regression
model_a = LinearRegression()
model_b = Ridge(alpha=10.0)

# Use the same split (KFold) for both models to make this a paired comparison
kf = KFold(n_splits=10, shuffle=True, random_state=42)

scores_a = cross_val_score(model_a, X, y, cv=kf, scoring="neg_mean_squared_error")
scores_b = cross_val_score(model_b, X, y, cv=kf, scoring="neg_mean_squared_error")

mse_a = -scores_a
mse_b = -scores_b

print("=== MSE (Mean Squared Error) from 10-fold cross validation ===")
print(f"Linear Regression: mean MSE={mse_a.mean():.4f}, std dev={mse_a.std():.4f}")
print(f"Ridge Regression  : mean MSE={mse_b.mean():.4f}, std dev={mse_b.std():.4f}")

# Paired t-test: tests whether the score difference between the two models,
# evaluated on the same folds, is statistically significant (null hypothesis: mean difference is 0)
t_stat, p_value = stats.ttest_rel(mse_a, mse_b)

print(f"\n=== Paired t-test (testing the MSE difference between models) ===")
print(f"t-statistic: {t_stat:.4f}")
print(f"p-value: {p_value:.4f}")

alpha = 0.05
if p_value < alpha:
    print(f"Since p-value < {alpha}, the performance difference between the two models is statistically significant")
else:
    print(f"Since p-value >= {alpha}, the performance difference between the two models is not statistically significant")

Execution Result:

=== MSE (Mean Squared Error) from 10-fold cross validation ===
Linear Regression: mean MSE=244.2298, std dev=87.0243
Ridge Regression  : mean MSE=357.3553, std dev=96.6626

=== Paired t-test (testing the MSE difference between models) ===
t-statistic: -2.7153
p-value: 0.0238
Since p-value < 0.05, the performance difference between the two models is statistically significant

6. Statistical Methods for A/B Testing

6.1 The Frequentist Approach to A/B Testing

In Chapter 4, we handled A/B testing within the Bayesian framework, directly computing "the probability that Variant B is better than Variant A." In this section, by contrast, we implement a frequentist A/B test using the long-established two-sample t-test. It's a widely used method for comparing a continuous metric (purchase amount, time on page, etc.) between two groups.

When we can't assume the two groups have equal variance, the robust choice in practice is Welch's t-test. The test statistic is computed as follows.

$$t = \frac{\bar{x}_A - \bar{x}_B}{\sqrt{s_A^2/n_A + s_B^2/n_B}}$$

Here $\bar{x}_A, \bar{x}_B$ are the sample means of each group, $s_A^2, s_B^2$ are the sample variances of each group, and $n_A, n_B$ are the sample sizes of each group. The degrees of freedom are computed using an approximation called the Welch–Satterthwaite equation.

6.2 Effect Size: What Statistical Significance Alone Doesn't Tell You

The p-value helps us judge whether a difference is likely due to chance, but it doesn't tell us how large the difference actually is in practical terms. When the sample size is very large, even a difference that is practically negligible can become statistically significant. This is why we also use effect size, and in particular Cohen's d, to evaluate the substantive magnitude of a difference.

$$d = \frac{\bar{x}_B - \bar{x}_A}{s_{\text{pooled}}}$$

As a general rule of thumb, $|d|\approx 0.2$ indicates a small effect, $0.5$ a medium effect, and $0.8$ or above a large effect (though this is only a rough guideline, and interpretation varies by field and context).

6.3 Python Implementation

We'll use Welch's t-test to examine whether a new layout (Variant B) for an e-commerce site increases the average purchase amount. We implement the standard formula by hand and cross-check it against the result from scipy.stats.

import numpy as np
from scipy import stats

# A/B test: examine whether a new e-commerce layout (Variant B) increases the average purchase amount
# Assume we observed the purchase amounts (in yen) of buyers under Variant A (existing layout) and Variant B (new layout)
np.random.seed(42)
n_a, n_b = 60, 65
group_a = np.random.normal(loc=4200, scale=900, size=n_a)   # Variant A: mean 4200 yen
group_b = np.random.normal(loc=4550, scale=950, size=n_b)   # Variant B: mean 4550 yen

def welch_t_test(x, y):
    """
    Manually implement Welch's t-test (a two-sample t-test that does not assume equal variance)

    Parameters:
    -----------
    x, y : ndarray
        The two groups of data to compare

    Returns:
    --------
    t_stat, dof, p_value : t-statistic, degrees of freedom, two-sided p-value
    """
    mean_x, mean_y = np.mean(x), np.mean(y)
    var_x, var_y = np.var(x, ddof=1), np.var(y, ddof=1)
    n_x, n_y = len(x), len(y)

    se_diff = np.sqrt(var_x / n_x + var_y / n_y)
    t_stat = (mean_x - mean_y) / se_diff

    # Approximate the degrees of freedom using the Welch-Satterthwaite equation
    dof = (var_x / n_x + var_y / n_y) ** 2 / (
        (var_x / n_x) ** 2 / (n_x - 1) + (var_y / n_y) ** 2 / (n_y - 1)
    )

    p_value = 2 * (1 - stats.t.cdf(abs(t_stat), dof))
    return t_stat, dof, p_value

t_stat, dof, p_value = welch_t_test(group_a, group_b)

print("=== Welch's t-test (manual implementation) ===")
print(f"Variant A: mean={group_a.mean():.2f} yen, std dev={group_a.std(ddof=1):.2f} yen, n={n_a}")
print(f"Variant B: mean={group_b.mean():.2f} yen, std dev={group_b.std(ddof=1):.2f} yen, n={n_b}")
print(f"\nt-statistic: {t_stat:.4f}")
print(f"Degrees of freedom (approx.): {dof:.2f}")
print(f"p-value (two-sided): {p_value:.4f}")

# Cross-check with scipy.stats.ttest_ind (equal_var=False gives Welch's test)
t_stat_scipy, p_value_scipy = stats.ttest_ind(group_a, group_b, equal_var=False)
print(f"\n=== Cross-check with scipy.stats.ttest_ind ===")
print(f"t-statistic: {t_stat_scipy:.4f}")
print(f"p-value: {p_value_scipy:.4f}")

# Effect size (Cohen's d): evaluates the substantive magnitude of the difference, not just its significance
pooled_std = np.sqrt(((n_a - 1) * group_a.var(ddof=1) + (n_b - 1) * group_b.var(ddof=1)) / (n_a + n_b - 2))
cohens_d = (group_b.mean() - group_a.mean()) / pooled_std

print(f"\nEffect size (Cohen's d): {cohens_d:.4f}")

alpha = 0.05
print(f"\n=== Conclusion ===")
if p_value < alpha:
    print(f"Since p-value ({p_value:.4f}) < significance level ({alpha}), the average purchase amounts of")
    print("Variant A and Variant B differ significantly. Switching to Variant B is worth considering.")
else:
    print(f"Since p-value ({p_value:.4f}) >= significance level ({alpha}), there is no statistically significant difference.")

Execution Result:

=== Welch's t-test (manual implementation) ===
Variant A: mean=4060.81 yen, std dev=817.67 yen, n=60
Variant B: mean=4553.63 yen, std dev=906.21 yen, n=65

t-statistic: -3.1960
Degrees of freedom (approx.): 122.94
p-value (two-sided): 0.0018

=== Cross-check with scipy.stats.ttest_ind ===
t-statistic: -3.1960
p-value: 0.0018

Effect size (Cohen's d): 0.5698

=== Conclusion ===
p-value (0.0018) < significance level (0.05), the average purchase amounts of
Variant A and Variant B differ significantly. Switching to Variant B is worth considering.

7. Summary and Next Steps

In this chapter, we examined — both in theory and through implementation — how the statistical knowledge covered throughout this series is put to use inside representative machine learning algorithms.

Revisiting the Learning Objectives

Let's revisit the learning objectives set out at the start of this chapter.

• ✅ Understand the statistical foundations of machine learning algorithms — Sections 1, 2, and 3 examined the statistical underpinnings of linear regression, logistic regression, and Naive Bayes
• ✅ Apply statistical knowledge to machine learning — confirmed in each section through implementations combining scikit-learn with numpy/scipy
• ✅ Quantify prediction uncertainty — Section 4 implemented the computation of predictive variance with a Gaussian process
• ✅ Perform model evaluation statistically — Sections 5 and 6 implemented cross validation, the paired t-test, the two-sample t-test, and effect size

Looking Back at the Series as a Whole

This series, "Introduction to Statistics for Machine Learning," concludes with this fifth chapter. We began with the fundamentals of descriptive statistics and probability in Chapter 1, broadened our view to probability distributions in Chapter 2, built the skeleton of statistical inference — estimation and hypothesis testing — in Chapter 3, learned the modern perspective of Bayesian statistics in Chapter 4, and in this final chapter confirmed how all of it connects to machine learning practice. Statistics is never a mere "opening act" for machine learning — it remains the foundation for understanding how models behave, correctly interpreting their results, and making decisions with confidence. Please carry the statistical perspective you've gained here forward into more advanced machine learning study and practice.

Practice Problems

import numpy as np
from sklearn.linear_model import LinearRegression

# Advertising spend and sales data (both in units of 10,000 yen)
ad_spend = np.array([10, 15, 20, 25, 30, 35, 40, 45, 50, 55]).reshape(-1, 1)
sales = np.array([120, 145, 158, 190, 210, 225, 260, 275, 300, 320])

model = LinearRegression()
model.fit(ad_spend, sales)

r_squared = model.score(ad_spend, sales)
print(f"Intercept: {model.intercept_:.4f}")
print(f"Slope: {model.coef_[0]:.4f}")
print(f"Coefficient of determination R^2: {r_squared:.4f}")

# Predicted sales when advertising spend is 42 (10,000 yen)
pred = model.predict(np.array([[42]]))
print(f"Predicted sales at advertising spend of 42: {pred[0]:.2f} (10,000 yen)")

# Output:
# Intercept: 74.7394
# Slope: 4.4788
# Coefficient of determination R^2: 0.9965
# Predicted sales at advertising spend of 42: 262.85 (10,000 yen)

import numpy as np
from sklearn.datasets import load_wine
from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB
from sklearn.metrics import accuracy_score

# Use the Wine dataset (classifying three cultivars from wine chemical composition)
wine = load_wine()
X, y = wine.data, wine.target

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=0, stratify=y
)

gnb = GaussianNB()
gnb.fit(X_train, y_train)
y_pred = gnb.predict(X_test)

accuracy = accuracy_score(y_test, y_pred)
print(f"Test accuracy: {accuracy:.4f}")

# Check the prior probabilities and the posterior probability of the first test sample
print("\nPrior probability for each class:")
for name, prior in zip(wine.target_names, gnb.class_prior_):
    print(f"  {name}: {prior:.4f}")

posterior = gnb.predict_proba(X_test[:1])[0]
posterior_str = ", ".join(f"{p:.4e}" for p in posterior)
print(f"\nPosterior probability of the first test sample: [{posterior_str}]")
print(f"Predicted class: {wine.target_names[y_pred[0]]}, true class: {wine.target_names[y_test[0]]}")

# Output:
# Test accuracy: 0.9630
#
# Prior probability for each class:
#   class_0: 0.3306
#   class_1: 0.4032
#   class_2: 0.2661
#
# Posterior probability of the first test sample: [1.0000e+00, 1.0339e-13, 1.9798e-38]
# Predicted class: class_0, true class: class_0

import numpy as np
from scipy import stats

# A/B test: comparing time spent on the page (in seconds) after opening two email subject line variants (A/B)
np.random.seed(1)
group_a = np.random.normal(loc=45, scale=12, size=40)  # Variant A (existing subject line)
group_b = np.random.normal(loc=52, scale=13, size=42)  # Variant B (new subject line)

t_stat, p_value = stats.ttest_ind(group_a, group_b, equal_var=False)

pooled_std = np.sqrt(((len(group_a) - 1) * group_a.var(ddof=1)
                       + (len(group_b) - 1) * group_b.var(ddof=1))
                      / (len(group_a) + len(group_b) - 2))
cohens_d = (group_b.mean() - group_a.mean()) / pooled_std

print(f"Variant A: mean time on page={group_a.mean():.2f} sec, n={len(group_a)}")
print(f"Variant B: mean time on page={group_b.mean():.2f} sec, n={len(group_b)}")
print(f"\nt-statistic: {t_stat:.4f}")
print(f"p-value: {p_value:.4f}")
print(f"Effect size (Cohen's d): {cohens_d:.4f}")

alpha = 0.05
if p_value < alpha:
    print(f"\nSince p-value < {alpha}, the difference is statistically significant")
else:
    print(f"\nSince p-value >= {alpha}, the difference is not statistically significant")

# Output:
# Variant A: mean time on page=43.92 sec, n=40
# Variant B: mean time on page=54.20 sec, n=42
#
# t-statistic: -3.9320
# p-value: 0.0002
# Effect size (Cohen's d): 0.8683
#
# Since p-value < 0.05, the difference is statistically significant