09. 単回帰分析をゼロから実装
テーマ: 最小二乗法（OLS）による広告費と売上の関係の推定
学習点: 数式の実装, t検定, 決定係数, 残差分析, scipy との照合
依存: SciPy ／ 難易度: 中級
実行方法
uv run 09_ols_simple.py
スクリプト冒頭の PEP 723 メタデータ（# /// script）により、必要なライブラリは
uv が自動的に仮想環境へ導入します。事前の pip install は不要です。
解説
何をするプログラムか
「広告費を増やせば売上はどれだけ伸びるのか」というマーケティングの典型的な問いを、単回帰分析（最小二乗法, OLS）で推定するスクリプトです。15 期分の広告費と売上（いずれも百万円）のデータから、売上 = α + β×広告費 + u というモデルの係数を推定し、係数の標準誤差・t 値・p 値・95% 信頼区間・決定係数 R² までを教科書の公式どおりに自前実装します。
推定だけで終わらず、scipy.stats.linregress と 10 桁一致することの照合、残差プロットによる線形仮定の目視チェック、Durbin-Watson 統計量による系列相関の確認まで行う、実証分析のミニチュア版です。
コードの読みどころ
ols()関数の中心は、偏差平方和sxxと偏差積和sxyをジェネレータ式で計算し、beta = sxy / sxx、alpha = my - beta * mxと推定量の公式をそのまま書き下している部分です。数式とコードが 1 対 1 に対応します。- 誤差分散は
s2 = rss / (n - 2)と自由度 n − 2 で割ります。α と β の 2 つを推定した分だけ自由度が減る、という統計学の約束がコード上の- 2に現れています。 - 検定は
p = 2 * (1 - stats.t.cdf(abs(r["t_b"]), df))で両側 p 値を、tcrit = stats.t.ppf(0.975, df)で信頼区間用の臨界値を求めます。SciPy は t 分布の確率計算にだけ使い、推定そのものは標準の算術で完結しています。 assert abs(r["beta"] - lr.slope) < 1e-9は、自前実装がscipy.stats.linregressと一致しなければ実行を止める検証です。「自分の実装を信頼できる実装と突き合わせる」作法の実例です。- 残差プロットは、各残差を最大絶対値
scaleで正規化して 41 文字幅の行の位置posに*を置くだけのテキスト描画です。グラフライブラリなしでもパターンの有無は十分見て取れます。 - Durbin-Watson 統計量は隣り合う残差の差の平方和を RSS で割ったもので、ジェネレータ式 1 つで計算しています。
理論的背景
OLS は残差平方和を最小にする直線を選ぶ方法で、傾きの推定量は β̂ = Σ(x−x̄)(y−ȳ) / Σ(x−x̄)² です。その標準誤差は SE(β̂) = √(s²/Σ(x−x̄)²)、s² = RSS/(n−2) で与えられ、誤差項の正規性・等分散・独立の仮定のもとで t = β̂/SE(β̂) は帰無仮説 β = 0 の下で自由度 n−2 の t 分布に従います。R² は売上の変動のうちモデルで説明できる割合、自由度調整済み R² は説明変数を増やすことへのペナルティを加えた版です。
実行結果の見方
β̂ = 4.9366（SE 0.0891）、t = 55.4 で p 値は表示桁上ゼロとなり、「広告費 100 万円の増加は売上約 494 万円の増加と関連する」ことが強く支持されます。R² = 0.9958 とあてはまりも極めて良好です。scipy との照合では slope・intercept・stderr が 10 桁まで一致しており、自前実装の正しさが確認できます。ただし残差プロットを見ると、x の小さい端と大きい端で負、中央で正という緩やかな山型のパターンがあり、Durbin-Watson 統計量も 0.549 と 2 を大きく下回ります。係数が有意でも残差診断までは合格とは限らない、という残差分析の重要性を示す結果です。
ソースコード
# /// script
# requires-python = ">=3.11"
# dependencies = [
# "scipy",
# ]
# ///
"""09: 単回帰分析をゼロから実装 -------------------------------------------
テーマ: 最小二乗法（OLS）による広告費と売上の関係の推定
学習点: 数式の実装, t検定, 決定係数, 残差分析, scipy との照合
根拠: β̂ = Σ(x-x̄)(y-ȳ) / Σ(x-x̄)^2、SE(β̂) = √(s^2 / Σ(x-x̄)^2)、
s^2 = RSS/(n-2)。t = β̂/SE(β̂) は帰無仮説 β=0 の下で自由度 n-2 の
t 分布に従う（誤差項の正規性・等分散・独立を仮定）。
"""
import math
from scipy import stats
def ols(x, y):
n = len(x)
mx, my = sum(x) / n, sum(y) / n
sxx = sum((xi - mx) ** 2 for xi in x)
sxy = sum((xi - mx) * (yi - my) for xi, yi in zip(x, y))
beta = sxy / sxx
alpha = my - beta * mx
fitted = [alpha + beta * xi for xi in x]
resid = [yi - fi for yi, fi in zip(y, fitted)]
rss = sum(r ** 2 for r in resid)
tss = sum((yi - my) ** 2 for yi in y)
s2 = rss / (n - 2) # 誤差分散の不偏推定量
se_b = math.sqrt(s2 / sxx)
se_a = math.sqrt(s2 * (1 / n + mx ** 2 / sxx))
r2 = 1 - rss / tss
adj_r2 = 1 - (1 - r2) * (n - 1) / (n - 2)
return {
"alpha": alpha, "beta": beta, "se_a": se_a, "se_b": se_b,
"t_b": beta / se_b, "r2": r2, "adj_r2": adj_r2,
"rss": rss, "n": n, "fitted": fitted, "resid": resid,
"sigma": math.sqrt(s2),
}
def main() -> None:
# 広告費(百万円) と 売上(百万円)
ad = [1.2, 2.0, 2.4, 3.1, 3.5, 4.0, 4.6, 5.2, 5.8, 6.5,
7.0, 7.7, 8.1, 8.9, 9.4]
sales = [22.1, 27.3, 29.8, 33.4, 36.0, 37.9, 41.2, 44.8, 47.1, 50.9,
53.2, 56.0, 57.3, 61.4, 63.0]
r = ols(ad, sales)
df = r["n"] - 2
p = 2 * (1 - stats.t.cdf(abs(r["t_b"]), df))
tcrit = stats.t.ppf(0.975, df)
print("=== 単回帰 売上 = α + β×広告費 + u ===")
print(f" n = {r['n']}, 自由度 = {df}")
print(f" α̂ = {r['alpha']:>8.4f} (SE {r['se_a']:.4f})")
print(f" β̂ = {r['beta']:>8.4f} (SE {r['se_b']:.4f}) "
f"t = {r['t_b']:.3f} p = {p:.3e}")
print(f" 95%信頼区間 β: "
f"[{r['beta'] - tcrit * r['se_b']:.4f}, "
f"{r['beta'] + tcrit * r['se_b']:.4f}]")
print(f" R² = {r['r2']:.4f} / 自由度調整済 R² = {r['adj_r2']:.4f}")
print(f" 残差標準誤差 σ̂ = {r['sigma']:.4f}")
print(f"\n 解釈: 広告費を100万円増やすと売上は約 {r['beta']:.2f} 百万円増加"
f"（p<0.001 で有意）。")
print("\n[scipy.stats.linregress との照合]")
lr = stats.linregress(ad, sales)
print(f" slope : 自前 {r['beta']:.10f} / scipy {lr.slope:.10f}")
print(f" intercept : 自前 {r['alpha']:.10f} / scipy {lr.intercept:.10f}")
print(f" stderr : 自前 {r['se_b']:.10f} / scipy {lr.stderr:.10f}")
assert abs(r["beta"] - lr.slope) < 1e-9
print("\n[残差プロット（テキスト）] 系統的パターンがなければ線形仮定は妥当")
scale = max(abs(v) for v in r["resid"])
for xi, res in zip(ad, r["resid"]):
pos = int(20 + res / scale * 18)
line = [" "] * 41
line[20] = "|"
line[pos] = "*"
print(f" x={xi:>4.1f} {''.join(line)} {res:+.3f}")
dw = (sum((r["resid"][i] - r["resid"][i - 1]) ** 2
for i in range(1, len(r["resid"]))) / r["rss"])
print(f"\n Durbin-Watson統計量 = {dw:.3f} (2付近なら系列相関なし)")
if __name__ == "__main__":
main()
実行結果
=== 単回帰 売上 = α + β×広告費 + u ===
n = 15, 自由度 = 13
α̂ = 17.9621 (SE 0.5218)
β̂ = 4.9366 (SE 0.0891) t = 55.419 p = 0.000e+00
95%信頼区間 β: [4.7442, 5.1291]
R² = 0.9958 / 自由度調整済 R² = 0.9955
残差標準誤差 σ̂ = 0.8652
解釈: 広告費を100万円増やすと売上は約 4.94 百万円増加（p<0.001 で有意）。
[scipy.stats.linregress との照合]
slope : 自前 4.9366333555 / scipy 4.9366333555
intercept : 自前 17.9620874383 / scipy 17.9620874383
stderr : 自前 0.0890787424 / scipy 0.0890787424
[残差プロット（テキスト）] 系統的パターンがなければ線形仮定は妥当
x= 1.2 * | -1.786
x= 2.0 * | -0.535
x= 2.4 *| -0.010
x= 3.1 |* +0.134
x= 3.5 | * +0.760
x= 4.0 |* +0.191
x= 4.6 | * +0.529
x= 5.2 | * +1.167
x= 5.8 | * +0.505
x= 6.5 | * +0.850
x= 7.0 | * +0.681
x= 7.7 * +0.026
x= 8.1 * | -0.649
x= 8.9 * | -0.498
x= 9.4 * | -1.366
Durbin-Watson統計量 = 0.549 (2付近なら系列相関なし)