結論

最小二乗法とは、データと式の縦のずれ(残差)の二乗の合計が最小になるように、式の係数を決める方法です。直線なら傾きと切片が公式で求まります。

この記事でわかること

  • 5点の手計算表と、傾き・切片の公式2つ
  • Python・SQL・JavaScriptの9通りの計算結果
  • 外れ値1点で傾きが2.0→0.5になる実測

2026年10月22日更新|NumPy 2.4.6・scikit-learn 1.9.1で実行確認

最小二乗法は、直線を仮に引いて残差を測り、その二乗の合計が一番小さくなる傾きと切片を探す、という考え方です。直線の場合は探さなくても、公式で一度に求まります。

flowchart TD
    A[データの点を集める] --> B[直線 y = ax + b を仮に引く]
    B --> C[各点と直線の縦のずれ<br>残差を求める]
    C --> D[残差を二乗して合計する]
    D --> E{合計が最小か}
    E -->|いいえ 傾きと切片を変える| B
    E -->|はい| F[最小二乗法の直線]

急いでいる人向けに、よく調べられることを先にまとめます。

知りたいこと 答え 詳しい節
傾き a の公式 Sxy ÷ Sxx(偏差の積の和 ÷ x の偏差の二乗和) 公式
切片 b の公式 ȳ − a × x̄(直線は平均の点を通る) 公式
5点の例の答え y = 0.9x + 1.3(R² = 0.81) 手計算
Python で1行 np.polyfit(x, y, 1) Python
行列で書くと XᵀX β = Xᵀy(正規方程式) 正規方程式
当てはまりの良さ R² = 1 − 残差の二乗和 ÷ 全体の変動 決定係数
弱いところ 外れ値。1点で傾きが2.0→0.5 外れ値
曲線に当てはめる 多項式は polyfit、それ以外は curve_fit 多項式・非線形

確認環境:Ubuntu 24.04・Python 3.11.15・NumPy 2.4.6・SciPy 1.17.1・scikit-learn 1.9.1・statsmodels 0.15.0・pandas 3.0.2・matplotlib 3.10.9・PostgreSQL 16.15・Node.js 22.22.2。5点の手計算、Python の7つの関数、SQL の集約関数、ライブラリを使わない JavaScript の実装で同じ答えになることを確かめ、外れ値・過学習・重み付き・非線形・正規方程式の精度・勾配降下法の実験を実行して、出力をそのまま載せています。関数の仕様は NumPy(2.5 のマニュアル)・SciPy(1.18 のマニュアル)・scikit-learn 1.9.1・statsmodels 0.15.0・PostgreSQL 16 の公式ドキュメントで、歴史は Wikipedia「最小二乗法」で確かめました(確認日 2026年10月22日)。

本記事の内容は、ご自由にお使いください。
ご利用の際は、出典として本ページへのリンクを記載いただけますようお願いします。

(記載例)出典:株式会社RJC「最小二乗法とは|公式・手計算の例とPythonでの求め方を実例で解説」

最小二乗法とは(考え方を図で)

この節では、最小二乗法が「何を最小にしているのか」が分かります。

最小二乗法(さいしょうにじょうほう、英語 least squares method)は、測ったデータに式を当てはめるとき、「実際の値と式の値の差」の二乗を全部足した値が最小になるように、式の係数を決める方法です。 この差を残差(ざんさ)と呼びます。

たとえば、処理するデータの件数 x(万件)と処理時間 y(秒)を5回測ったとします。この記事では、ずっとこの5点を例に使います。

回 処理件数 x(万件) 処理時間 y(秒)
1回目 1 2
2回目 2 3
3回目 3 5
4回目 4 4
5回目 5 6

点はおおよそ右上がりですが、一直線には並びません。最小二乗法で求めた直線は y = 0.9x + 1.3 です。次の図の棒の上端が測った値、線が直線で、棒の上端と線の差が残差です。

---
config:
  xyChart:
    width: 480
    height: 300
---
xychart-beta
    title "5点(棒)と直線 y = 0.9x + 1.3(線)"
    x-axis "処理件数(万件)" [1, 2, 3, 4, 5]
    y-axis "処理時間(秒)" 0 --> 7
    bar [2, 3, 5, 4, 6]
    line [2.2, 3.1, 4.0, 4.9, 5.8]

最小にする値は、残差の二乗和 S です。式で書くと次のとおりです(Σ は「全部の点について足す」という意味です)。

S(a, b) = Σ ( yi − (a·xi + b) )²
        = (点1の残差)² +(点2の残差)² + … +(点nの残差)²

a と b を少しずらすと S は必ず大きくなる、という点が最小二乗法の答えです。 5点で、b を 1.3 に固定して a だけを動かしたときの S を計算しました。

---
config:
  xyChart:
    width: 480
    height: 280
---
xychart-beta
    title "傾き a を動かしたときの残差の二乗和 S(b = 1.3)"
    x-axis "傾き a" [0.7, 0.8, 0.9, 1.0, 1.1]
    y-axis "S" 0 --> 4.5
    line [4.10, 2.45, 1.90, 2.45, 4.10]

S は a について「おわんの形」(2次式)で、底の a = 0.9 で最小の 1.90 になりました。b を動かした場合も、1.1・1.2・1.3・1.4・1.5 で S は 2.10・1.95・1.90・1.95・2.10 となり、b = 1.3 が底でした。

最小二乗法が使われる場面

1予測の式を作る

件数と処理時間、広告費と問い合わせ数など、2つの量の関係を直線にして、別の件数のときの値を見積もる。

2実験の値を決める

ばねの伸びから定数を出す、冷めていく温度から冷める速さを出すなど、式の中の定数を測定から求める。

3回帰分析の土台

統計の線形回帰(単回帰・重回帰)の係数は、ふつう最小二乗法で求める。statsmodels の結果にも「Least Squares」と出る。

4機械学習の損失

平均二乗誤差(MSE)を小さくする学習は、残差の二乗和を小さくすることと同じ。ニューラルネットワークの回帰でも使う。

Wikipedia によると、最小二乗法を初めて発表したのはルジャンドルで1805年です。ガウスは1809年に出版し、1795年には考えていたと主張したため、発明者ははっきりしていません。

最小二乗法の公式(直線の傾きと切片)

この節では、直線 y = ax + b の傾き a と切片 b を求める公式が分かります。

傾き a は「x と y の偏差の積の和」を「x の偏差の二乗の和」で割った値、切片 b は「y の平均 − a × x の平均」です。 偏差は「その値 − 平均」のことです。

x̄ = x の平均,  ȳ = y の平均
Sxx = Σ (xi − x̄)²              … x の偏差の二乗和
Sxy = Σ (xi − x̄)(yi − ȳ)       … x と y の偏差の積の和

傾き a = Sxy ÷ Sxx
切片 b = ȳ − a × x̄

Sxx と Sxy を n(点の数)で割ると、x の分散と、x と y の共分散になります。そのため、傾き a = 共分散 ÷ x の分散 と書くこともできます(割る n は分子と分母で消えるので、答えは同じです)。

偏差を使わず、合計(Σ)だけで書いた形もあります。電卓や表で計算するときに便利です。

a = ( n·Σxy − Σx·Σy ) ÷ ( n·Σx² − (Σx)² )
b = ( Σy − a·Σx ) ÷ n

公式から分かる性質が3つあります。

  • 直線は必ず平均の点 (x̄, ȳ) を通る。 切片の式を変形すると ȳ = a·x̄ + b になるためです。5点の例では (3, 4) を通ります。
  • 残差の合計は0になる。 切片を入れた直線では、上にずれた点と下にずれた点がちょうど打ち消し合います(下の手計算で 0 を確かめます)。
  • x がすべて同じ値だと傾きが決まらない。 Sxx = 0 になり、割り算ができません。

原点を通る直線(切片なし)の公式

y = ax のように、x = 0 のとき y = 0 と決まっている場合は、切片を入れずに a だけを求めます。

a = Σxy ÷ Σx²

5点では a = 69 ÷ 55 = 1.254545 でした。残差の二乗和は 3.4364 で、切片ありの直線(1.9000)より大きくなります。原点を通ると決まっている理由(物理の法則など)が無いなら、切片ありで求めます。 原点を通る直線の R² の注意は、後の「よくある間違い」の節にあります。

最小二乗法の手計算の例(5点)

この節では、5点のデータから、傾き・切片・残差・決定係数を手で計算する手順が分かります。

1

平均を出す

x の平均 x̄ = (1+2+3+4+5) ÷ 5 = 3、y の平均 ȳ = (2+3+5+4+6) ÷ 5 = 4。

x̄ = 3・ȳ = 4

2

偏差を出す

各点で x − x̄ と y − ȳ を計算します(下の表の3列目と4列目)。

平均からの差

3

Sxx と Sxy を足す

(x − x̄)² を足すと Sxx = 10、(x − x̄)(y − ȳ) を足すと Sxy = 9。

Sxx = 10・Sxy = 9

4

傾きと切片を出す

a = 9 ÷ 10 = 0.9、b = 4 − 0.9 × 3 = 1.3。直線は y = 0.9x + 1.3 です。

公式に入れる

5

残差と R² で確かめる

予測値と残差を出し、残差の二乗和 1.90 と決定係数 R² = 0.81 を求めます。

当てはまりを見る

ステップ2と3の表です。

点 (x, y) x − x̄ y − ȳ (x − x̄)² (x − x̄)(y − ȳ)
(1, 2) −2 −2 4 4
(2, 3) −1 −1 1 1
(3, 5) 0 1 0 0
(4, 4) 1 0 1 0
(5, 6) 2 2 4 4
合計 0 0 Sxx = 10 Sxy = 9

ステップ5の表です。予測値 ŷ(ワイハット)は直線の式に x を入れた値、残差 e は y − ŷ です。

x y 予測値 ŷ = 0.9x + 1.3 残差 e = y − ŷ e²
1 2 2.2 −0.2 0.04
2 3 3.1 −0.1 0.01
3 5 4.0 1.0 1.00
4 4 4.9 −0.9 0.81
5 6 5.8 0.2 0.04
合計 0.0 1.90

y の全体の変動 Syy = Σ(y − ȳ)² = 4 + 1 + 1 + 0 + 4 = 10 なので、決定係数 R² = 1 − 1.90 ÷ 10 = 0.81 です。

同じ計算を Python で書き、表の値を出力させました(ライブラリは使っていません)。

x = [1, 2, 3, 4, 5]          # 処理件数(万件)
y = [2, 3, 5, 4, 6]          # 処理時間(秒)
n = len(x)
xm, ym = sum(x) / n, sum(y) / n
sxx = sum((xi - xm) ** 2 for xi in x)
sxy = sum((xi - xm) * (yi - ym) for xi, yi in zip(x, y))
a = sxy / sxx
b = ym - a * xm
sse = sum((yi - (a * xi + b)) ** 2 for xi, yi in zip(x, y))
sst = sum((yi - ym) ** 2 for yi in y)
print(f"Sxx = {sxx}, Sxy = {sxy}")
print(f"傾き a = {a}, 切片 b = {b:.1f}")
print(f"残差の二乗和 S = {sse:.2f}, R² = {1 - sse / sst:.2f}")
print(f"残差の合計 = {sum(yi - (a*xi+b) for xi, yi in zip(x, y)):.10f}")
Sxx = 10.0, Sxy = 9.0
傾き a = 0.9, 切片 b = 1.3
残差の二乗和 S = 1.90, R² = 0.81
残差の合計 = 0.0000000000

合計だけで書いた公式でも同じ答えになりました。

Σx=15 Σy=20 Σx²=55 Σxy=69
和の形の公式: a = (nΣxy - ΣxΣy)/(nΣx² - (Σx)²) = (345 - 300)/(275 - 225) = 0.9, b = 1.3

この直線で 10万件のときを予測すると 0.9 × 10 + 1.3 = 10.3 秒です。ただし、測った範囲(1〜5万件)の外の予測は、関係がそのまま続くという仮定の上に立っています。範囲の外では、後の多項式の節のように、大きく外れることがあります。

最小二乗法の公式の導き方(偏微分)

この節では、傾きと切片の公式が、残差の二乗和を偏微分して0とおくことで出てくる流れが分かります。

S(a, b) は a と b の2次式で、下に凸(おわんの形)なので、a で偏微分した値と b で偏微分した値が両方0になる点が最小です。 偏微分は「ほかの文字を定数とみなして、1つの文字で微分する」ことです。

S(a, b) = Σ ( yi − a·xi − b )²

① b で偏微分して 0 とおく
   ∂S/∂b = −2 Σ ( yi − a·xi − b ) = 0
   → Σyi = a·Σxi + n·b
   → b = ȳ − a·x̄                       … 切片の公式

② a で偏微分して 0 とおく
   ∂S/∂a = −2 Σ xi ( yi − a·xi − b ) = 0
   → Σxi·yi = a·Σxi² + b·Σxi

③ ② に ① の b を入れて a について解く
   Σxi·yi − n·x̄·ȳ = a ( Σxi² − n·x̄² )
   左辺 = Σ(xi − x̄)(yi − ȳ) = Sxy,  右辺のかっこ = Σ(xi − x̄)² = Sxx
   → a = Sxy ÷ Sxx                     … 傾きの公式

①の「Σyi = a·Σxi + n·b」と、②の「Σxi·yi = a·Σxi² + b·Σxi」の2本の連立方程式を、正規方程式と呼びます。5点の値を入れると、次のようになります。

20 = 15a + 5b
69 = 55a + 15b
→ a = 0.9, b = 1.3

①が「残差の合計 = 0」、②が「残差と x の積の合計 = 0」を表していることも分かります。前の節の性質(残差の合計が0)は、ここから来ています。

正規方程式を行列で書く(重回帰・多項式に広げる)

この節では、正規方程式を行列で書いた形と、それを NumPy で解く方法、そして直接解くときの注意が分かります。

データを行列 X、係数をベクトル β にまとめると、正規方程式は XᵀX β = Xᵀy と1本で書けます。 X の各行は1つの点で、列は「x」と「1」(切片の分)です。

    | 1  1 |        | 2 |
    | 2  1 |        | 3 |          | a |
X = | 3  1 |,   y = | 5 |,   β =   | b |
    | 4  1 |        | 4 |
    | 5  1 |        | 6 |

XᵀX β = Xᵀy  →  β = (XᵀX)⁻¹ Xᵀy

NumPy で XᵀX と Xᵀy を作って解きました。

import numpy as np
x = np.array([1, 2, 3, 4, 5], dtype=float)
y = np.array([2, 3, 5, 4, 6], dtype=float)
X = np.column_stack([x, np.ones_like(x)])
XtX, Xty = X.T @ X, X.T @ y
print("XᵀX =", XtX.tolist(), " Xᵀy =", Xty.tolist())
print("np.linalg.solve:", np.linalg.solve(XtX, Xty))
XᵀX = [[55.0, 15.0], [15.0, 5.0]]  Xᵀy = [69.0, 20.0]
np.linalg.solve: [0.9 1.3]

XᵀX の中身は [[Σx², Σx], [Σx, n]]、Xᵀy は [Σxy, Σy] で、前の節の連立方程式そのものです。この形のまま、列を増やすと重回帰(説明変数が2つ以上)、列を x², x³ にすると多項式の当てはめになります。

正規方程式を直接解くと精度が落ちることがある

XᵀX を作ると、計算の誤差の増えやすさ(条件数)がおよそ2乗になります。 そのため、NumPy の np.linalg.lstsq(特異値分解)や scikit-learn は、XᵀX を作らずに解きます。Wikipedia も、逆行列を明示的に求めるのは通常よい方法ではないと書いています。

係数がすべて 1 の多項式の値(雑音なし)を、0〜10 の50点で当てはめ直し、求めた係数と正解(すべて 1)とのずれを比べました。

次数2: cond(X)=1.4e+02 cond(XᵀX)=2.0e+04  係数の誤差 正規方程式=4.5e-14 QR=2.4e-15 lstsq=2.4e-14
次数4: cond(X)=2.7e+04 cond(XᵀX)=7.4e+08  係数の誤差 正規方程式=3.9e-10 QR=9.1e-13 lstsq=2.6e-12
次数6: cond(X)=6.7e+06 cond(XᵀX)=4.6e+13  係数の誤差 正規方程式=1.7e-06 QR=1.2e-10 lstsq=6.8e-10
次数8: cond(X)=1.9e+09 cond(XᵀX)=3.4e+18  係数の誤差 正規方程式=3.5e-02 QR=3.8e-08 lstsq=2.1e-08
x をそろえた次数8: cond(X)=7.3e+02 cond(XᵀX)=5.3e+05
3.5%正規方程式の係数のずれ(次数8)正解 1 に対して 0.035 ずれた
2e-8lstsq の係数のずれ(次数8)XᵀX を作らない解き方
1/260万x をそろえたときの cond(X)1.9e+09 → 7.3e+02

次数8では、正規方程式で解いた係数が 0.035 ずれ、lstsq は 2.1e-08 のずれで済みました。 x を「平均0・標準偏差1」にそろえると、条件数は 1.9e+09 から 7.3e+02 まで下がります。数式で理解するときは正規方程式、プログラムで解くときは lstsq などの関数、と使い分けます。行列の作り方は、関連記事「NumPyの使い方 基本一覧|配列の作成・演算・形状変更」で扱っています。

なぜ二乗するのか(絶対値ではない理由)

この節では、残差を「二乗」で足す理由と、絶対値で足す方法との違いが分かります。

残差をそのまま足すと、プラスとマイナスが打ち消し合い、どんな直線でも合計0にできてしまいます(平均の点を通る直線なら、どれも合計0です)。そこで、符号をなくすために二乗か絶対値を使います。

残差の符号をなくす2つの方法

絶対値の和(最小絶対偏差)外れ値に強い。ただし微分できない点があり、公式1本では解けない(線形計画法などで解く)
→
二乗の和(最小二乗法)微分して0とおくと公式で一度に解ける。答えは1つに決まる。大きなずれを重く見るので外れ値に弱い

誤差が正規分布に従うと考えられるとき、最小二乗法の答えは最尤推定(もっともありそうな値)と一致する

二乗を使う理由は、次の4つにまとめられます。

理由 中身
計算が簡単 2次式なので偏微分して0とおけば、公式(正規方程式)で解ける
答えが1つに決まる x がすべて同じ値でなければ、最小になる a・b は1組だけ
確率の考えと合う 誤差が正規分布なら、最小二乗法の答えは最尤推定と同じ
線形の推定で最良 誤差の平均が0・分散が同じ・互いに無相関なら、線形で偏りのない推定の中で分散が最小(ガウス=マルコフの定理)

その代わり、ずれが2倍になると罰は4倍になるので、遠く離れた1点に直線が強く引っぱられます。 実際にどれくらい動くかは、後の「外れ値の影響」の節で実測します。

Pythonで最小二乗法を計算する(NumPy・scikit-learn・SciPy・statsmodels)

この節では、Python の主なライブラリで最小二乗法の直線を求める書き方と、同じ5点での出力が分かります。

1行で済ませるなら np.polyfit(x, y, 1)、機械学習の流れに乗せるなら scikit-learn の LinearRegression、標準誤差や p 値まで見るなら statsmodels の OLS が向いています。

import numpy as np
from numpy.polynomial import Polynomial
from scipy import stats
from sklearn.linear_model import LinearRegression
import statsmodels.api as sm

x = np.array([1, 2, 3, 4, 5], dtype=float)
y = np.array([2, 3, 5, 4, 6], dtype=float)

# 1) np.polyfit(次数1): 高い次数の係数から順に返る
a, b = np.polyfit(x, y, 1)

# 2) Polynomial.fit(NumPy が新しいコードに勧める API)
p = Polynomial.fit(x, y, 1)
print(p.convert().coef)          # 低い次数から [b, a]

# 3) np.linalg.lstsq(行列 X を自分で作る)
X = np.column_stack([x, np.ones_like(x)])
coef, ssr, rank, sv = np.linalg.lstsq(X, y)

# 4) scipy.stats.linregress
r = stats.linregress(x, y)

# 5) scikit-learn(X は2次元: 行=データ、列=特徴量)
model = LinearRegression().fit(x.reshape(-1, 1), y)
print(model.coef_, model.intercept_, model.score(x.reshape(-1, 1), y))

# 6) statsmodels(定数項は自分で足す)
res = sm.OLS(y, sm.add_constant(x)).fit()

実行した結果です。

numpy 2.4.6 | scipy 1.17.1 | scikit-learn 1.9.1 | statsmodels 0.15.0

[1] np.polyfit       : a = 0.9  b = 1.3
[2] Polynomial.fit   : [1.3 0.9]   ← [b, a] の順
    convert() しないと: [4.  1.8]  window= [-1.  1.]  domain= [1. 5.]
[3] np.linalg.lstsq  : [0.9 1.3]  残差二乗和 = [1.9]  rank = 2  特異値 = [7.6912 0.9194]
[4] XᵀX = [[55.0, 15.0], [15.0, 5.0]]  Xᵀy = [69.0, 20.0]
    np.linalg.solve   : [0.9 1.3]
[5] stats.linregress : slope=0.9 intercept=1.3 rvalue=0.9000 pvalue=0.0374 stderr=0.2517 intercept_stderr=0.8347
[6] LinearRegression : coef_ = [0.9]  intercept_ = 1.3  score(R²) = 0.81
    predict(x=10)    : [10.3]
[7] statsmodels OLS  : [1.3 0.9]  R² = 0.81  調整済み R² = 0.7467

どの方法でも a = 0.9、b = 1.3 になりました。 違いは、係数の並び順と、一緒に返ってくる情報です。

関数 返る係数の順 一緒に分かること 注意
np.polyfit(x, y, 1) [a, b](高い次数から) full=True で残差・ランク、cov=True で共分散 公式は新しいコードには Polynomial.fit を勧めている
Polynomial.fit(x, y, 1) .convert().coef で [b, a] 数値が安定しやすい convert() を忘れると、x を −1〜1 に写した係数のまま
np.linalg.lstsq(X, y) X の列の順 残差二乗和・ランク・特異値 切片の列(1 の列)を自分で足す
stats.linregress(x, y) slope・intercept r・p 値・標準誤差 説明変数は1つだけ
LinearRegression() coef_・intercept_ score() で R²、predict() で予測 X は2次元(reshape(-1, 1))
sm.OLS(y, X) X の列の順 標準誤差・t 値・p 値・信頼区間・AIC 切片は add_constant で足さないと入らない

Polynomial.fit の convert() しない係数 [4, 1.8] は、x の範囲 1〜5 を −1〜1 に写した変数での係数です(x = 3 のとき 4、x の幅2あたり1.8増える)。NumPy のドキュメントは、この写し方のおかげで数値が安定しやすいと説明しています。

scikit-learn の LinearRegression は、ドキュメントによると内部で scipy.linalg.lstsq を呼ぶ通常の最小二乗法です。positive=True を付けると、係数を0以上に制限した最小二乗法(scipy.optimize.nnls)になります。

CSV を読んで直線を引く(pandas・matplotlib)

CSV のデータから直線を求め、グラフに残差の線も描くコードです。実行すると y = 0.900x + 1.300 と表示され、点・直線・残差(点線)を描いた画像が保存されました。

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")                      # 画面の無い環境でファイルに書き出す
import matplotlib.pyplot as plt

df = pd.read_csv("data.csv")               # 列は「件数」「時間」
a, b = np.polyfit(df["件数"], df["時間"], 1)
print(f"y = {a:.3f}x + {b:.3f}")

xs = np.linspace(df["件数"].min(), df["件数"].max(), 100)
fig, ax = plt.subplots(figsize=(5, 3.5))
ax.scatter(df["件数"], df["時間"], label="data")
ax.plot(xs, a * xs + b, label=f"y = {a:.2f}x + {b:.2f}")
ax.vlines(df["件数"], df["時間"], a * df["件数"] + b, linestyles="dotted", label="residual")
ax.set_xlabel("x"); ax.set_ylabel("y"); ax.legend()
fig.savefig("fit.png", dpi=120, bbox_inches="tight")

CSV の読み込みと列の取り出し方は、関連記事「pandasの使い方 逆引き81選|pandas 3.0の実行結果つき」で扱っています。グラフの日本語が文字化けするときの直し方は、関連記事「matplotlibの使い方と日本語文字化け対策|OS別・Colab・3.12対応」で扱っています。セルごとに試しながら計算する環境は、関連記事「JupyterLabの使い方|ショートカット・マジックコマンド一覧」で扱っています。

Python以外で最小二乗法を計算する(SQL・JavaScript)

この節では、データベースの SQL と、ライブラリを使わない JavaScript で同じ直線を求める方法が分かります。

PostgreSQL の regr_slope・regr_intercept

PostgreSQL には、最小二乗法の直線を求める集約関数があります。 引数は (Y, X) の順で、従属変数(y)が先です。どちらかが NULL の行は計算から除かれます(PostgreSQL 16 のドキュメント)。

CREATE TABLE 計測 (件数 numeric, 時間 numeric);
INSERT INTO 計測 VALUES (1,2),(2,3),(3,5),(4,4),(5,6);
SELECT regr_slope(時間, 件数)     AS 傾き,
       regr_intercept(時間, 件数) AS 切片,
       regr_r2(時間, 件数)        AS 決定係数,
       corr(時間, 件数)           AS 相関係数,
       regr_count(時間, 件数)     AS 件数
FROM 計測;
 傾き | 切片 | 決定係数 | 相関係数 | 件数 
------+------+----------+----------+------
  0.9 |  1.3 |     0.81 |      0.9 |    5
(1 row)

途中の値(Sxx・Sxy・Syy・平均)も関数で出せます。手計算の表と同じ値でした。

 sxx | sxy | syy | x平均 | y平均 
-----+-----+-----+-------+-------
  10 |   9 |  10 |     3 |     4
(1 row)

ウィンドウ関数の形(OVER ())で使うと、行ごとの予測値と残差を一緒に出せます。

SELECT 件数, 時間,
       regr_slope(時間, 件数) OVER () * 件数 + regr_intercept(時間, 件数) OVER () AS 予測,
       時間 - (regr_slope(時間, 件数) OVER () * 件数 + regr_intercept(時間, 件数) OVER ()) AS 残差
FROM 計測 ORDER BY 件数;
 件数 | 時間 | 予測 |         残差         
------+------+------+----------------------
    1 |    2 |  2.2 | -0.20000000000000018
    2 |    3 |  3.1 | -0.10000000000000009
    3 |    5 |    4 |                    1
    4 |    4 |  4.9 |  -0.9000000000000004
    5 |    6 |  5.8 |  0.20000000000000018
(5 rows)

regr 系の関数は double precision(浮動小数点数)で計算するので、残差の末尾に小さな誤差が出ます。表示するときは round(…::numeric, 2) で丸めます。

JavaScript(ライブラリなし)

公式をそのまま書いた関数です。Node.js 22 で実行しました。

// 最小二乗法で直線 y = a*x + b を求める(ライブラリなし)
export function linearFit(xs, ys) {
  if (xs.length !== ys.length || xs.length < 2) throw new Error('点が2つ以上必要です');
  const n = xs.length;
  const xm = xs.reduce((s, v) => s + v, 0) / n;
  const ym = ys.reduce((s, v) => s + v, 0) / n;
  let sxx = 0, sxy = 0, syy = 0;
  for (let i = 0; i < n; i++) {
    const dx = xs[i] - xm, dy = ys[i] - ym;
    sxx += dx * dx; sxy += dx * dy; syy += dy * dy;
  }
  if (sxx === 0) throw new Error('x がすべて同じ値なので傾きが決まりません');
  const a = sxy / sxx;
  const b = ym - a * xm;
  const sse = xs.reduce((s, x, i) => s + (ys[i] - (a * x + b)) ** 2, 0);
  return { a, b, r2: 1 - sse / syy };
}

console.log(linearFit([1, 2, 3, 4, 5], [2, 3, 5, 4, 6]));
try { linearFit([3, 3, 3], [1, 2, 3]); } catch (e) { console.log('Error:', e.message); }
{ a: 0.9, b: 1.2999999999999998, r2: 0.8099999999999999 }
Error: x がすべて同じ値なので傾きが決まりません

b が 1.2999999999999998 になるのは、0.9 × 3 などを2進数の小数で計算するときの丸め誤差です。画面に出すときは toFixed(3) などで丸めます。平均を引いてから掛ける書き方(偏差の形)は、「n·Σxy − Σx·Σy」のように大きな数どうしを引く書き方より、桁落ちしにくくなります。

決定係数R²と当てはまりの見方(標準誤差・p値)

この節では、求めた直線がどれくらいデータに合っているか、係数がどれくらい確かかを読む方法が分かります。

決定係数 R² は「y のばらつきのうち、直線で説明できた割合」で、1 − 残差の二乗和 ÷ 全体の変動 で求めます。 5点では 1 − 1.90 ÷ 10 = 0.81 で、y のばらつきの81%を直線で説明できています。説明変数が1つの直線では、R² は相関係数 r の2乗と同じです(r = 0.9 → 0.81)。

statsmodels の summary() は、係数ごとの確かさもまとめて表示します。5点の出力の主要部分です。

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.810
Model:                            OLS   Adj. R-squared:                  0.747
Method:                 Least Squares   F-statistic:                     12.79
No. Observations:                   5   AIC:                             13.35
Df Residuals:                       3   BIC:                             12.57
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const(b)       1.3000      0.835      1.558      0.217      -1.356       3.956
x(a)           0.9000      0.252      3.576      0.037       0.099       1.701
==============================================================================
出力の項目 意味 5点の値
R-squared 決定係数 R² 0.810
Adj. R-squared 自由度調整済み R²。説明変数を増やすだけでは上がらないように補正 0.747
std err 係数の標準誤差(係数のぶれの大きさ) 傾き 0.252
t・P>|t| 「本当は係数が0」と仮定したときの t 値と p 値 傾き 3.576・0.037
[0.025 0.975] 係数の95%信頼区間 傾き 0.099〜1.701

傾きの標準誤差は、手でも出せます。残差の二乗和を自由度(点の数 − 係数の数 = 3)で割り、Sxx で割って平方根をとります。

s² = 1.90 ÷ 3 = 0.6333
傾きの標準誤差 = √(0.6333 ÷ 10) = 0.2517
t = 0.9 ÷ 0.2517 = 3.576
0.81決定係数 R²相関係数 0.9 の2乗
0.037傾きの p 値5点なので信頼区間は広い
0.10〜1.70傾きの95%信頼区間点を増やすと狭くなる

傾きは 0.9 と求まっても、5点では95%信頼区間が 0.10〜1.70 と広く、「1万件増えると0.1〜1.7秒増える」くらいしか言えません。 R² が高いことと、係数が確かなことは別の話です。statsmodels は、5点では正規性の検定(omnibus)ができないという警告(omni_normtest is not valid with less than 8 observations)も出しました。

R² は、学習に使っていないデータで測ると0より小さくなることもあります。scikit-learn のドキュメントも、score() の R² は最大 1.0 で、平均を返すだけのモデルより悪ければ負になると書いています。

外れ値の影響(1点で直線がどれだけ動くか)

この節では、外れ値が最小二乗法の直線をどれくらい動かすかと、外れ値に強い方法の違いが、実測で分かります。

5点のうち1点を打ち間違えた場合

x = 5 の 6 を 16 と打ち間違えて計算しました。

データ 傾き a 切片 b R² x = 10 の予測
元の5点 0.900 1.300 0.8100 10.30
x = 5 を 16 に打ち間違え 2.900 −2.700 0.6469 26.30

1点の打ち間違いで傾きが3倍以上になり、10万件の予測は10.3秒から26.3秒に変わりました。

外れ値の位置で影響が変わる(30点)

本当の関係が y = 2x + 1 のデータ30点(雑音つき)に、外れ値を1点足しました。1つは x の真ん中(x = 5、y = 40)、もう1つは x の範囲の外(x = 25、y = 0)です。最小二乗法と、外れ値に強い4つの方法(scikit-learn)で傾きと切片を求めました。

from sklearn.linear_model import (LinearRegression, HuberRegressor, RANSACRegressor,
                                  TheilSenRegressor, QuantileRegressor)
X = x.reshape(-1, 1)
LinearRegression().fit(X, y)                       # 最小二乗法
HuberRegressor().fit(X, y)                         # 小さい残差は二乗、大きい残差は絶対値
QuantileRegressor(quantile=0.5, alpha=0).fit(X, y) # 絶対値の和(最小絶対偏差)
TheilSenRegressor(random_state=0).fit(X, y)        # 2点ずつの傾きの中央値
RANSACRegressor(random_state=0).fit(X, y)          # 外れ値を除いた点で当てはめ直す
--- 30点(外れ値なし)
  最小二乗法(OLS)         a=  2.011  b=  0.796
  Huber 回帰           a=  2.008  b=  0.786
  最小絶対偏差(LAD)        a=  2.037  b=  0.512
  Theil-Sen          a=  2.009  b=  0.554
  RANSAC             a=  2.011  b=  0.796
--- 真ん中(x=5)に y=40 を1点足す
  最小二乗法(OLS)         a=  2.011  b=  1.736
  Huber 回帰           a=  2.007  b=  0.854
  最小絶対偏差(LAD)        a=  2.037  b=  0.512
  Theil-Sen          a=  2.016  b=  0.568
  RANSAC             a=  2.011  b=  0.796
--- 右端の外(x=25)に y=0 を1点足す
  最小二乗法(OLS)         a=  0.500  b=  7.675
  Huber 回帰           a=  1.884  b=  1.347
  最小絶対偏差(LAD)        a=  1.968  b=  0.782
  Theil-Sen          a=  1.967  b=  0.735
  RANSAC             a=  2.011  b=  0.796

x=25 に外れ値1点を足したときの傾き本当の傾きは 2.0。棒の長さは 2.2 を全幅とする

外れ値なし OLS2.011
最小二乗法0.500
Huber 回帰1.884
Theil-Sen1.967
最小絶対偏差1.968
RANSAC2.011

読み取れることは次の3つです。

  • 真ん中の外れ値は、傾きをほとんど変えず(2.011)、切片だけを 0.796 → 1.736 に動かしました。 x の平均に近い点は、てこの支点の近くにあるので、傾きへの影響が小さくなります。
  • 範囲の外の外れ値は、最小二乗法の傾きを 2.011 → 0.500 にしました。 てこの支点から遠い点は、直線を強く回します。
  • 外れ値に強い方法は、どれも 1.88〜2.01 に収まりました。外れ値が無いときは、最小二乗法とほぼ同じ答えでした。

外れ値を見つける(てこ比・Cook の距離)

statsmodels で、足した点の「てこ比(leverage。x が平均から離れている度合い)」と「Cook の距離(その点を除くと係数がどれだけ変わるか)」を出しました。

res = sm.OLS(y, sm.add_constant(x)).fit()
infl = res.get_influence()
lev = infl.hat_matrix_diag            # てこ比
cook = infl.cooks_distance[0]         # Cook の距離
stud = infl.resid_studentized_external
足した点 てこ比(平均 0.065) Cook の距離 スチューデント化残差
真ん中(x = 5、y = 40) 0.032 0.466 27.48
範囲の外(x = 25、y = 0) 0.605 21.549 −30.77

どちらの場合も、Cook の距離が目安の 4/n(= 0.129)を超えたのは、足した1点だけでした。範囲の外の点は、てこ比が平均の約9倍で、Cook の距離も46倍大きく、直線への影響の大きさが数値に出ています。 見つけた点は、すぐに消すのではなく、打ち間違い・測り間違いかどうかを元の記録で確かめます。

多項式の当てはめと過学習

この節では、曲線(多項式)に最小二乗法を使う方法と、次数を上げすぎると起きる過学習が、数値で分かります。

多項式 y = c0 + c1·x + c2·x² + … も、係数については1次式なので、同じ最小二乗法(正規方程式)で解けます。 NumPy では np.polyfit(x, y, 次数) か Polynomial.fit(x, y, 次数) です。

本当の関係が 2次式 y = 0.5x² − 2x + 3 のデータを10点(雑音つき)作り、次数を変えて当てはめました。学習に使った10点の誤差(RMSE:残差の二乗平均の平方根)と、新しく作った200点での誤差を比べています。

from numpy.polynomial import Polynomial
for deg in [1, 2, 3, 5, 7, 9]:
    p = Polynomial.fit(x_tr, y_tr, deg).convert()
    print(deg, rmse(p(x_tr), y_tr), rmse(p(x_te), y_te), np.abs(p.coef).max())
次数 学習の10点の RMSE 新しい200点の RMSE 係数の絶対値の最大
1 1.838 1.683 1.1
2 0.694 0.922 3.5
3 0.665 0.946 4.0
5 0.360 1.624 23.6
7 0.349 1.843 133.1
9 0.000 646.563 100248.7

次数を上げるほど学習の点への誤差は小さくなり、9次では10点を全部通りました(残差の最大 5.57e-09)。ところが新しい点への誤差は 646.563 まで増えました。 これが過学習(学習に使ったデータの雑音まで覚えてしまい、新しいデータで外れる状態)です。新しい点の誤差が最小だったのは、本当の関係と同じ2次でした。2次で求めた係数は [3.513, −2.76, 0.625](本当は [3, −2, 0.5])です。

範囲の外では差がさらに大きくなります。x = 6.5(学習の点は 0.29〜5.58)の予測は、2次式が 11.96(本当は 11.125)、9次式は −28944.9 でした。

flowchart TD
    A[次数が低すぎる<br>学習も新しい点も外れる] -->|次数を上げる| B[ちょうどよい<br>新しい点の誤差が最小]
    B -->|さらに上げる| C[次数が高すぎる<br>学習だけ合う 過学習]

過学習を防ぐ方法の1つが、係数が大きくなることに罰を足すリッジ回帰(L2 正則化)です。最小にする値が「残差の二乗和 + alpha × 係数の二乗和」になります。同じ9次の多項式(x をそろえてから)で alpha を変えました。

alpha 学習の RMSE 新しい点の RMSE
0.000001 0.362 1.446
0.01 0.533 1.054
0.1 0.765 0.937
1.0 0.933 1.061

alpha = 0.1 で、新しい点の誤差は 0.937 になり、2次式(0.922)に近づきました。次数や alpha は、学習に使わないデータでの誤差(交差検証など)を見て決めます。

重み付き最小二乗法(WLS)

この節では、測定ごとに精度が違うときに使う重み付き最小二乗法と、NumPy・statsmodels での重みの渡し方が分かります。

重み付き最小二乗法は、精度の高い点の残差を重く、低い点の残差を軽く数えて、二乗和を最小にする方法です。 誤差の標準偏差 σ(シグマ)が分かっているなら、二乗した残差に 1/σ² を掛けて足します。

S = Σ wi ( yi − a·xi − b )²,   wi = 1 / σi²

ライブラリごとに、重みの渡し方が違います。

関数 渡すもの 意味
np.polyfit(x, y, 1, w=…) w = 1/σ 二乗する前の残差に掛ける(ドキュメントの説明)
sm.WLS(y, X, weights=…) weights = 1/σ² 分散の逆数に比例する重み
LinearRegression().fit(X, y, sample_weight=…) 1/σ² 二乗した残差に掛ける
curve_fit(f, x, y, sigma=…) sigma = σ 残差を σ で割ってから二乗する

5点のうち x = 3 の測定だけ誤差が10倍(σ = 2.0、ほかは 0.2)だったとして計算しました。

重みなし        : [0.9 1.3]
w=1/sigma      : [0.9    1.0531]
statsmodels WLS: [0.9    1.0531]  ← weights は 1/sigma²

np.polyfit に 1/σ、statsmodels に 1/σ² を渡すと、同じ答えになりました。 傾きが 0.9 のままなのは、x = 3 が x の平均の位置にあり、もともと傾きには効いていないためです。精度の低い 5 という値の影響が減り、切片が 1.3 → 1.0531 に下がりました。

誤差が x に比例して大きくなるデータ(200点、本当は a = 3、b = 2)では、次のようになりました。同じ条件でデータを1000回作り直し、推定値のばらつきも測っています。

方法 傾き(標準誤差) 切片(標準誤差) 1000回の傾きのばらつき
最小二乗法(重みなし) 3.0184(0.0552) 1.8571(0.3361) 0.0538
重み付き(1/σ²) 3.0010(0.0383) 1.9466(0.1198) 0.0358

重みを付けると、傾きのばらつきは約3分の2、切片のばらつきは約半分(0.2242 → 0.1129)になりました。 重みなしでも平均的には正しい値の近くに来ますが、ぶれが大きくなります。

非線形最小二乗法(SciPyのcurve_fit)

この節では、指数関数や正弦波のように、係数について1次式でないモデルに最小二乗法を使う方法と、つまずきやすい点が分かります。

係数が式の中で掛け算や指数の肩に入るモデルは、正規方程式で一度に解けないので、初期値から少しずつ係数を直して残差の二乗和を下げていきます(反復法)。 SciPy の curve_fit は、範囲の制限が無いとき Levenberg-Marquardt 法('lm')、範囲を付けると 'trf' を使います(SciPy のドキュメント)。

冷めていく温度 T(t) = C + A·exp(−k·t)(本当は A = 70、k = 0.12、C = 22)に、雑音を足した16点で当てはめました。

import numpy as np
from scipy.optimize import curve_fit

model = lambda t, A, k, C: C + A * np.exp(-k * t)
popt, pcov = curve_fit(model, t, T, p0=[60, 0.1, 20])   # p0 は初期値
perr = np.sqrt(np.diag(pcov))                            # 係数の標準誤差
p0=[60,0.1,20] : [71.2017  0.1156 21.0311]  標準誤差 [0.598  0.0027 0.4748]  本当は [70, 0.12, 22]
関数の評価回数 nfev = 17  mesg: Both actual and predicted relative reductions in the sum of squares
  are at most 0.000000

初期値で答えが変わる(局所解)

p0 を省くと、初期値はすべて 1 になります。 上の温度のモデルでは省いても同じ答えになりましたが、正弦波 y = A·sin(ωt + φ)(本当は A = 2、ω = 3、φ = 0.5)では、エラーを出さずに別の答えで止まりました。

sine = lambda t, A, w, ph: A * np.sin(w * t + ph)
for p0 in [None, [2, 3, 0]]:
    p, _ = curve_fit(sine, ts, ys, p0=p0)
    r = ys - sine(ts, *p)
    print(f"正弦波 p0={p0}: A, ω, φ = {p.round(4)}  残差二乗和={r @ r:.3f}")
正弦波 p0=None: A, ω, φ = [159.3558   0.003    3.1358]  残差二乗和=108.446  (本当は 2, 3, 0.5)
正弦波 p0=[2, 3, 0]: A, ω, φ = [2.0137 2.9938 0.5032]  残差二乗和=1.889  (本当は 2, 3, 0.5)

反復法は、初期値の近くの「谷」(局所解)で止まることがあります。 残差の二乗和を比べると、正しい答えは 1.889、止まった答えは 108.446 でした。グラフに描いて当てはまりを目で確かめる、初期値を変えて何度か解き残差の二乗和を比べる、といった確認をします。

計算の回数の上限に届くと、次のエラーになりました。

RuntimeError: Optimal parameters not found: Number of calls to function has reached maxfev = 5.

対数をとって直線にする方法との違い

y = A·exp(−k·t) は、両辺の対数をとると log y = log A − k·t の直線になります。この形で直線に当てはめる方法と、curve_fit を比べました(本当は A = 50、k = 0.15。雑音は y に直接足しています)。

本当は A=50, k=0.15 / 対数をとって直線: A=106.581, k=0.2334(使えた点 15/16)/ curve_fit: A=48.949, k=0.1464

対数をとる方法は A = 106.581 と大きくずれ、0 以下になった1点は対数がとれずに使えませんでした。 対数をとると、y の小さい点(雑音の割合が大きい点)の残差が大きく数えられるためです。雑音が y に足される形なら、curve_fit でそのまま当てはめるほうが本当の値に近くなりました。

外れ値に強い非線形の当てはめ(least_squares)

scipy.optimize.least_squares は、loss で残差の数え方を変えられます。温度のデータの1点に +25 の外れ値を入れました。

from scipy.optimize import least_squares
res = lambda p: model(t, *p) - T_out
least_squares(res, x0=[60, 0.1, 20], loss="soft_l1", f_scale=1.0)
least_squares loss=linear  : [76.8124  0.0926 16.6912]  status=2
least_squares loss=soft_l1 : [71.4896  0.1134 20.7137]  status=2
least_squares loss=huber   : [71.4047  0.1143 20.8509]  status=2

linear(普通の最小二乗法)は k が 0.0926 までずれ、soft_l1 と huber は外れ値の無いときの答え(71.2017・0.1156・21.0311)に近い値になりました。ドキュメントは、soft_l1 を外れ値に強い当てはめの最初の選択肢として挙げています。

最小二乗法と回帰分析・機械学習の関係

この節では、最小二乗法・回帰分析・機械学習の言葉の関係と、使い分けが分かります。

最小二乗法

係数を決める「計算の方法」

残差の二乗和を最小にする

直線・多項式・非線形の式に使える

回帰分析

変数の関係を調べる「分析の枠組み」

係数の推定・検定・予測・診断まで含む

係数の推定にふつう最小二乗法を使う

回帰分析は「y を x で説明する式を作って調べる」分析全体の名前で、その係数を決める計算に最小二乗法を使います。 回帰分析の中でも、最尤法・ベイズ推定・頑健な推定など、最小二乗法以外の方法を使うことがあります。

名前 何をするか 最小二乗法との関係
単回帰分析 説明変数1つの直線 この記事の5点の例そのもの
重回帰分析 説明変数が2つ以上 正規方程式の X の列を増やす
多項式回帰 x², x³ の列を足す 係数について1次式なので同じ方法
リッジ回帰 係数の二乗和に罰を足す 最小二乗法+正則化
ロジスティック回帰 「起きる確率」を予測 最小二乗法ではなく最尤法で解く
ニューラルネットワークの回帰 層を重ねた非線形の式 平均二乗誤差を勾配降下法で小さくする

重回帰(説明変数が2つ)

処理件数 x1 と同時に動くジョブ数 x2 から処理時間を予測する形(本当は y = 0.8·x1 + 0.5·x2 + 1.0、40点)です。

X = np.column_stack([x1, x2])
m = LinearRegression().fit(X, y2)
A = np.column_stack([X, np.ones(n)])
np.linalg.lstsq(A, y2)[0]
LinearRegression: a1, a2 = [0.8126 0.5173]  b = 0.8579  R² = 0.9529
np.linalg.lstsq : [0.8126 0.5173 0.8579]
statsmodels     : [0.8579 0.8126 0.5173]  調整済み R² = 0.9503

3つの方法で同じ係数になりました。statsmodels は定数項を先頭の列に足したので、b が先に並んでいます。

勾配降下法でも同じ答えにたどり着く

ニューラルネットワークの学習は、損失(平均二乗誤差 MSE)を傾き(勾配)の逆向きに少しずつ下げる勾配降下法で行います。5点の直線を、公式を使わずに勾配降下法(学習率 0.02、a = b = 0 から開始)で求めました。

a, b, lr = 0.0, 0.0, 0.02
for step in range(1, 5001):
    e = (a * x5 + b) - y5          # 予測 − 実測
    loss = np.mean(e ** 2)         # 平均二乗誤差(MSE)
    a -= lr * 2 * np.mean(e * x5)  # ∂loss/∂a
    b -= lr * 2 * np.mean(e)       # ∂loss/∂b
  step    1: a=0.5520 b=0.1600 MSE=18.0000
  step   10: a=1.1506 b=0.3876 MSE=0.5336
  step  100: a=1.0371 b=0.8049 MSE=0.4252
  step  500: a=0.9091 b=1.2672 MSE=0.3802
  step 1000: a=0.9003 b=1.2989 MSE=0.3800
  step 2000: a=0.9000 b=1.3000 MSE=0.3800
  step 5000: a=0.9000 b=1.3000 MSE=0.3800
  最小二乗法の答え: a=0.9 b=1.3 MSE=0.3800

2000回で、公式と同じ a = 0.9、b = 1.3 にたどり着きました(MSE 0.38 = 残差の二乗和 1.90 ÷ 5)。 直線なら公式で一度に解けますが、層を重ねたニューラルネットワークは公式で解けないので、この反復で係数(重み)を決めます。ニューラルネットワークの学習の流れは、関連記事「ニューラルネットワークとは|仕組み・学習・AttentionをNumPyで実演」で扱っています。

どの方法を選ぶか(判断表)

状況 選ぶ方法 Python の関数
2つの量の関係を直線で表したい 最小二乗法(単回帰) np.polyfit(x, y, 1)・stats.linregress
係数の確かさ(p 値・信頼区間)も知りたい 最小二乗法+統計量 sm.OLS(y, sm.add_constant(X))
説明変数が2つ以上 重回帰 LinearRegression・sm.OLS
曲がった関係(多項式で表せる) 多項式回帰 Polynomial.fit(x, y, 次数)
測定ごとに精度が違う 重み付き最小二乗法 polyfit(w=1/σ)・sm.WLS
指数・正弦波など式の形が決まっている 非線形最小二乗法 curve_fit・least_squares
外れ値が混ざっている 頑健な回帰 HuberRegressor・TheilSenRegressor・RANSACRegressor
x と y の両方に同じくらい誤差がある 直交回帰(全最小二乗法) 特異値分解で求める(下の節)
次数が高い・変数が多く過学習する リッジ回帰など Ridge

最小二乗法のよくある間違いとエラーの直し方

この節では、最小二乗法でよくある考え違いと、Python で出るエラーの原文と直し方が分かります。

y を x で回帰した直線と、x を y で回帰した直線は違う

最小二乗法は「縦(y の方向)のずれ」を最小にするので、x と y を入れ替えると別の直線になります。 5点で x を y で回帰すると x = 0.9y − 0.6 になり、これを y について解くと y = 1.1111x + 0.6667 でした。y = 0.9x + 1.3 とは違います。

点から直線への垂直な距離を最小にする直交回帰(全最小二乗法)は、特異値分解で求められます。5点では a = 1.0000、b = 1.0000 で、3つの直線の傾きは 0.9・1.0・1.111 と並びました。

xc, yc = x - x.mean(), y - y.mean()
_, _, vt = np.linalg.svd(np.column_stack([xc, yc]))
nx, ny = vt[-1]                 # 直線の法線ベクトル
a_tls = -nx / ny
b_tls = y.mean() - a_tls * x.mean()

予測したい量を y に置き、x の測定の誤差が小さいときは通常の最小二乗法を使います。

原点を通る直線の R² は、ライブラリで値が違う

5点を原点を通る直線(a = 1.254545)で当てはめ、R² を出しました。

LinearRegression(fit_intercept=False).coef_ = [1.254545]  score(R²) = 0.6564
statsmodels(定数項なし)R² = 0.9618 ← 中心化しない R²(Σy² を分母に使う)

同じ直線なのに、R² が 0.6564 と 0.9618 に分かれました。 statsmodels は定数項が無いとき、全体の変動に Σ(y − ȳ)² ではなく Σy² を使います(1 − 3.4364 ÷ 90 = 0.9618)。切片ありのモデルと R² を比べるときは、同じ定義でそろえます。

チェックリスト

  • 1散布図を描いてから当てはめたか曲がった関係に直線を当てはめても、R² だけでは気づきにくい
  • 2外れ値を確かめたか範囲の外の1点で傾きが 2.0 → 0.5 になった。Cook の距離で探し、元の記録を見る
  • 3測った範囲の外を予測していないか9次式は x = 6.5 で −28944.9 を返した。範囲の外は関係が続く保証がない
  • 4係数の並び順を確かめたかpolyfit は [a, b]、Polynomial と statsmodels は [b, a] の順
  • 5相関と因果を分けたか直線が引けても、x を変えれば y が変わるとは限らない

エラー・警告の原文と直し方

どれも、この記事の環境で実際に出たものです。

エラー・警告 原因 直し方
ValueError: Expected 2D array, got 1D array instead: scikit-learn に1次元の x を渡した x.reshape(-1, 1) で「n 行 1 列」にする
ValueError: Cannot calculate a linear regression if all x values are identical linregress の x がすべて同じ x にばらつきのあるデータを使う
RankWarning: Polyfit may be poorly conditioned 点の数に対して次数が高い、x がすべて同じ、x の値が大きい 次数を下げる、Polynomial.fit を使う、x を平均0にそろえる
TypeError: expected x and y to have same length x と y の個数が違う 欠けた行を確かめてそろえる
ValueError: Input y contains NaN. scikit-learn に欠損値を渡した dropna() などで欠損を除くか埋める
np.polyfit が [nan nan] を返す y に NaN が入っている(エラーにならない) np.isfinite で NaN の行を除いてから計算
LinAlgError: Incompatible dimensions lstsq の X の行数と y の長さが違う X.shape と y.shape を確かめる
RuntimeError: Optimal parameters not found: Number of calls to function has reached maxfev = 5. curve_fit が上限回数までに収束しない 初期値 p0 を近づける、maxfev を増やす、x をそろえる
ValueWarning: omni_normtest is not valid with less than 8 observations statsmodels の summary で点が8未満 警告。点を増やすと正規性の検定ができる

scikit-learn の1次元のエラーは、続けて次の文が出ます。

Expected 2D array, got 1D array instead:
array=[1. 2. 3. 4. 5.].
Reshape your data either using array.reshape(-1, 1) if your data has a single feature or array.reshape(1, -1) if it contains a single sample.

よくある質問

Q. 最小二乗法とは簡単に言うと何ですか? A. データの点と式(直線など)の縦のずれ(残差)を二乗して足した値が、いちばん小さくなるように式の係数を決める方法です。5点の例では y = 0.9x + 1.3 になりました。

Q. 最小二乗法の公式は? A. 傾き a = Sxy ÷ Sxx(x と y の偏差の積の和 ÷ x の偏差の二乗和)、切片 b = ȳ − a × x̄ です。合計だけで書くと a = (n·Σxy − Σx·Σy) ÷ (n·Σx² − (Σx)²) です。

Q. 最小二乗法はなぜ二乗するのですか? A. 残差の符号をなくすためです。二乗なら微分して0とおくだけで公式が出て、答えが1つに決まり、誤差が正規分布なら最尤推定とも一致します。その代わり外れ値に弱くなります。

Q. 最小二乗法の公式はどうやって導きますか? A. 残差の二乗和 S(a, b) を a と b でそれぞれ偏微分して0とおくと、2本の連立方程式(正規方程式)が出ます。それを解くと傾きと切片の公式になります。「公式の導き方(偏微分)」の節に途中式があります。

Q. 最小二乗法と回帰分析の違いは? A. 最小二乗法は係数を決める計算の方法、回帰分析は変数の関係を式で調べる分析全体です。回帰分析の係数の推定に、ふつう最小二乗法を使います。

Q. Pythonで最小二乗法を計算するには? A. 直線なら a, b = np.polyfit(x, y, 1) です。scikit-learn の LinearRegression、SciPy の stats.linregress、statsmodels の OLS でも同じ答えになり、statsmodels は標準誤差と p 値も出します。

Q. 最小二乗法を行列で書くと? A. XᵀX β = Xᵀy(正規方程式)で、β = (XᵀX)⁻¹ Xᵀy です。プログラムで解くときは、XᵀX を作らない np.linalg.lstsq のほうが精度が落ちにくく、8次の多項式で誤差が 0.035 と 2.1e-08 の差になりました。

Q. 2次関数や曲線にも使えますか? A. 使えます。多項式は係数について1次式なので np.polyfit(x, y, 2) で解けます。指数関数などは scipy.optimize.curve_fit で、初期値から反復して解きます。

Q. 決定係数R²はどう求めますか? A. 1 − 残差の二乗和 ÷ y の全体の変動 です。5点では 1 − 1.90 ÷ 10 = 0.81 でした。直線(説明変数1つ)では相関係数の2乗と同じです。

Q. 最小二乗法の欠点は? A. 外れ値に弱いことです。範囲の外に1点足すと、傾きが 2.011 から 0.500 に変わりました。外れ値が心配なら、Huber 回帰・Theil-Sen・RANSAC などを使います。

Q. 原点を通る直線にするには? A. a = Σxy ÷ Σx² で求めます。scikit-learn では LinearRegression(fit_intercept=False) です。R² の定義がライブラリで違うので、比べるときは注意します。

Q. 最小二乗法は誰が考えましたか? A. Wikipedia によると、最初に発表したのはルジャンドル(1805年)です。ガウスは1809年に出版し、1795年には考えていたと主張しました。

まとめ

  • 最小二乗法は、残差(データと式の縦のずれ)の二乗和を最小にするように係数を決める方法。
  • 直線は 傾き a = Sxy ÷ Sxx、切片 b = ȳ − a·x̄。5点の例は y = 0.9x + 1.3、R² = 0.81。直線は平均の点を通り、残差の合計は0。
  • 公式は残差の二乗和を偏微分して0とおくと出る。行列では XᵀX β = Xᵀy(正規方程式)。プログラムでは lstsq などで解くと精度が落ちにくい。
  • Python は np.polyfit・LinearRegression・stats.linregress・sm.OLS、SQL は regr_slope・regr_intercept で同じ答えになった。
  • 外れ値に弱く、範囲の外の1点で傾きが 2.011 → 0.500。次数を上げすぎると過学習(9次で新しい点の誤差 646.563)。
  • 精度が違う測定は重み付き、式の形が決まった曲線は curve_fit(初期値に注意)、外れ値は頑健な回帰を使う。

参考・出典

確認日はいずれも 2026年10月22日です。

  • Wikipedia「最小二乗法」:https://ja.wikipedia.org/wiki/最小二乗法
  • Wikipedia「ガウス=マルコフの定理」:https://ja.wikipedia.org/wiki/ガウス=マルコフの定理
  • NumPy「numpy.polyfit」:https://numpy.org/doc/stable/reference/generated/numpy.polyfit.html
  • NumPy「numpy.polynomial.polynomial.Polynomial.fit」:https://numpy.org/doc/stable/reference/generated/numpy.polynomial.polynomial.Polynomial.fit.html
  • NumPy「numpy.linalg.lstsq」:https://numpy.org/doc/stable/reference/generated/numpy.linalg.lstsq.html
  • SciPy「scipy.stats.linregress」:https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.linregress.html
  • SciPy「scipy.optimize.curve_fit」:https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html
  • SciPy「scipy.optimize.least_squares」:https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html
  • scikit-learn「LinearRegression」:https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html
  • scikit-learn「Linear Models」(Ridge・Huber・Theil-Sen・RANSAC・Quantile):https://scikit-learn.org/stable/modules/linear_model.html
  • statsmodels「OLS」:https://www.statsmodels.org/stable/generated/statsmodels.regression.linear_model.OLS.html
  • statsmodels「WLS」:https://www.statsmodels.org/stable/generated/statsmodels.regression.linear_model.WLS.html
  • statsmodels「OLSInfluence」:https://www.statsmodels.org/stable/generated/statsmodels.stats.outliers_influence.OLSInfluence.html
  • PostgreSQL 16 Documentation「9.21. Aggregate Functions」:https://www.postgresql.org/docs/16/functions-aggregate.html

本記事の内容は、ご自由にお使いください。
ご利用の際は、出典として本ページへのリンクを記載いただけますようお願いします。

(記載例)出典:株式会社RJC「最小二乗法とは|公式・手計算の例とPythonでの求め方を実例で解説」

株式会社RJC ― SI事業・SES事業・AI駆動開発。RJCは一緒に成長を楽しめる会社です。

WE ARE HIRING

RJCで一緒に開発しながら、
成長を楽しみませんか?

RJCは、Web・モバイル・AIを活用した開発プロジェクトで、テックリードやPM・PMOも活躍するシステム開発会社です。会社を知る、待遇を確かめる、話を聞いてみる。気になるところ見てみてください!

  • 127日年間休日
  • 12時間平均残業時間
  • 毎日ガチャ遊びココロも大切にする福利厚生。アマギフなどの賞品ラインナップ!

ほかにも、チケットレストラン、書籍読み放題、2年ごとの慰労報奨(休暇 or 金一封)、11期連続の黒字決算。

ABOUT RJC RJCがどんな会社か知る 考え方、研修、働き方、福利厚生、社員の前職まで。RJCのことが丸わかり! RJC丸わかりページへ JOB DESCRIPTION 仕事内容・待遇を見てみる 仕事内容、給与・待遇、選考の流れ。経験者も未経験も!応募前に知りたいこと、まとめました! 募集要項を見る ENTRY エントリーする エントリーは1〜2分・履歴書不要です。まずは話を聞いてみたい、という方でも歓迎です! エントリーフォームへ

RJCで一緒に開発しながら、 成長を楽しみませんか?