← 目次 ｜ ← 前: 07 ｜ 次: 09 →

08. 記述統計をゼロから実装

テーマ: 平均・分散・分位点・歪度・尖度を自前実装し標準ライブラリと照合

学習点: sorted, sum, statistics モジュール, 不偏推定量, アサーション

依存: 標準ライブラリのみ ／ 難易度: 基礎

実行方法

uv run 08_descriptive_stats.py

スクリプト冒頭の PEP 723 メタデータ（# /// script）により、必要なライブラリは uv が自動的に仮想環境へ導入します。事前の pip install は不要です。

解説

何をするプログラムか

データ分析の第一歩である記述統計を、ライブラリに頼らず自分の手で実装するスクリプトです。題材は 20 人分の年収データ（万円）で、大半が 300〜700 万円台に集まる一方、980 万円と 2,400 万円という高所得者を含む「右に歪んだ分布」になっています。所得・売上・資産といった経済データに典型的な形です。

平均・中央値・不偏分散・標準偏差・変動係数・歪度・超過尖度・分位点を計算し、IQR（四分位範囲）による外れ値検出を行ったうえで、自前実装が標準ライブラリの statistics モジュールと一致することをアサーションで検証します。最後に、最大値 1 つを除くだけで平均が 15.6% も動くのに中央値はほとんど動かないことを示し、「所得の代表値には中央値」という統計実務の定石を裏づけます。

コードの読みどころ

理論的背景

不偏分散が偏差平方和を n ではなく n − 1 で割る（ベッセル補正）のは、偏差を「真の平均」ではなく「標本平均」から測ることで平方和が系統的に小さくなるのを補正し、E[Σ(x − x̄)²] = (n − 1)σ² という関係から母分散の不偏推定量を得るためです。歪度は分布の非対称性（正なら右裾が長い）、超過尖度は正規分布を 0 とした裾の厚さの指標で、標準化した偏差の 3 乗・4 乗の平均として定義されます。

実行結果の見方

平均 605.75 万円に対して中央値は 475.00 万円と大きく下回っており、高所得側の外れ値が平均を引き上げていることが分かります。歪度 3.421、超過尖度 11.167 という大きな正の値も、右裾が長く裾の厚い分布であることを裏づけます。IQR 基準では 980 と 2400 が外れ値と判定されました。最後のブロックでは、最大値 2,400 万円を 1 つ除くだけで平均が 605.8 → 511.3 万円（−15.6%）へ動く一方、中央値は 475.0 → 470.0 万円とほぼ不変です。政府統計の所得分布で中央値が併記される理由を、自分で計算した数値として確認できます。

ソースコード

# /// script
# requires-python = ">=3.11"
# dependencies = []
# ///
"""08: 記述統計をゼロから実装 ---------------------------------------------
テーマ: 平均・分散・分位点・歪度・尖度を自前実装し標準ライブラリと照合
学習点: sorted, sum, statistics モジュール, 不偏推定量, アサーション
根拠: 標本分散を n-1 で割る（ベッセル補正）のは、母分散の不偏推定量に
      するため。E[Σ(x-x̄)^2] = (n-1)σ^2 が成り立つ。
"""
import math
import statistics as st


def mean(xs): return sum(xs) / len(xs)


def variance(xs, ddof: int = 1):
    """ddof=1 で不偏分散、ddof=0 で標本分散。"""
    m = mean(xs)
    return sum((x - m) ** 2 for x in xs) / (len(xs) - ddof)


def quantile(xs, q: float):
    """線形補間による分位点（numpy の既定 'linear' と同じ定義）。"""
    s = sorted(xs)
    if len(s) == 1:
        return s[0]
    pos = (len(s) - 1) * q
    lo = math.floor(pos)
    hi = math.ceil(pos)
    return s[lo] + (s[hi] - s[lo]) * (pos - lo)


def skewness(xs):
    """歪度（母集団定義）。正なら右に裾が長い。"""
    m, s = mean(xs), math.sqrt(variance(xs, ddof=0))
    return sum(((x - m) / s) ** 3 for x in xs) / len(xs)


def kurtosis_excess(xs):
    """超過尖度。正規分布で 0 になるよう 3 を引く。"""
    m, s = mean(xs), math.sqrt(variance(xs, ddof=0))
    return sum(((x - m) / s) ** 4 for x in xs) / len(xs) - 3


def main() -> None:
    # 年収データ（万円）— 少数の高所得者を含む右に歪んだ分布
    income = [320, 350, 380, 400, 410, 425, 430, 450, 460, 470,
              480, 500, 520, 540, 560, 610, 680, 750, 980, 2400]

    print(f"n = {len(income)}")
    print(f"平均    : {mean(income):>9.2f} 万円")
    print(f"中央値  : {quantile(income, 0.5):>9.2f} 万円  <- 外れ値に頑健")
    print(f"最頻階級: 400万円台")
    print(f"不偏分散: {variance(income):>9.2f}")
    print(f"標準偏差: {math.sqrt(variance(income)):>9.2f} 万円")
    print(f"変動係数: {math.sqrt(variance(income)) / mean(income):>9.3f}")
    print(f"歪度    : {skewness(income):>9.3f}  (>0: 右裾が長い)")
    print(f"超過尖度: {kurtosis_excess(income):>9.3f}  (>0: 裾が厚い)")

    print("\n[分位点]")
    for q in [0.10, 0.25, 0.50, 0.75, 0.90]:
        print(f"  {q:>5.0%} 分位: {quantile(income, q):>8.2f} 万円")
    iqr = quantile(income, 0.75) - quantile(income, 0.25)
    fence = quantile(income, 0.75) + 1.5 * iqr
    print(f"  IQR = {iqr:.2f} / 外れ値の上側境界 (Q3+1.5*IQR) = {fence:.2f}")
    print(f"  外れ値: {[x for x in income if x > fence]}")

    print("\n[標準ライブラリとの照合]")
    assert abs(mean(income) - st.mean(income)) < 1e-9
    assert abs(variance(income) - st.variance(income)) < 1e-9
    assert abs(variance(income, 0) - st.pvariance(income)) < 1e-9
    assert abs(quantile(income, 0.5) - st.median(income)) < 1e-9
    print("  すべて statistics モジュールと一致 (assert 通過)")

    print("\n[平均と中央値の乖離] 最大値を除くと")
    trimmed = income[:-1]
    print(f"  平均 {mean(income):.1f} -> {mean(trimmed):.1f} 万円 "
          f"({mean(trimmed) / mean(income) - 1:+.1%})")
    print(f"  中央値 {quantile(income, .5):.1f} -> {quantile(trimmed, .5):.1f} 万円")
    print("  ※ 所得分布の代表値に中央値が使われる理由がここにある。")


if __name__ == "__main__":
    main()

実行結果

n = 20
平均    :    605.75 万円
中央値  :    475.00 万円  <- 外れ値に頑健
最頻階級: 400万円台
不偏分散: 201719.14
標準偏差:    449.13 万円
変動係数:     0.741
歪度    :     3.421  (>0: 右裾が長い)
超過尖度:    11.167  (>0: 裾が厚い)

[分位点]
    10% 分位:   377.00 万円
    25% 分位:   421.25 万円
    50% 分位:   475.00 万円
    75% 分位:   572.50 万円
    90% 分位:   773.00 万円
  IQR = 151.25 / 外れ値の上側境界 (Q3+1.5*IQR) = 799.38
  外れ値: [980, 2400]

[標準ライブラリとの照合]
  すべて statistics モジュールと一致 (assert 通過)

[平均と中央値の乖離] 最大値を除くと
  平均 605.8 -> 511.3 万円 (-15.6%)
  中央値 475.0 -> 470.0 万円
  ※ 所得分布の代表値に中央値が使われる理由がここにある。

← 目次 ｜ ← 前: 07 ｜ 次: 09 →