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

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
$$

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

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

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

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

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()を使えば、直線フィット自体は簡単に実行できます。

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