はじめに

scipy.optimize.curve_fitの戻り値であるpcovは何かと思い検索すると、分散共分散行列であると出てくる。
名前から分散と共分散が入ってる行列なんだろうなと想像される。
だからと言ってそんな認識でいてはいけないと思い調べてみる。

データの準備

お決まりのirisを使う。

import seaborn as sns
import matplotlib.pyplot as plt

df = sns.load_dataset('iris')
sns.pairplot(df)
plt.show()

f:id:freedom-_-q:20210703214411p:plain

相関してそうなpetal_lengthとpetal_widthを使ってフィッティングする。

とりあえずどんな数字か見てみよう

from scipy.optimize import curve_fit

f = lambda x, a, b: a*x + b
popt, pcov = curve_fit(f, df['petal_length'], df['petal_width'])
print(pcov)
#[[ 9.18230745e-05 -3.45071113e-04]
# [-3.45071113e-04  1.58101581e-03]]
plt.imshow(pcov, cmap='jet')
plt.colorbar()
plt.show()

f:id:freedom-_-q:20210703215652p:plain

分散共分散行列の構造

分散共分散行列 $\Sigma$ は以下のようになっていて、対角成分に分散、その他に共分散が入っている。

$\Sigma =\begin{pmatrix} \sigma _{x}^{2} & \sigma _{xy} \\ \sigma _{xy} & \sigma _{y}^{2} \end{pmatrix}$

対角成分のルートを取れば標準偏差になるはず。

perr = np.sqrt(np.diag(pcov))
print(perr)
# array([9.18230745e-05, 1.58101581e-03])

...ふーん、て感じ。

決定係数を求める

フィッティングしたのだから、うまくできているか気になる。
というわけで指標となる決定係数を求める。

素直に求める方法

$R^{2}=1-\dfrac{\sum \left( y_{i}-f\left( x_{i}\right) \right) ^{2}}{\sum \left( y_{i}-\mu y\right) ^{2}}$

r2 = 1 - (np.sum((df['petal_width'] - f(df['petal_length'], *popt))**2) /
          np.sum((df['petal_width'] - np.mean(df['petal_width']))**2))
print(r2)
#0.9271098389904927

sklearn

実務ではこれ一択でしょ。

from sklearn.metrics import r2_score

r2 = r2_score(df['petal_width'], f(df['petal_length'], *popt))
print(r2)
#0.9271098389904927

決定係数を2乗する方法

最小二乗法による直線フィッティングの場合，相関係数の二乗と決定係数は一致する

だそうなので*1、決定係数 $\rho$ を求め2乗する。

$\rho =\dfrac{\sigma _{xy}}{\sigma _{x}\sigma _{y}}$

rho = pcov[0,1]/np.multiply(*perr)
print(rho**2)
#0.8202177539001069

なんか違う値出てきちゃうんだけどなんで？？

*1:https://manabitimes.jp/math/1016

freedom-_-qの勉強履歴

メモ書きが主になるかと思います。勉強強制のために一日一記事目指してます。頭良くないので間違いが多々あるかと思います。

分散共分散行列pcovの理解と可視化