A from-scratch simple linear regression (one x, one y) Python class
The whole point of building this by hand instead of just calling scikit-learn is to SEE the machinery: how the line is fitted, where the error numbers come from, and why we can say "we're 95% confident the slope is between A and B" instead of just reporting a single number.
import numpy as np
from scipy import stats
import matplotlib
matplotlib.use("TkAgg") # must come before importing pyplot
import matplotlib.pyplot as plt
class linear_regression:
"""
A from-scratch simple linear regression (one x, one y).
Usage:
lr = linear_regression(x, y)
lr.fit() # runs every calculation, in the right order
lr.summarize() # prints a human-readable report
lr.predict(60) # predict y for a new x, with a confidence range
lr.plot() # draws the data + fitted line + interval band
"""
def __init__(self, x, y, alpha=0.05):
# alpha is the "risk" we're willing to accept when we say we're
# confident about something. alpha = 0.05 means "95% confident"
# (1 - alpha). This is the standard default used almost everywhere.
if len(x) != len(y):
raise Exception("x and y length don't match")
self.count = len(x)
self.number_of_predictors = 1 # Only a single x
# Degrees of freedom = "how many independent pieces of information
# are left over after we've used some of the data to estimate the
# model itself". We used up 2 pieces of info to estimate the line
# (1 for the slope, 1 for the intercept), so we subtract those.
# This number matters because it controls how "wide" or "narrow"
# our t-distribution is - fewer data points -> wider, more
# cautious intervals.
self.degrees_of_freedom = len(x) - 1 - self.number_of_predictors
self.alpha = alpha
# critical_t is the "how many standard errors away from the middle
# do we need to go to capture 95% of the possibilities" number.
# It comes from the t-distribution rather than the normal
# distribution because we're estimating things from a limited
# sample, not from the whole population - the t-distribution
# automatically gets wider (more forgiving) for small samples.
self.critical_t = stats.t.ppf(1 - self.alpha / 2, self.degrees_of_freedom)
self.x = x
self.y = y
# Everything below is set to None on purpose. None of these values
# exist until fit() actually runs the calculations. Declaring them
# here isn't strictly required by Python, but it means anyone
# reading this class (future me included) can see, in one place,
# everything this model is capable of producing - before any math
# has actually happened.
self.mean_x = None
self.mean_y = None
self.var_x = None
self.var_y = None
self.difference_x = []
self.difference_y = []
self.cov = None
self.sd_x = None
self.sd_y = None
self.R = None
self.R_squared = None
self.slope = None
self.intercept = None
self.residuals = None
self.abs_residuals = None
self.sqr_residuals = None
self.mae = None
self.mse = None
self.rmse = None
self.adjusted_R_squared = None
self.RSS = None
self.y_sum_of_squares = None
self.x_sum_of_squares = None
self.slope_error = None
self.intercept_error = None
self.t_statistic_for_slope = None
self.p_value_for_slope = None
self.slope_bounds = [None, None]
self.intercept_bounds = [None, None]
self.predicted_bounds = [None, None]
def calc_mean(self):
# The average of x and the average of y. Nothing fancy - this is
# the "center point" that the whole regression line is built
# around. The fitted line is guaranteed to pass through
# (mean_x, mean_y).
self.mean_x = np.mean(self.x)
self.mean_y = np.mean(self.y)
return self.mean_x, self.mean_y
def calc_variance(self):
# Variance = "on average, how far do the values spread out from
# their own mean, squared". We square the distances so that
# points above and below the mean don't cancel each other out.
self.var_x = np.var(self.x)
self.var_y = np.var(self.y)
return self.var_x, self.var_y
def calc_difference(self):
# For every point, how far is it from the mean? We need these
# "distance from center" numbers as building blocks for
# covariance, standard deviation, and correlation below.
self.difference_x = [x - self.mean_x for x in self.x]
self.difference_y = [y - self.mean_y for y in self.y]
return self.difference_x, self.difference_y
def calc_cov(self):
# Covariance - quantifies the joint variability between the independent variable(s) and the dependent variable.
# In plain terms: when x goes up, does y also tend to go up (positive
# covariance), go down (negative), or is there no pattern (near zero)?
self.cov = sum([x * y for x, y in zip(self.difference_x, self.difference_y)]) / self.count
return self.cov
def calc_sd(self):
# Standard deviation is just the square root of variance - it
# brings the number back into the same units as the original data
# (variance is in "squared" units, which isn't intuitive to read).
variance_x = sum(d ** 2 for d in self.difference_x) / self.count
variance_y = sum(d ** 2 for d in self.difference_y) / self.count
self.sd_x = variance_x ** 0.5
self.sd_y = variance_y ** 0.5
return self.sd_x, self.sd_y
def calc_R(self):
# R (the correlation coefficient) rescales covariance into a
# fixed range between -1 and 1, which makes it comparable across
# completely different datasets. +1 = perfect positive
# relationship, -1 = perfect negative, 0 = no linear relationship.
self.R = self.cov / (self.sd_x * self.sd_y)
return self.R
def calc_Rsquared(self):
# To test whether the model is able to describe the data, we use R².
# If the model accurately describes the data, the value of R² will be 1.
# If R² is 0.75, it means that 75% of the variation in Y is strictly driven by X values. The remaining 25% comes from other factors.
# The problem with R² is that we can get high values of it even though our model is not successful.
# For example, through 2 points only one straight line equation can be passed, so R² for a model based on 2 samples
# will be 1, but in reality the sample is too small for us to learn about the world from it. Another problem is that the more independent variables there are,
# the more R² may increase because in the best case the parameters that do not describe the model at all will not add to R², and in the worst case they will make it higher.
self.R_squared = self.R ** 2
return self.R_squared
def calc_slope(self):
# The slope answers "for every +1 change in x, how much does y
# change on average?". Dividing covariance (how x and y move
# together) by the variance of x (how spread out x is) gives the
# best-fitting line's steepness. This is the "least squares"
# slope - the one that minimizes the total squared error.
self.slope = self.cov / self.var_x
return self.slope
def calc_intercept(self):
# Once we know the slope, the intercept is whatever value makes
# the line pass through the mean point (mean_x, mean_y) - see the
# note in calc_mean() above.
self.intercept = self.mean_y - self.slope * self.mean_x
return self.intercept
def calc_residuals(self):
# Residuals - errors of the model
# For every real data point, how far off was our line's guess?
# residual = actual_y - predicted_y. We also keep an absolute
# version (ignores direction, just "how big was the miss") and a
# squared version (punishes big misses much more than small ones)
# because different error metrics below need different versions.
self.residuals = [self.y[idx] - (self.slope * self.x[idx] + self.intercept) for idx in range(self.count)]
self.abs_residuals = [np.abs(res) for res in self.residuals]
self.sqr_residuals = [res ** 2 for res in self.residuals]
return self.residuals, self.abs_residuals, self.sqr_residuals
def calc_mae(self):
# MAE (Mean Absolute Error): Measures the average size of the mistakes.
# Easy to read: "on average, our predictions are off by MAE units".
self.mae = np.mean(self.abs_residuals)
return self.mae
def calc_mse(self):
# MSE (Mean Squared Error): Punishes larger mistakes heavily by squaring them.
# Harder to read directly (it's in squared units) but it's the
# quantity the least-squares line-fitting is actually minimizing.
self.mse = np.mean(self.sqr_residuals)
return self.mse
def calc_rmse(self):
# RMSE (Root Mean Squared Error): Brings the MSE back into the original currency/units.
# A good "typical error size" number to put in front of someone
# who isn't going to think in squared units.
self.rmse = np.sqrt(self.mse)
return self.rmse
def calc_adjusted_R_squared(self):
# R² can only stay flat or go up as you add more predictors, even useless ones, because a linear regression's least-squares fit can always exploit an extra dimension to reduce leftover error,
# even if that predictor has zero real relationship to y. So plain R² is a bad way to compare models with different numbers of predictors since it always rewards throwing more variables in,
# regardless of whether they help.
# Adjusted R² corrects for this by penalizing for each additional predictor, so it only increases when a new variable improves the fit by more than chance would predict.
# Note to future me: number_of_predictors is always 1 in this class
# right now, so adjusted_R_squared will stay very close to R_squared.
# The penalty only becomes meaningful once this class is extended
# to handle multiple x variables (multiple regression).
self.adjusted_R_squared = 1 - ((1 - self.R_squared) * (self.count - 1) / (self.count - self.number_of_predictors - 1))
return self.adjusted_R_squared
def calc_rss(self):
# RSS (Residual Sum of Squares)
# Just the sum (not the average) of all the squared errors. This
# is a stepping stone toward the "residual standard error" below -
# think of it as "total leftover error the model couldn't explain".
self.RSS = np.sum(self.sqr_residuals)
return self.RSS
def calc_y_sum_of_squares(self):
# Residual Standard Error
# This is like RMSE's more "honest" cousin: instead of dividing
# total error by count (n), we divide by degrees_of_freedom
# (n - 2), which corrects for the fact that we already used 2 of
# our data points' worth of information to estimate the slope and
# intercept. This is the error measure we use for all the
# confidence/prediction interval math below - not RMSE.
self.y_sum_of_squares = np.sqrt(self.RSS / self.degrees_of_freedom)
return self.y_sum_of_squares
def calc_x_sum_of_squares(self):
# SS_xx (Sum of Squares for x)
# Measures how spread out our x values are. We'll need this in
# the slope/intercept/prediction error formulas below - the more
# spread out x is, the more confidently we can pin down the slope.
self.x_sum_of_squares = np.sum([(x-self.mean_x) ** 2 for x in self.x])
return self.x_sum_of_squares
def calc_slope_error(self):
# Standard Error of Slope (SE_a)
# How much would our estimated slope likely wobble if we re-ran
# this experiment with a new random sample? Smaller = more
# confident about the slope we calculated.
self.slope_error = self.y_sum_of_squares / np.sqrt(self.x_sum_of_squares)
return self.slope_error
def calc_intercept_error(self):
# Standard Error of the intercept
# Same idea as slope error, but for the intercept.
self.intercept_error = self.y_sum_of_squares * np.sqrt(1 / self.count + self.mean_x ** 2 / self.x_sum_of_squares)
return self.intercept_error
def calc_t_statistic_for_slope(self):
# A "signal to noise" ratio for the slope: how many standard
# errors away from zero is our slope? If this number is large
# (in either direction), it's a sign the slope probably isn't
# just random noise - x probably really does affect y.
self.t_statistic_for_slope = self.slope / self.slope_error
return self.t_statistic_for_slope
def calc_p_value_for_slope(self):
# Survival function (sf) is 1 - CDF. We multiply by 2 for a two-tailed test.
# The p-value translates the t-statistic above into "if x had NO
# real effect on y at all, what's the probability we'd see a
# slope this extreme just by random chance?". Small p-value
# (conventionally < 0.05) = the slope is probably real, not luck.
self.p_value_for_slope = 2 * stats.t.sf(np.abs(self.t_statistic_for_slope), self.degrees_of_freedom)
return self.p_value_for_slope
def calc_slope_bounds(self):
# Turns "here's our single best-guess slope" into "here's a range
# we're 95% confident the TRUE slope falls within". Same shape as
# every other confidence interval: estimate +/- (critical value *
# standard error).
self.slope_bounds = [self.slope - self.critical_t * self.slope_error, self.slope + self.critical_t * self.slope_error]
return self.slope_bounds
def calc_intercept_bounds(self):
# Same idea as calc_slope_bounds, but for the intercept.
self.intercept_bounds = [self.intercept - self.critical_t * self.intercept_error, self.intercept + self.critical_t * self.intercept_error]
return self.intercept_bounds
def fit(self):
# This is the "run everything" method. The ORDER here isn't
# arbitrary - each calc_* method reads values that an earlier one
# produced (e.g. calc_slope needs calc_cov and calc_variance to
# have already run). Calling fit() means future me never has to
# remember that dependency chain by hand.
self.calc_mean()
self.calc_variance()
self.calc_difference()
self.calc_cov()
self.calc_sd()
self.calc_R()
self.calc_Rsquared()
self.calc_slope()
self.calc_intercept()
self.calc_residuals()
self.calc_mae()
self.calc_mse()
self.calc_rmse()
self.calc_adjusted_R_squared()
self.calc_rss()
self.calc_y_sum_of_squares()
self.calc_x_sum_of_squares()
self.calc_slope_error()
self.calc_intercept_error()
self.calc_t_statistic_for_slope()
self.calc_p_value_for_slope()
self.calc_slope_bounds()
self.calc_intercept_bounds()
return self # returning self lets us chain: linear_regression(x, y).fit()
def summarize(self):
# A human-readable report of everything the model calculated.
# This is the "don't make me go dig through 20 attributes" method.
lines = []
lines.append("=" * 50)
lines.append("Linear Regression Summary")
lines.append("=" * 50)
lines.append(f"{'Observations:':<25}{self.count}")
lines.append(f"{'Degrees of freedom:':<25}{self.degrees_of_freedom}")
lines.append(f"{'Alpha:':<25}{self.alpha}")
lines.append("-" * 50)
lines.append("-" * 50)
lines.append('The core metrics\nThese metrics measure the "errors" (residuals) of your model.')
lines.append(f"{' MAE:':<25}{self.mae:.4f}")
lines.append(f"{' MSE:':<25}{self.mse:.4f}")
lines.append(f"{' RMSE:':<25}{self.rmse:.4f}")
lines.append("-" * 50)
lines.append(f"{'R:':<25}{self.R:.4f}")
lines.append(f"{' R-squared:':<25}{self.R_squared:.4f}")
lines.append(f"{' Adjusted R-squared:':<25}{self.adjusted_R_squared:.4f}")
lines.append("Standard R² can falsely increase if you add useless variables. Adjusted R² penalizes the model for the number of predictors used (sample size).")
lines.append(f"{'Intercept:':<25}{self.intercept:.4f}")
lines.append("The y-intercept (b) represents the predicted y when x = 0.")
lines.append(f"{' 95% CI:':<25}[{self.intercept_bounds[0]:.4f}, {self.intercept_bounds[1]:.4f}]")
lines.append(f"{' Std. Error:':<25}{self.intercept_error:.4f}")
lines.append("-" * 50)
lines.append(f"{'Slope:':<25}{self.slope:.4f}")
lines.append(f"{' 95% CI:':<25}[{self.slope_bounds[0]:.4f}, {self.slope_bounds[1]:.4f}]")
lines.append('Instead of a single "guess," a confidence interval gives you a range (e.g., 95% confident) where the true market slope lies.')
lines.append(f"{' Std. Error:':<25}{self.slope_error:.4f}")
lines.append(f"{' t-statistic:':<25}{self.t_statistic_for_slope:.4f}")
lines.append(f"{' p-value:':<25}{self.p_value_for_slope:.4g}")
lines.append("Hypothesis test against the assumption (H₀) that the slope is 0")
lines.append("=" * 50)
report = "\n".join(lines)
print(report)
return report
def to_dict(self):
# Same information as summarize(), but as plain nested
# dictionaries/numbers instead of formatted text - meant for
# code to consume (e.g. json.dumps this straight into a PHP/JS
# frontend), not for a person to read on a terminal.
return {
"n": self.count,
"degrees_of_freedom": self.degrees_of_freedom,
"alpha": self.alpha,
"intercept": {
"value": float(self.intercept),
"std_error": float(self.intercept_error),
"ci": [float(b) for b in self.intercept_bounds],
},
"slope": {
"value": float(self.slope),
"std_error": float(self.slope_error),
"ci": [float(b) for b in self.slope_bounds],
"t_statistic": float(self.t_statistic_for_slope),
"p_value": float(self.p_value_for_slope),
},
"fit": {
"r": float(self.R),
"r_squared": float(self.R_squared),
"adjusted_r_squared": float(self.adjusted_R_squared),
},
"error": {
"mae": float(self.mae),
"mse": float(self.mse),
"rmse": float(self.rmse),
},
}
def plot(self, show=True):
# Draws: the raw data points, the fitted line, and a shaded band
# showing the prediction interval. The band gets WIDER the
# farther x is from mean_x - see the note in predict() below for
# why. That widening is the whole visual point of this chart:
# it shows we're less sure about predictions far from our data.
x_line = np.linspace(min(self.x), max(self.x), 200)
fitted = [self.intercept + self.slope * xi for xi in x_line]
bounds = [self.predict(xi) for xi in x_line]
lower = [b[0] for b in bounds]
upper = [b[2] for b in bounds]
fig, ax = plt.subplots(figsize=(8, 6))
ax.scatter(self.x, self.y, color="black", s=25, label="Data", zorder=3)
ax.plot(x_line, fitted, color="crimson", linewidth=2, label="Fitted line")
ax.fill_between(x_line, lower, upper, color="crimson", alpha=0.15,
label=f"{int((1 - self.alpha) * 100)}% prediction interval")
ax.set_title(f"Linear Regression (R² = {self.R_squared:.3f}, p = {self.p_value_for_slope:.4g})")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.legend()
fig.tight_layout()
if show:
plt.show()
return fig, ax
def predict(self, x):
# Quantiles (Prediction Intervals vs. Confidence Intervals)
# When evaluating a specific x value, statisticians look at two kinds of intervals:
# 1. Confidence Interval (For the Average): Where the average y of all xs that the model was trained on lies.
# 2. Prediction Interval: Predicted value y for a specific given x.
#
# This method returns a Prediction Interval, which is intentionally
# WIDER than a confidence interval, because predicting one new,
# individual point is inherently less certain than estimating an
# average.
#
# Note the "(x - self.mean_x) ** 2" term inside the square root:
# this is why the interval gets wider the further x is from
# mean_x. Makes sense intuitively - we trust our line most near
# the center of the data we actually observed, and trust it less
# the further we extrapolate away from it.
predicted = self.intercept + self.slope * x
standard_error = self.y_sum_of_squares * np.sqrt(
1 + 1 / self.count + (x - self.mean_x) ** 2 / self.x_sum_of_squares
)
lower_predicted = predicted - self.critical_t * standard_error
upper_predicted = predicted + self.critical_t * standard_error
return lower_predicted, predicted, upper_predicted
# Test cases
x = [60, 70, 76, 85, 100, 110]
y = [124, 133, 96, 155, 185, 105]
lr = linear_regression(x, y)
summary = lr.fit().summarize()
print(summary)
lr.plot()
print(lr.predict(60))