最小二乗法は何を最小にしているのか？直線フィットをダミーデータとPythonで図解する

散布図に直線を当てはめる「直線フィット」は、統計解析やデータサイエンスで頻繁に登場します。

Pythonでは、たとえばnumpy.polyfit()を使うと、わずか1行で直線を求められます。

しかし、

• その直線は、どのような基準で選ばれているのでしょうか
• 「最小二乗法」は、何を最小にしているのでしょうか
• 観測値と直線の差は、どのように計算されるのでしょうか

と聞かれると、少し説明が難しいかもしれません。

この記事では、ダミーデータを作成し、直線、予測値、残差、残差平方和をPythonで可視化します。

Google Colabでコードを順番に実行しながら、最小二乗法が行っていることを図で理解することを目指します。

この記事で分かること

この記事では、次の内容を扱います。

• 説明変数、目的変数、予測値、残差という用語
• 最小二乗法が最小にする量
• 残差を二乗する理由
• 傾きと切片を変えたときの残差平方和の変化
• numpy.polyfit()による直線フィット
• 残差プロットの基本的な読み方
• 最小二乗法を使うときの注意点

先に結論

$n$個のデータ点$(x_i, y_i)$に対して、次の直線を考えます。

$$
\hat{y}_i = ax_i + b
$$

ここで、

• $a$：直線の傾き
• $b$：直線の切片
• $y_i$：実際に観測された値
• $\hat{y}_i$：直線から予測された値

です。

観測値と予測値の差を残差と呼びます。

$$
e_i = y_i - \hat{y}_i
$$

最小二乗法は、残差を二乗して足し合わせた量、

RSS(a,b)=\sum_{i=1}^{n}\left\{y_i-(ax_i+b)\right\}^2

を最小にする傾き$a$と切片$b$を選びます。

したがって、最小二乗法が最小にしているものは、

観測値と予測値の差そのものではなく、残差の二乗和

です。

この残差の二乗和を、残差平方和、英語ではResidual Sum of Squares、略してRSSと呼びます。

Google Colabでの実行方法

Google Colabで新しいノートブックを作成し、以下のコードセルを上から順番に実行してください。

1. ライブラリを読み込む

今回使用するライブラリは、次の3つです。

• NumPy：数値計算
• pandas：表形式データの表示
• Matplotlib：グラフの作成

from pathlib import Path

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# 乱数を再現可能にする
SEED = 42
rng = np.random.default_rng(SEED)

# 図の保存先
output_dir = Path("least_squares_figures")
output_dir.mkdir(exist_ok=True)

# 図の基本設定
plt.rcParams["figure.figsize"] = (8, 5)
plt.rcParams["figure.dpi"] = 120

2. ダミーデータを作る

今回は、$x$と$y$の間におおよそ直線的な関係があるダミーデータを作ります。

データを生成する式は、次のとおりです。

$$
y_i=
2.5x_i + 1.0 + \varepsilon_i
$$

$\varepsilon_i$は、測定値のばらつきに相当するランダムなノイズです。

今回は、平均0、標準偏差2.0の正規分布からノイズを発生させます。

$$
\varepsilon_i
\sim
\mathcal{N}(0, 2.0^2)
$$

# データ点の数
n = 20

# 0から10までのxを20点作る
x = np.linspace(0, 10, n)

# ダミーデータを作るための真のパラメータ
true_slope = 2.5
true_intercept = 1.0
noise_std = 2.0

# ランダムなノイズ
noise = rng.normal(0, noise_std, size=n)

# y = 2.5x + 1.0 + ノイズ
y = true_slope * x + true_intercept + noise

# 表として確認する
df = pd.DataFrame({"x": x, "y": y})
display(df.head())

# 散布図
plt.figure()
plt.scatter(x, y, s=50)
plt.xlabel("x")
plt.ylabel("y")
plt.title("Dummy data")
plt.grid(alpha=0.3)
plt.tight_layout()

plt.savefig(
    output_dir / "fig01_dummy_data.png",
    dpi=160,
)
plt.show()

散布図を見ると、$x$が大きくなるにつれて$y$も大きくなる傾向があります。

ただし、すべての点が一本の直線上に並んでいるわけではありません。

これは、データにランダムなノイズを加えているためです。

ダミーデータなので、今回はデータを生成した「真の直線」が、

$$
y = 2.5x + 1.0
$$

であることを知っています。

一方、実際の研究データでは、通常、この真の関係は分かりません。

そこで、観測されたデータから、もっともよく当てはまる直線を推定します。

3. 直線フィットで使う用語

最小二乗法を理解する前に、基本的な用語を整理します。

用語	記号	意味
説明変数	$x_i$	予測に使用する変数
目的変数	$y_i$	説明したい、または予測したい観測値
回帰直線	$\hat{y}=ax+b$	データに当てはめる直線
予測値	$\hat{y}_i$	回帰直線から計算された値
残差	$e_i$	観測値と予測値の差
残差平方	$e_i^2$	残差を二乗した値
残差平方和	$RSS$	すべての残差平方の合計

観測値と予測値

データとして実際に得られた値が、観測値$y_i$です。

一方、直線から計算される値が、予測値$\hat{y}_i$です。

$$
\hat{y}_i = ax_i+b
$$

残差

観測値と予測値の差を残差と呼びます。

$$
e_i = y_i-\hat{y}_i
$$

データ点が直線より上にある場合、残差は正になります。

データ点が直線より下にある場合、残差は負になります。

厳密には、統計学でいう「誤差」と「残差」は少し異なります。

• 誤差：観測値と未知の真の値との差
• 残差：観測値と推定したモデルの予測値との差

真の値は通常分からないため、実際のデータ解析では残差を調べます。

4. なぜ残差を二乗するのか

残差をそのまま足し合わせると、正の残差と負の残差が打ち消し合ってしまいます。

たとえば、残差が$+3$と$-3$であれば、その合計は0です。

しかし、どちらの点も直線から3だけ離れています。

そこで、残差を二乗します。

$$
(+3)^2 = 9
$$

$$
(-3)^2 = 9
$$

二乗すれば、正負にかかわらず非負の値になります。

さらに、大きな残差は強く評価されます。

たとえば、残差が1から2へ2倍になると、残差平方は1から4へ4倍になります。

$$
1^2 = 1
$$

$$
2^2 = 4
$$

このため、最小二乗法は大きく外れたデータ点の影響を受けやすいという性質も持ちます。

残差を二乗する主な理由は、次のように整理できます。

1. 正と負の残差が打ち消し合わない
2. 大きな残差をより強く評価できる
3. 数学的に微分しやすく、解を求めやすい

また、誤差が独立な正規分布に従うと仮定すると、最小二乗法は最尤推定とも対応します。

5. numpy.polyfit()で直線を当てはめる

NumPyのpolyfit()を使って、ダミーデータに直線を当てはめます。

deg=1は、1次式、すなわち直線を当てはめることを意味します。

# 1次式を最小二乗法で当てはめる
slope_hat, intercept_hat = np.polyfit(
    x,
    y,
    deg=1,
)

# 各xに対する予測値
y_hat = slope_hat * x + intercept_hat

# 残差
residuals = y - y_hat

# 残差平方
squared_residuals = residuals**2

# 残差平方和
rss = np.sum(squared_residuals)

# 平均二乗誤差
mse = np.mean(squared_residuals)

# 二乗平均平方根誤差
rmse = np.sqrt(mse)

# 決定係数
sst = np.sum((y - np.mean(y))**2)
r2 = 1 - rss / sst

print(f"estimated slope     = {slope_hat:.3f}")
print(f"estimated intercept = {intercept_hat:.3f}")
print(f"RSS                 = {rss:.3f}")
print(f"MSE                 = {mse:.3f}")
print(f"RMSE                = {rmse:.3f}")
print(f"R^2                 = {r2:.3f}")

実行結果は、次のようになります。

estimated slope     = 2.596
estimated intercept = 0.452
RSS                 = 55.834
MSE                 = 2.792
RMSE                = 1.671
R^2                 = 0.957

推定された直線は、次のようになります。

$$
\hat{y}=
2.596x + 0.452
$$

ダミーデータを生成した真の直線は、

$$
y = 2.5x+1.0
$$

でした。

推定値は真の値に近いものの、完全には一致しません。

これは、データ点の数が有限であり、さらにランダムなノイズが含まれているためです。

RSS、MSE、RMSEの違い

RSSは、残差平方をすべて足し合わせた値です。

$$
RSS = \sum_{i=1}^{n} e_i^2
$$

MSEは、RSSをデータ点の数$n$で割った単純平均です。

$$
MSE = \frac{1}{n}\sum_{i=1}^{n}e_i^2
$$

RMSEは、MSEの平方根です。

$$
RMSE = \sqrt{MSE}
$$

RSSとMSEは、目的変数$y$の単位を二乗した単位になります。

一方、RMSEは$y$と同じ単位になるため、誤差の大きさを直感的に解釈しやすい指標です。

なお、ここでは予測誤差を要約する目的で$RSS/n$をMSEとしています。

切片を含む単回帰モデルから誤差分散を不偏推定する場合には、自由度を考慮して$RSS/(n-2)$を使います。両者は目的が異なるため、区別が必要です。

6. 複数の候補となる直線を比較する

最小二乗法がどの直線を選んでいるのかを見るために、3本の直線を比較します。

• Candidate A：傾き1.5、切片4.0
• Candidate B：傾き2.0、切片1.0
• Least-squares fit：最小二乗法で求めた直線

それぞれについてRSSを計算します。

candidates = [
    ("Candidate A", 1.5, 4.0),
    ("Candidate B", 2.0, 1.0),
    (
        "Least-squares fit",
        slope_hat,
        intercept_hat,
    ),
]

x_line = np.linspace(
    x.min(),
    x.max(),
    200,
)

plt.figure()
plt.scatter(
    x,
    y,
    s=50,
    label="data",
)

for name, a, b in candidates:
    candidate_y_hat = a * x + b

    candidate_rss = np.sum(
        (y - candidate_y_hat) ** 2
    )

    plt.plot(
        x_line,
        a * x_line + b,
        linewidth=2,
        label=f"{name}: RSS={candidate_rss:.1f}",
    )

plt.xlabel("x")
plt.ylabel("y")
plt.title("Candidate lines and their RSS")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()

plt.savefig(
    output_dir / "fig02_candidate_lines.png",
    dpi=160,
)
plt.show()

このコードでは、次のRSSが得られます。

直線	傾き	切片	RSS
Candidate A	1.5	4.0	352.1
Candidate B	2.0	1.0	239.9
Least-squares fit	2.596	0.452	55.8

最小二乗法で求めた直線のRSSが、もっとも小さくなっています。

ここで重要なのは、最小二乗法が「すべてのデータ点を通る直線」を探しているわけではないことです。

ノイズを含むデータでは、通常、すべての点を一本の直線で通ることはできません。

そこで、全データ点からの残差平方を合計し、その合計がもっとも小さくなる直線を選びます。

7. 傾きと切片を変えたときのRSSを見る

直線は、傾き$a$と切片$b$の2つのパラメータで決まります。

$$
\hat{y}=ax+b
$$

そこで、傾きと切片を少しずつ変えながら、すべての組み合わせについてRSSを計算してみます。

# 調べる傾きと切片の範囲
slope_values = np.linspace(
    1.0,
    4.0,
    241,
)

intercept_values = np.linspace(
    -4.0,
    6.0,
    241,
)

# 傾きと切片の格子を作る
SLOPE, INTERCEPT = np.meshgrid(
    slope_values,
    intercept_values,
)

# すべての傾き・切片について予測値を計算する
predicted = (
    SLOPE[..., np.newaxis] * x
    + INTERCEPT[..., np.newaxis]
)

# すべての組み合わせについてRSSを計算する
rss_grid = np.sum(
    (y - predicted) ** 2,
    axis=-1,
)

# 等高線として表示するRSS
levels = [
    60,
    80,
    100,
    150,
    250,
    400,
    700,
    1000,
    1500,
]

plt.figure(figsize=(8, 6))

contours = plt.contour(
    SLOPE,
    INTERCEPT,
    rss_grid,
    levels=levels,
)

plt.clabel(
    contours,
    inline=True,
    fontsize=8,
    fmt="%d",
)

# 最小二乗法で求めた位置
plt.scatter(
    slope_hat,
    intercept_hat,
    s=100,
    marker="x",
    linewidths=3,
    label="least-squares minimum",
)

plt.xlabel("slope a")
plt.ylabel("intercept b")
plt.title("RSS landscape")
plt.grid(alpha=0.2)
plt.legend()
plt.tight_layout()

plt.savefig(
    output_dir / "fig03_rss_landscape.png",
    dpi=160,
)
plt.show()

横軸は傾き$a$、縦軸は切片$b$です。

それぞれの等高線上では、RSSが同じ値になっています。

内側へ進むほどRSSは小さくなり、図中のx印が最小二乗法によって求めた位置です。

今回の最小値は、およそ次の位置にあります。

$$
a = 2.596
$$

$$
b = 0.452
$$

等高線が斜めに伸びた楕円形になっている点にも注目できます。

傾きを少し大きくして切片を少し小さくすると、データが存在する範囲内では似たような直線になることがあります。

このため、傾きと切片は完全に独立に決まるのではなく、互いに関係しながらRSSを変化させます。

最小二乗法は、この2次元の「RSSの地形」のもっとも低い場所を探していると考えられます。

8. 残差を図で確認する

観測値と予測値の差を、図に描いてみます。

plt.figure()

# 観測値
plt.scatter(
    x,
    y,
    s=50,
    label="observed y",
)

# 最小二乗法で求めた直線
plt.plot(
    x_line,
    slope_hat * x_line + intercept_hat,
    linewidth=2,
    label="fitted line",
)

# 観測値と予測値を結ぶ縦線
plt.vlines(
    x,
    y_hat,
    y,
    linewidth=1.5,
    alpha=0.8,
    label="residuals",
)

# 予測値
plt.scatter(
    x,
    y_hat,
    s=25,
    marker="x",
    label="predicted y",
)

plt.xlabel("x")
plt.ylabel("y")
plt.title("Fitted line and residuals")
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()

plt.savefig(
    output_dir / "fig04_residuals.png",
    dpi=160,
)
plt.show()

丸印が観測値$y_i$、回帰直線上のx印が予測値$\hat{y}_i$です。

両者を結ぶ縦線の長さが、残差の絶対値に相当します。

$$
e_i = y_i-\hat{y}_i
$$

データ点が直線より上にあれば、残差は正です。

データ点が直線より下にあれば、残差は負です。

ここで、通常の最小二乗法が最小化しているのは、$y$方向の縦の距離です。

データ点から直線までの最短距離、すなわち直線に垂直な距離を最小化しているわけではありません。

通常の回帰分析では、$x$を与えたときの$y$を説明または予測することを考えるためです。

$x$にも無視できない測定誤差がある場合は、直交距離回帰など、別の方法が適することがあります。

9. 各データ点の残差平方を見る

RSSは、各データ点の残差平方を足し合わせたものです。

$$
RSS=
e_1^2
+
e_2^2
+
\cdots
+
e_n^2
$$

各データ点がRSSにどのくらい寄与しているのかを、棒グラフで確認します。

plt.figure(figsize=(9, 5))

plt.bar(
    np.arange(n),
    squared_residuals,
)

plt.xlabel("data point index")
plt.ylabel("squared residual")
plt.title(
    f"Contribution of each point to RSS "
    f"(total RSS = {rss:.1f})"
)
plt.grid(axis="y", alpha=0.3)
plt.tight_layout()

plt.savefig(
    output_dir / "fig05_squared_residuals.png",
    dpi=160,
)
plt.show()

各棒の高さが、そのデータ点の残差平方です。

棒の高さをすべて足し合わせると、RSSの約55.8になります。

今回の結果では、インデックス4、すなわち5番目のデータ点の寄与が比較的大きくなっています。

残差を二乗するため、直線から大きく離れた点はRSSへ強く影響します。

これは最小二乗法の便利な性質である一方、外れ値の影響を受けやすい理由でもあります。

10. 傾きと切片はどのように計算されるのか

単回帰の最小二乗法には、傾きと切片を直接求める公式があります。

$x$と$y$の平均を、それぞれ$\bar{x}$と$\bar{y}$とすると、傾きは次のように求められます。

$$
\hat{a}=
\frac{
\sum_{i=1}^{n}
(x_i-\bar{x})(y_i-\bar{y})
}{
\sum_{i=1}^{n}
(x_i-\bar{x})^2
}
$$

切片は次の式で求められます。

$$
\hat{b}=\bar{y}-\hat{a}\bar{x}
$$

これらの式は、RSSを$a$と$b$で偏微分し、両方を0と置くことで導けます。

$$
\frac{\partial RSS}{\partial a}=0
$$

$$
\frac{\partial RSS}{\partial b}=0
$$

公式をそのままPythonで計算して、np.polyfit()の結果と比較します。

# xとyの平均
x_bar = np.mean(x)
y_bar = np.mean(y)

# 最小二乗法の公式による傾き
slope_manual = np.sum(
    (x - x_bar) * (y - y_bar)
) / np.sum(
    (x - x_bar) ** 2
)

# 最小二乗法の公式による切片
intercept_manual = (
    y_bar
    - slope_manual * x_bar
)

print(
    f"manual slope     = "
    f"{slope_manual:.6f}"
)

print(
    f"manual intercept = "
    f"{intercept_manual:.6f}"
)

print(
    "np.polyfitと一致:",
    np.allclose(
        [slope_hat, intercept_hat],
        [
            slope_manual,
            intercept_manual,
        ],
    ),
)

実行結果は、次のようになります。

manual slope     = 2.596461
manual intercept = 0.451831
np.polyfitと一致: True

np.polyfit()の結果と、公式から計算した結果が一致しました。

切片の式、

$$
\hat{b}=\bar{y}-
\hat{a}\bar{x}
$$

を変形すると、

$$
\bar{y}=
\hat{a}\bar{x}
+
\hat{b}
$$

となります。

これは、切片を含む最小二乗法の回帰直線が、必ず平均の点$(\bar{x},\bar{y})$を通ることを表しています。

11. 残差プロットを確認する

直線を描いただけでは、モデルがデータの特徴を適切に表しているかを十分に判断できません。

そこで、横軸に$x$、縦軸に残差を置いた残差プロットを作ります。

plt.figure()

plt.scatter(
    x,
    residuals,
    s=50,
)

plt.axhline(
    0,
    linestyle="--",
)

plt.xlabel("x")
plt.ylabel("residual")
plt.title("Residual plot")
plt.grid(alpha=0.3)
plt.tight_layout()

plt.savefig(
    output_dir / "fig06_residual_plot.png",
    dpi=160,
)
plt.show()

今回の残差は、0の周囲におおむねランダムに分布しています。

少なくとも、この図からは明らかな曲線的パターンは見られません。

残差プロットに次のようなパターンが現れた場合は、注意が必要です。

残差のパターン	考えられる問題
U字形や逆U字形	関係が直線ではない可能性
右に進むほど広がる	残差の分散が一定でない可能性
一定の周期がある	時系列的な依存や周期成分の可能性
一部だけ極端に大きい	外れ値や影響力の強い点の可能性

最小二乗法は、どのようなデータに対しても計算自体はできます。

しかし、直線を計算できることと、その直線が適切なモデルであることは別の問題です。

12. 最小二乗法は万能ではない

最小二乗法は非常に便利ですが、いくつか注意点があります。

外れ値の影響を受けやすい

残差を二乗するため、大きな残差は強く評価されます。

たとえば、残差が2倍になると、残差平方は4倍になります。

外れ値がある場合、回帰直線がその点に強く引っ張られることがあります。

直線関係を仮定している

今回のモデルは、

$$
\hat{y}=ax+b
$$

という直線です。

実際の関係が曲線である場合、直線を当てはめてもデータの構造を適切に表せないことがあります。

高い決定係数だけでは十分ではない

今回の決定係数は約0.957でした。

これは、このダミーデータに対して直線が比較的よく当てはまっていることを示します。

しかし、決定係数が高いことだけで、

• 因果関係がある
• モデルが正しい
• 外挿にも使える
• 残差に問題がない

とは判断できません。

散布図、残差プロット、データの取得方法、解析目的を合わせて確認する必要があります。

回帰分析だけで因果関係は示せない

$x$と$y$に直線的な関係が見られても、それだけで$x$が$y$を引き起こしたとは言えません。

第三の変数や実験条件の違いが、両方に影響している可能性があります。

13. 残差平方以外を使う方法もある

最小二乗法は、残差平方和を最小化します。

$$
\sum_{i=1}^{n}e_i^2
$$

一方、残差の絶対値を足し合わせる方法もあります。

$$
\sum_{i=1}^{n}|e_i|
$$

これはL1損失や絶対誤差と呼ばれます。

残差を二乗しないため、一般に大きな外れ値の影響を最小二乗法より受けにくくなります。

また、データ点ごとに信頼度や測定精度が異なる場合は、重み付き最小二乗法を使うことがあります。

$$
\sum_{i=1}^{n}w_i e_i^2
$$

$w_i$は、各データ点に与える重みです。

最小二乗法は重要な基本手法ですが、データの性質や解析目的に応じて、損失関数や回帰方法を選ぶ必要があります。

14. 試してみると理解が深まる変更

コードの一部を変更して、結果がどう変わるか試してみてください。

ノイズを小さくする

noise_std = 0.5

データ点が真の直線の近くに集まり、推定値も真の傾きと切片に近づきやすくなります。

ノイズを大きくする

noise_std = 5.0

データのばらつきが大きくなり、推定された直線も不安定になります。

データ点を増やす

n = 100

ノイズがランダムであれば、データ点を増やすことで真の関係を推定しやすくなります。

外れ値を追加する

yを作成した直後に、次のコードを追加します。

y[4] += 15

その後、直線フィットをやり直し、傾き、切片、RSSがどのように変化するかを確認してみてください。

曲線的なデータを作る

次のように、2次関数からデータを作ります。

y = 0.4 * x**2 + noise

直線を当てはめたあと、残差プロットにどのようなパターンが現れるかを確認してみてください。

まとめ

最小二乗法による直線フィットでは、次の直線を考えます。

$$
\hat{y}_i=ax_i+b
$$

観測値と予測値の差が残差です。

$$
e_i=y_i-\hat{y}_i
$$

最小二乗法は、残差平方和、

RSS=\sum_{i=1}^{n}\left\{y_i-(ax_i+b)\right\}^2

を最小にする傾き$a$と切片$b$を選びます。

今回のポイントをまとめると、次のようになります。

• 最小二乗法は、残差ではなく残差の二乗和を最小化する
• 残差を二乗することで、正負の打ち消しを防げる
• 大きな残差は、RSSへ強く影響する
• 傾きと切片は、RSSがもっとも小さくなる組み合わせとして選ばれる
• 最小二乗法で直線を計算できても、その直線が適切とは限らない
• 散布図だけでなく、残差プロットも確認することが重要である

np.polyfit()を使えば、直線フィット自体は簡単に実行できます。

しかし、その背後で何が最小化されているのかを理解しておくと、結果をより適切に解釈できるようになります。