データサイエンス講座(統計編) 10〜17 多変数の記述統計学（かめ@米国データサイエンティスト）のノート

2026-04-12statistics Python NumPy pandas Seaborn sklearn

「データサイエンス講座(統計編)」（かめ@米国データサイエンティストさん）のノートです。一変数の記述統計学がテーマだった前回の続きになります。第10回「共分散」から第17回「決定係数」まで、内容は 多変数の記述統計学 となります。

今回も前回に引き続き、多変数の統計量をスクラッチ実装していきます。多変数ということで、少し慣れが必要ですが記述の見通しを良くなるので、ベクトルと行列 を使っていきます。

さて今回の個人的ポイントは、NumPyの実装では アダマール積 と ブロードキャスト の組み合わせが有効だということです。ついつい慣れ親しんだ普通の行列積を考えがちですが、実装時には必ずともそれが最適とは限らないのです。→ 「11 相関係数をわかりやすく解説: 相関行列」

かめ@米国データサイエンティストさんの講座のお供に、今回も是非よろしくお願いします！JupyterLabのノートブックはこちらにあります。

10 絶対にわかる共分散

絶対にわかる共分散【統計学入門⑩】 - 米国データサイエンティストのブログ

分散共分散行列 (variance-covariance matrix)
- numpy.cov(m, y=None, bias=False, ddof=None): m と y は共分散を計算したいデータ（同じ型）、biasは不偏性の指定、更に ddof で bias の指定は上書きされる
- pandas.DataFrame.cov(ddof=1)

共分散(covariance)

s_{xy} := \frac{1}{n-1} \sum_i (x_i - \overline{x}) (y_i - \overline{y})

補足: ライブラリのデフォルトに合わせて $n - 1$ で割ることにする。

補足（重要！）: 変数 $x$ の偏差を並べたベクトル $\boldsymbol{x'} = (x_i - \overline{x})_i$ 、同様に変数 $y$ の偏差を並べたベクトル $\boldsymbol{y'}$ とした時、右側の $\sum_i (x_i - \overline{x}) (y_i - \overline{y})$ の部分は、ベクトルの内積 $\lang \boldsymbol{x'}, \boldsymbol{y'} \rang$ であることに注意しておこう。これが後の分散共分散行列や相関係数の理解に役立つ。

s_{xy} = \frac{1}{n-1} \lang \boldsymbol{x'}, \boldsymbol{y'} \rang

def my_cov(x: np.ndarray, y: np.ndarray) -> np.ndarray:
    return np.sum((x - np.mean(x)) * (y - np.mean(y))) / (len(x) - 1)

実行例

>>> weight = np.array([42, 46, 53, 56, 58, 61, 62, 63, 65, 67, 73])
>>> height = np.array([138, 150, 152, 163, 164, 167, 165, 182, 180, 180, 183])
>>> my_cov(weight, height)
np.float64(127.54545454545455)

分散共分散行列 (variance-covariance matrix)

$p$ 個の変数 $x_j$ ( $1 \le j \le p$ ) に対して、変数の（共）分散を成分に持つ $p$ 次正方行列 $S$

S := (s_{jk})_{j, k}\\

を 分散共分散行列 と呼ぶ。但し $s_{jk}$ は、変数 $x_j$ と変数 $x_k$ の（共）分散とする。 $j = k$ の場合、つまり対角成分は、変数 $x_j$ の分散 $s_j^2$ である。

ところで、標本サイズが $n$ の時、全てのデータを並べた行列を

X = (x_{ij})_{1 \le i \le n, 1 \le j \le p}

とする。但し $x_{ij}$ は、 $i$ 番目の個体の $j$ 番目の変数の値を表す。つまり pandas のデータフレームと同じ形式である。

この時、分散共分散行列 $S$ は

S = \frac{1}{n-1} {}^t\!X_C X_C \\ (X_C := (E_n - \frac{1}{n} J_n) X)

で計算できる。但し $E_n$ は $n$ 次単位行列、 $J_n$ は全成分が $1$ からなる $n$ 次正方行列である。また、 ${}^t\!X_C$ は転置行列。

weight と height の2変数の場合で具体的に見ていこう。X は各変数のデータを列ベクトルとして並べた行列。

>>> X = np.column_stack((weight, height))
>>> X
array([[ 42, 138],
       [ 46, 150],
       [ 53, 152],
       [ 56, 163],
       [ 58, 164],
       [ 61, 167],
       [ 62, 165],
       [ 63, 182],
       [ 65, 180],
       [ 67, 180],
       [ 73, 183]])

numpy.column_stack — NumPy Manual

$J_n$ の前に、まずは $n$ 次の行ベクトル j から始める。

>>> n = X.shape[0]
>>> j = np.ones((1, n))
>>> j
array([[1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]])

numpy.ones — NumPy Manual

$\frac{1}{n} j X$ は、各変数の平均を並べた行ベクトルになる。

>>> j @ X / n
array([[ 58.72727273, 165.81818182]])
>>> np.mean(weight), np.mean(height)
(np.float64(58.72727272727273), np.float64(165.8181818181818))

補足: @ 演算子は、行列同士の掛け算

Python, NumPyで行列の演算（逆行列、行列式、固有値など） | note.nkmk.me

j を $n$ 行並べた J

>>> J = np.ones((n, n))
>>> J
array([[1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]])

を掛けると、

>>> J @ X / n
array([[ 58.72727273, 165.81818182],
       [ 58.72727273, 165.81818182],
       [ 58.72727273, 165.81818182],
       [ 58.72727273, 165.81818182],
       [ 58.72727273, 165.81818182],
       [ 58.72727273, 165.81818182],
       [ 58.72727273, 165.81818182],
       [ 58.72727273, 165.81818182],
       [ 58.72727273, 165.81818182],
       [ 58.72727273, 165.81818182],
       [ 58.72727273, 165.81818182]])

各変数の平均を $n$ 行分並べた行列になる。つまり先ほどの $X$ から $\frac{1}{n} J_n X$ を引いた行列 $X_C$ は、 $X$ の成分から各変数の平均を引いた「偏差ベクトル」 $\boldsymbol{x_j'}$ を並べた行列になる。

>>> Xc = X - J @ X / n
>>> Xc
array([[-16.72727273, -27.81818182],
       [-12.72727273, -15.81818182],
       [ -5.72727273, -13.81818182],
       [ -2.72727273,  -2.81818182],
       [ -0.72727273,  -1.81818182],
       [  2.27272727,   1.18181818],
       [  3.27272727,  -0.81818182],
       [  4.27272727,  16.18181818],
       [  6.27272727,  14.18181818],
       [  8.27272727,  14.18181818],
       [ 14.27272727,  17.18181818]])

補足: NumPyではブロードキャストが行われるため、実際には行ベクトルjをかけても同じ結果が得られる。実装時には、こちらのほうが効率が良い。

>>> X - j @ X / n
array([[-16.72727273, -27.81818182],
       [-12.72727273, -15.81818182],
       [ -5.72727273, -13.81818182],
       [ -2.72727273,  -2.81818182],
       [ -0.72727273,  -1.81818182],
       [  2.27272727,   1.18181818],
       [  3.27272727,  -0.81818182],
       [  4.27272727,  16.18181818],
       [  6.27272727,  14.18181818],
       [  8.27272727,  14.18181818],
       [ 14.27272727,  17.18181818]])

従って ${}^t\!X_C X_C$ は偏差ベクトルの内積を並べた行列になる。

{}^t\!X_C X_C = (\lang \boldsymbol{x_j'}, \boldsymbol{x_k'} \rang )_{j, k}

先ほど共分散で述べた通り、偏差ベクトルの内積から $1-n$ を割ると共分散になるので、

\frac{1}{1-n} {}^t\!X_C X_C = (\frac{1}{1-n} \lang \boldsymbol{x_j'}, \boldsymbol{x_k'} \rang )_{j, k} = (s_{jk})_{j, k} = S

になることが分かる。

>>> Xc.T @ Xc / (m - 1)
array([[ 82.81818182, 127.54545455],
       [127.54545455, 218.76363636]])

numpy.ndarray.T — NumPy Manual

以上を関数にする。

def my_cov(*args: np.ndarray) -> np.ndarray:
    X = np.column_stack(args)
    n = X.shape[0]

    j = np.ones((1, n))
    # ブロードキャストされるので正方行列Jをかけた場合と同じ結果になる
    Xc = X - j @ X / n
    return Xc.T @ Xc / (n - 1)

複数個の変数を受け取れるように、可変長引数 *args を利用している。

Pythonの可変長引数（*args, **kwargs）の使い方 | note.nkmk.me

実行例

>>> my_cov(weight, height)
array([[ 82.81818182, 127.54545455],
       [127.54545455, 218.76363636]])

numpy.cov() の実行例1（2つのベクトルを渡す）

>>> np.cov(weight, height)
array([[ 82.81818182, 127.54545455],
       [127.54545455, 218.76363636]])

numpy.cov() の実行例2（データ行列を渡す）

>>> np.cov(X.T) # 転置
array([[ 82.81818182, 127.54545455],
       [127.54545455, 218.76363636]])

注意: numpy.cov() に渡すデータは、行に変数を、列にサンプルの個体を並べる形式を取る。従って、実行例2のように pandas と同様のデータ行列の場合、転置しなければならない。

ちなみに、引数 bias=True にすると、 $n - 1$ ではなく $n$ で割る。

>>> np.cov(weight, height, bias=True)
array([[ 75.2892562 , 115.95041322],
       [115.95041322, 198.87603306]])

更に ddof で bias の指定は上書きされる。

>>> np.cov(weight, height, ddof=0)
array([[ 75.2892562 , 115.95041322],
       [115.95041322, 198.87603306]])

numpy.cov — NumPy Manual

pandas.DataFrame.cov() の実行例

>>> df = pd.DataFrame({"weight": weight, "height": height})
>>> df
 	weight 	height
0 	42 	138
1 	46 	150
2 	53 	152
3 	56 	163
4 	58 	164
5 	61 	167
6 	62 	165
7 	63 	182
8 	65 	180
9 	67 	180
10 	73 	183
>>> df.cov()
    	weight 	height
weight 	82.818182 	127.545455
height 	127.545455 	218.763636

pandas.DataFrame.cov — pandas documentation

11 相関係数をわかりやすく解説

相関係数をわかりやすく解説【統計学入門11】 - 米国データサイエンティストのブログ

相関行列 (correlation matrix)
- numpy.corrcoef(x, y=None)
- pandas.DataFrame.corr()
行列操作
- 行列の対角成分、対角行列: numpy.diag(v, k=0)
- ndarrayに次元を追加: ndarray[:, np.newaxis]

相関係数 (correlation coefficient)

先ほど共分散で述べたように、共分散は偏差ベクトルの内積を $n-1$ で割ったものであった。同様に、標準偏差 $s_x$ は偏差ベクトル $\boldsymbol{x'}$ の ノルム $\| \boldsymbol{x'} \|$ を割ったものになる。

s_{x} = \sqrt{ \frac{1}{n-1} \lang \boldsymbol{x'}, \boldsymbol{x'} \rang} = \frac{1}{\sqrt{n-1}} \| \boldsymbol{x'} \| \\

従って、コーシー・シュワルツの不等式 より、不等式

- s_x s_y \le s_{xy} \le s_x s_y

が成り立つ。そこで 相関係数 を

r_{xy} := \frac{s_{xy}}{s_x s_y} (= \frac{\lang \boldsymbol{x'}, \boldsymbol{y'} \rang}{\| \boldsymbol{x'} \| \| \boldsymbol{y'} \|} )

と定義すると、相関係数の範囲は

-1 \le r_{xy} \le 1

となる。相関係数の最小と最大に関して次のことが言える。

$r=1 \iff$ 偏差ベクトル $\boldsymbol{x'}$ と $\boldsymbol{y'}$ の向きが同じ $\iff$ 変数 $x$ と変数 $y$ が 傾き正 で一直線上に並ぶ
$r=-1 \iff$ 偏差ベクトル $\boldsymbol{x'}$ と $\boldsymbol{y'}$ の向きが逆 $\iff$ 変数 $x$ と変数 $y$ が 傾き負 で一直線上に並ぶ

def my_corrcoef(x: np.ndarray, y: np.ndarray) -> np.ndarray:
	cov_mat = np.cov(x, y)
	s_x_y = cov_mat[1, 0]
	s_x, s_y = np.sqrt(np.diag(cov_mat))
    return s_x_y / (s_x * s_y)

4行目では共分散行列から対角成分のベクトルを取り出して、平方根を計算している。numpy.diag() 関数は、行列から対角成分の抽出にも、対角行列の生成にも使える。

numpy.diag — NumPy v2.4 Manual

実行例

>>> weight = np.array([42, 46, 53, 56, 58, 61, 62, 63, 65, 67, 73])
>>> height = np.array([138, 150, 152, 163, 164, 167, 165, 182, 180, 180, 183])
>>> my_corrcoef(weight, height)
np.float64(0.9475771359123567)

相関行列 (correlation matrix)

分散共分散行列と同様に、 $p$ 個の変数 $x_j$ ( $1 \le j \le p$ ) に対して、変数の相関係数を成分に持つ $p$ 次正方行列 $R$

R := (r_{jk})_{j, k}\\

を 相関行列 と呼ぶ。但し $r_{jk}$ は、変数 $x_j$ と変数 $x_k$ の相関係数とする。対角成分は全て $1$ になる。

さて、分散共分散行列 $S$ の対角成分 $s_{jj}$ （変数 $x_j$ の分散 $s_j^2$ ）からなる $p$ 次対角行列を

D := diag(S) = diag(s_{jj})_j

とする。更に、 $D$ の対角成分の平方根の逆数（つまり標準偏差 $s_j$ の逆数）を成分にもつ対角行列を

D^{-1/2} := diag(d_j^{-1/2})_j = diag(1/ \sqrt{s_{jj}})_j = diag(1/ s_j)_j

としたとき、相関行列は

C = D^{-1/2} S D^{-1/2}

で計算できる。対角行列を 左から かけると、行列 $S$ の対応する 行がスカラー倍 される。一方、対角行列を 右から かけると、行列 $S$ の対応する 列がスカラー倍 される。

def my_corrcoef(*args: np.ndarray) -> np.ndarray:
    X = np.column_stack(args)
    S = np.cov(X.T)

    D_inv_std = np.diag(np.diag(S) ** (-1 / 2))
    return D_inv_std @ S @ D_inv_std

実行例

>>> my_corrcoef(weight, height)
array([[1.        , 0.94757714],
       [0.94757714, 1.        ]])

但し、この実装は np.diag() 関数を二重に使って、ゼロ埋め行列を生成するので効率が悪い。

アダマール積 * を利用した、別の実装を考える。 $D^{-1/2}$ の対角成分からなるベクトル inv_std を計算する。

>>> inv_std = np.diag(S) ** (-1 / 2)
>>> inv_std
array([0.10988468, 0.06761023])

inv_std を S にアダマール積でかけると、

>>> np.broadcast_to(inv_std, S.shape)
array([[0.10988468, 0.06761023],
       [0.10988468, 0.06761023]])

で確認できるように、自動的にブロードキャストされて、S の各列へのスカラー倍になる。

次に inv_std を転置して、列ベクトルに変換する。

>>> inv_std[:, np.newaxis]
array([[0.10988468],
       [0.06761023]])

NumPy配列ndarrayに次元を追加するnp.newaxis, np.expand_dims() | note.nkmk.me

これを S にアダマール積でかけると、

>>> np.broadcast_to(inv_std[:, np.newaxis], S.shape)
array([[0.10988468, 0.10988468],
       [0.06761023, 0.06761023]])

で確認できるように、自動的にブロードキャストされて、S の行へのスカラー倍になる。従ってこれらアダマール積でかけ合わせると、

>>> inv_std[:, np.newaxis] * S * inv_std
array([[1.        , 0.94757714],
       [0.94757714, 1.        ]])

相関行列が得られる。アダマール積を利用することで、効率が改善している。

以上を関数にまとめると、次のようになる。

def my_corrcoef(*args: np.ndarray) -> np.ndarray:
    X = np.column_stack(args)
    S = np.cov(X.T)

    inv_std = np.diag(S) ** (-1 / 2)
    return inv_std[:, np.newaxis] * S * inv_std

numpy.corrcoef() の実行例

>>> np.corrcoef(weight, height)
array([[1.        , 0.94757714],
       [0.94757714, 1.        ]])

numpy.corrcoef — NumPy v2.4 Manual

pandas.DataFrame.corr() の実行例

>>> df = pd.DataFrame({"weight": weight, "height": height})
>>> df
 	weight 	height
0 	42 	138
1 	46 	150
2 	53 	152
3 	56 	163
4 	58 	164
5 	61 	167
6 	62 	165
7 	63 	182
8 	65 	180
9 	67 	180
10 	73 	183
>>> df.corr()
    	weight 	height
weight 	1.000000 	0.947577
height 	0.947577 	1.000000

pandas.DataFrame.corr — pandas 3.0.2 documentation

12 これだけは知っておいた方がいい相関係数のポイント3つ

これだけは知っておいた方がいい相関係数のポイント3つ【統計学入門12】 - 米国データサイエンティストのブログ

実験: 相関係数の強弱を確かめる

相関係数の強弱を確かめるために、実際に乱数を生成してプロットする実験を行う。

その前に、講座で与えられた確率変数 $X, Y$ が、所望の相関係数 $r$ になることを確かめる。まず確率変数 $A, E1, E2$ が互いに独立に標準正規分布に従うと仮定する。それぞれの共分散は0である。

$r \ge 0$ のとき、確率変数を

\begin{aligned} X &:= \sqrt{r} A + \sqrt{1-r} E_1 \\ Y &:= \sqrt{r} A + \sqrt{1-r} E_2 \end{aligned}

と定めると、分散は

V(X) = rV(A) + 2\sqrt{r} \sqrt{1-r} Cov(A, E_1) + (1-r)V(E_1) = 1 \\ V(Y) = rV(A) + 2\sqrt{r} \sqrt{1-r} Cov(A, E_2) + (1-r)V(E_2) = 1

であり、共分散は

\begin{aligned} Cov(X, Y) &= Cov(\sqrt{r} A + \sqrt{1-r} E_1, \sqrt{r} A + \sqrt{1-r} E_2) \\ &= rV(A) + \sqrt{r} \sqrt{1-r} Cov(A, E_2) + \sqrt{1-r} \sqrt{r} Cov(E_1, A) + (1-r) Cov(E_1, E_2) \\ &= r \end{aligned}

となるので、 $X, Y$ の相関係数は $r$ となる。

一方、 $r \lt 0$ のときは、

\begin{aligned} X &:= -\sqrt{-r} A + \sqrt{1+r} E_1 \\ Y &:= \sqrt{-r} A + \sqrt{1+r} E_2 \end{aligned}

と、正の平方根になるように定める。

V(X) = -rV(A) - 2\sqrt{-r} \sqrt{1+r} Cov(A, E_1) + (1+r)V(E_1) = 1 \\ V(Y) = -rV(A) + 2\sqrt{-r} \sqrt{1+r} Cov(A, E_2) + (1+r)V(E_2) = 1 \\

\begin{aligned} Cov(X, Y) &= Cov(-\sqrt{-r} A + \sqrt{1+r} E_1, \sqrt{-r} A + \sqrt{1+r} E_2) \\ &= rV(A) - \sqrt{-r} \sqrt{1+r} Cov(A, E_2) + \sqrt{1+r} \sqrt{-r} Cov(E_1, A) + (1+r) Cov(E_1, E_2) \\ &= r \end{aligned}

相関係数は $r$ になる。

相関係数を変化させながら一覧でプロットしてみると、まさに百聞は一見にしかず。

相関の大小の基準

0.7を超えたら相関は強い
0.2を下回るようなら相関はほぼない

13 回帰分析を図でわかりやすく解説！条件付き平均と最小二乗法って？

回帰分析を図でわかりやすく解説！条件付き平均と最小二乗法って？【統計学入門13】 - 米国データサイエンティストのブログ

散布図と回帰直線: seaborn.regplot(data=None, *, x=None, y=None)
- seaborn.regplot(data=df, x="column name", y="column name"): DataFrameの場合
- seaborn.regplot(x=array_x, x=array_y): ndarrayの場合

補足（Seabornの関数の呼び出し方法）: v0.11.0 で位置引数が非推奨になったことに伴い、v0.12 以降は記事のように sns.regplot(height, weight) と位置引数での呼び出しができなくなった。data が唯一の位置引数になり、他は キーワード引数で指定 する必要がある。

Release v0.12.0 · mwaskom/seaborn

呼び出し方法1: DataFrameを渡す（xとyは列名） sns.regplot(data=df, x="height", y="weight")

呼び出し方法2: ndarray を渡す sns.regplot(x=height, y=weight)

seaborn.regplot — seaborn 0.13.2 documentation

14 回帰直線の数式をイメージで理解する．回帰係数と相関係数の関係

回帰直線の数式をイメージで理解する．回帰係数と相関係数の関係【統計学入門14】 - 米国データサイエンティストのブログ

線形回帰: sklearn.linear_model.LinearRegression()

最小2乗法

残差平方和を最小にする $a, b$ を求める。

\argmin_{a, b} \sum_i (y_i - (a + bx_i))^2

$a, b$ に関して偏微分すると、正規方程式

\begin{aligned} \sum_i y_i &= na + b \sum x_i \\ \sum_i x_i y_i &= a \sum_i x_i + b \sum x_i^2 \end{aligned}

が得られる。この連立方程式を解くと、

\begin{aligned} a &= \overline{y} - b \overline{x} \\ b &= \frac{s_{xy}}{s_x^2} = r \frac{s_{y}}{s_x} \end{aligned}

回帰係数が得られる。回帰係数 $b$ に関しては、講座にあるように $r s_{y} / s_x$ が直感的で理解しやすいが、実際の計算は $s_{xy} / s_x^2$ の方が単純。

回帰係数の計算法1 ( $s_{xy} / s_x^2$ )

>>> weight = np.array([42, 46, 53, 56, 58, 61, 62, 63, 65, 67, 73])
>>> height = np.array([138, 150, 152, 163, 164, 167, 165, 182, 180, 180, 183])
>>> S = np.cov(height, weight)
>>> b = S[0, 1] / S[0, 0]
>>> b
np.float64(0.583028590425532)

回帰係数の計算法2 ( $r s_{y} / s_x$ )

>>> r = np.corrcoef(height, weight)[0, 1]
>>> s_x = np.std(height)
>>> s_y = np.std(weight)
>>> b = r * s_y / s_x
>>> b
np.float64(0.5830285904255318)

計算法1 ( $s_{xy} / s_x^2$ ) で関数にまとめると以下のようになる。

def my_linear_regression(x: np.ndarray, y: np.ndarray):
    S = np.cov(x, y)
    b = S[0, 1] / S[0, 0]
    a = np.mean(y) - b * np.mean(x)
    return (a, b)

実行例

>>> my_linear_regression(height, weight)
(np.float64(-37.949468085106396), np.float64(0.583028590425532))

sklearn.linear_model.LinearRegression() の実行例

>>> from sklearn.linear_model import LinearRegression
>>> # 説明変数: 軸を追加して行列（列ベクトル）にする
>>> X = height[:, np.newaxis]
>>> # 目的変数
>>> y = weight
>>> # インスタンス作成
>>> reg = LinearRegression()
>>> # 学習
>>> reg.fit(X, y)
>>> # 学習結果
>>> print(f"b={reg.coef_}")
b=[0.58302859]
>>> print(f"a={reg.intercept_}")
a=-37.94946808510644
>>> # 予測
>>> X = np.array([[175]])
>>> y = reg.predict(X)
>>> print(X, y)
[[175]] [64.08053524]

fit() メソッド: 学習
predict() メソッド: 予測

LinearRegression — scikit-learn 1.8.0 documentation

16 決定係数R2を超わかりやすく解説

決定係数R2を超わかりやすく解説【統計学入門16】 - 米国データサイエンティストのブログ

決定係数のイメージに関しては、かめさんが講座で十分説明してくれているので、ここでは数学的な補足を少しだけ。

（補足）平方和の分解と決定係数

データ $y_i$ の偏差は

y_i - \overline{y} = (\hat{y}_i - \overline{y}) + (y_i - \hat{y}_i)

予測値 $\hat{y}_i$ の偏差（右辺左側）と残差（右辺右側）に分けられる。最小2乗法では、予測値 $\hat{y}_i$ は説明変数 $x_i$ が生成する空間への直交射影になるので、右辺は左辺の 直交分解 になっている（詳しいことは今回は省略）。従って、次式のように平方和の分解が成り立つ。

\sum_i (y_i - \overline{y})^2 = \sum_i (y_i - \hat{y}_i)^2 + \sum_i (\hat{y}_i - \overline{y})^2

27-5. 決定係数と重相関係数 | 統計学の時間 | 統計WEB

左辺（「全変動」と呼ぶ）に対して、右辺左側（「回帰変動」と呼ぶ）は説明変数で 説明される変動、右辺右側（「残差変動」と呼ぶ）は説明変数で 説明されない変動 と考えられる。

そこで、全変動のうち 説明変数で説明される割合 を 決定係数 $R^2$ と定義する。

R^2 := 1 - \frac{\sum_i (y_i - \hat{y}_i)^2}{\sum_i (y_i - \overline{y})^2} ( = \frac{\sum_i (\hat{y}_i - \overline{y})^2}{\sum_i (y_i - \overline{y})^2} )

（補足）決定係数と相関係数の関係

最小2乗法では $\hat{y}_i - \overline{y} = b(\hat{x}_i - \overline{x})$ なので、

R^2 = \frac{\sum_i (\hat{y}_i - \overline{y})^2}{\sum_i (y_i - \overline{y})^2} = \frac{b^2 \sum_i (x_i - \overline{x})^2}{\sum_i (y_i - \overline{y})^2}

と分子を変形する。更に $b = r s_y / s_x$ を代入すると、

\frac{b^2 \sum_i (x_i - \overline{x})^2}{\sum_i (y_i - \overline{y})^2} = r^2 \frac{s_y^2}{s_x^2} . \frac{s_x^2}{s_y^2} = r^2

つまり、最小2乗法では決定係数は相関係数の平方になる。

17 決定係数R2で回帰式の精度を評価する

決定係数R2で回帰式の精度を評価する【統計学入門17】 - 米国データサイエンティストのブログ

決定係数: sklearn.metrics.r2_score(y_true, y_pred)

先ほど定義した決定係数 $R^2$ を実装する。引数 y_true はデータ y を、y_pred は予測値 $\hat{y}$ を表す。

def my_r2_score(y_true: np.ndarray, y_pred: np.ndarray):
    S_t = np.sum((y_true - np.mean(y_true)) ** 2)
    S_e = np.sum((y_true - y_pred) ** 2)
    return 1 - S_e / S_t

講座では学習データと評価データに分けて、評価データで決定係数を計算している。

実行例

>>> # 評価データ
>>> weight_test = np.array([43, 45, 50, 58, 58, 60, 66, 63, 65, 70, 72])
>>> height_test = np.array([140, 148, 152, 163, 160, 170, 163, 177, 177, 185, 180])
>>> # 予測値
>>> X_test = height[:, np.newaxis]
>>> y_test_pred = reg.predict(X_test)
>>> # 決定係数
>>> my_r2_score(weight_test, y_test_pred)
np.float64(0.8569626575763964)

sklearn.metrics.r2_score() の実行例

>>> from sklearn.metrics import r2_score
>>> r2_score(weight_test, y_test_pred)
0.8569626575763964

r2_score — scikit-learn 1.8.0 documentation