【python】最小二乗法の信頼区間と予測区間を導出する方法

はじめに

過去の記事でpythonで最小二乗法を使って実験データの近似直線や近似曲線を求める方法について説明しました。

今回は最小二乗法で求めた近似値の信頼区間および予測区間を導出し、グラフに重ね合わせる手法について説明しようと思います。

はじめに
最小二乗法を使った1次関数の近似
信頼区間と予測区間
- 信頼区間の導出
- 予測区間の導出
コード実装
- 信頼区間をプロット
  - 信頼区間を塗りつぶし
- 予測区間をプロット
  - 信頼区間、予測区間を塗りつぶし
詳しく学びたい方は
全コード

最小二乗法を使った1次関数の近似

今回は1次関数への近似を題材として取り上げます。
過去の記事と同じデータを使ってまずは近似直線を描いてみましょう。

import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
sns.set()

# 測定データ
x = np.array([0, 5, 10, 15, 20, 25, 30])  # 時間
y = np.array([20, 24.2, 25.5, 30.5, 32.4, 36.5, 39.3]) # 温度

# 最小二乗近似
A = np.vstack((x, np.ones(len(x)))).T
v = y.T
u, residuals, rank, s = np.linalg.lstsq(A, v, rcond=None)
y_hat = x * u[0] + u[1]

# 測定データのプロット
fig, ax = plt.subplots()
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.plot(x, y, 'o', label='measure data')
ax.plot(x, y_hat, 'r', label='least square')
ax.legend(loc='best')
plt.show()

結果は以下です。

信頼区間と予測区間

信頼区間とは、母集団の真の値がどの範囲にあると推定できるかを示すことです。
つまりこのケースでいうと近似直線の取りうる範囲とも言えます。

予測区間とは、ある母集団の標本値がどの範囲にあると推定できるかを示すものです。
このケースでいうとサンプルデータが近似曲線から離れうる範囲のことです。

この二つは混同しやすいので注意が必要です。

信頼区間の導出

回帰直線は以下を考えます。

$y = ax + b$

ここで、サンプルデータと近似式の残差が正規分布に従っていることが前提になります。
その場合、信頼区間は下記の式で表すことができます。

$ax + b \pm t(N - 2, \alpha)\sqrt{\left(\dfrac{1}{N} + \dfrac{(x - \bar{x})^{2}}{\sum(x_{i} - \bar{x})^{2}} \right)V_{e}}$

ここで、 $N$ はサンプルサイズを示しています。

信頼区間は自由度 $N-2$ の $t$ 分布に従うことがわかっており、 $t(N - 2, \alpha)$ は自由度 $N-2$ の $t$ 分布における信頼度 $1- \alpha$ での $t$ 値を表しています。 $\alpha$ は優位水準を表しています。大体、信頼度は99% $(\alpha = 0.01)$ や95% $(\alpha = 0.05)$ が使われます。

また、 $V_{e}$ は誤差の分散であり、以下の式で求められます。

$V_{e} = \dfrac{\sum y_{i}^{2} - \sum \hat{y_{i}}^{2}}{(N - 2)}$

予測区間の導出

こちらもサンプルデータと近似式の残差が正規分布に従っていることが前提になりますが、下記式で表すことができます。
$ax + b \pm t(N - 2, \alpha)\sqrt{\left(1 + \dfrac{1}{N} + \dfrac{(x - \bar{x})^{2}}{\sum(x_{i} - \bar{x})^{2}} \right)V_{e}}$

コード実装

それでは上記式をもとに、信頼区間、予測区間それぞれ求めてグラフに重ねてみます。

信頼区間をプロット

まずは信頼区間を求めるためのt値及び $V_{e}$ , $\bar{x}$ の値を求めます。
優位水準 $\alpha$ は0.05とします(95%信頼区間)。

from scipy.stats import t
alpha = 0.05 # 優位水準
n = len(x)
tvalue = t.ppf(1 - alpha / 2, n - 2) # 自由度n-2の時のt値(alphaを1/2することに注意)
ve = (np.sum(y ** 2) - np.sum(y_hat ** 2)) / (n - 2)
x_bar = np.mean(x)

ここまでで準備は整いましたので、実際に信頼区間を求めてプロットしてみましょう。

# 信頼区間
conf_hi = y_hat + tvalue * np.sqrt(1 / n + (x - x_bar) ** 2 / np.sum((x - x_bar) ** 2)) * ve
conf_lo = y_hat - tvalue * np.sqrt(1 / n + (x - x_bar) ** 2 / np.sum((x - x_bar) ** 2)) * ve

# 測定データのプロット
fig, ax = plt.subplots()
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.plot(x, y, 'o', label='measure data')
ax.plot(x, y_hat, 'r', label='least square')
ax.plot(x, conf_hi, 'k:', label='conf_95%')
ax.plot(x, conf_lo, 'k:')
ax.legend(loc='best')
plt.show()

結果は下記の通りです。

信頼区間を塗りつぶし

ついでに、信頼区間を塗りつぶしてみましょう。
塗りつぶしにはax.fill_between()メソッドが便利です。

# 測定データのプロット
fig, ax = plt.subplots()
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.plot(x, y, 'o', label='measure data')
ax.plot(x, y_hat, 'r', label='least square')
ax.plot(x, conf_hi, 'k:', label='conf_95%')
ax.plot(x, conf_lo, 'k:')
ax.fill_between(x, conf_hi, conf_lo, facecolor='lightcoral', alpha=0.3)
ax.legend(loc='best')
plt.show()

結果は以下です。

95%信頼区間なので、95%の信頼度で、近似直線がこの塗りつぶしのエリア内に入りますという意味合いです。

予測区間をプロット

続いて予測区間をプロットしてみます。t値及び $V_{e}$ , $\bar{x}$ の値は信頼区間を求めた際に導出したので省きます。
優位水準 $\alpha$ は0.05とします(95%予測区間)。

# 予測区間
pred_hi = y_hat + tvalue * np.sqrt(1 + 1 / n + (x - x_bar) ** 2 / np.sum((x - x_bar) ** 2)) * ve
pred_lo = y_hat - tvalue * np.sqrt(1 + 1 / n + (x - x_bar) ** 2 / np.sum((x - x_bar) ** 2)) * ve

# 測定データのプロット
fig, ax = plt.subplots()
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.plot(x, y, 'o', label='measure data')
ax.plot(x, y_hat, 'r', label='least square')
ax.plot(x, pred_hi, 'k--', label='pred_95%')
ax.plot(x, pred_lo, 'k--')
ax.legend(loc='best')
plt.show()

結果は以下の通り。

予測区間なので、各サンプルデータが95%の信頼度でこの区間内に入るということを示しています。

信頼区間、予測区間を塗りつぶし

先ほどの信頼区間と合わせて表示してみましょう

# 測定データのプロット
fig, ax = plt.subplots()
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.plot(x, y, 'o', label='measure data')
ax.plot(x, y_hat, 'r', label='least square')
ax.plot(x, conf_hi, 'k:', label='conf_95%')
ax.plot(x, conf_lo, 'k:')
ax.fill_between(x, conf_hi, conf_lo, facecolor='lightcoral', alpha=0.5)
ax.plot(x, pred_hi, 'k--', label='pred_95%')
ax.plot(x, pred_lo, 'k--')
ax.fill_between(x, pred_hi, conf_hi, facecolor='lightcoral', alpha=0.2)
ax.fill_between(x, conf_lo, pred_lo, facecolor='lightcoral', alpha=0.2)
ax.legend(loc='best')
plt.show()

結果は以下の通り。

詳しく学びたい方は

この本がおすすめです。

リンク

PyQでも統計に関する多くのコースがありますので、おすすめです。

全コード

終わりに、全コードを載せておきます。

import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
from scipy.stats import t
sns.set()

# 測定データ
x = np.array([0, 5, 10, 15, 20, 25, 30])
y = np.array([20, 24.2, 25.5, 30.5, 32.4, 36.5, 39.3])

# 最小二乗近似
A = np.vstack((x, np.ones(len(x)))).T
v = y.T
u, residuals, rank, s = np.linalg.lstsq(A, v, rcond=None)
y_hat = x * u[0] + u[1]

alpha = 0.05 # 優位水準
n = len(x)
tvalue = t.ppf(1 - alpha / 2, n - 2) # t値(alphaを1/2することに注意)
ve = (np.sum(y ** 2) - np.sum(y_hat ** 2)) / (n - 2)
x_bar = np.mean(x)

# 信頼区間
conf_hi = y_hat + tvalue * np.sqrt(1 / n + (x - x_bar) ** 2 / np.sum((x - x_bar) ** 2)) * ve
conf_lo = y_hat - tvalue * np.sqrt(1 / n + (x - x_bar) ** 2 / np.sum((x - x_bar) ** 2)) * ve

# 予測区間
pred_hi = y_hat + tvalue * np.sqrt(1 + 1 / n + (x - x_bar) ** 2 / np.sum((x - x_bar) ** 2)) * ve
pred_lo = y_hat - tvalue * np.sqrt(1 + 1 / n + (x - x_bar) ** 2 / np.sum((x - x_bar) ** 2)) * ve

# 測定データのプロット
fig, ax = plt.subplots()
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.plot(x, y, 'o', label='measure data')
ax.plot(x, y_hat, 'r', label='least square')
ax.plot(x, conf_hi, 'k:', label='conf_95%')
ax.plot(x, conf_lo, 'k:')
ax.fill_between(x, conf_hi, conf_lo, facecolor='lightcoral', alpha=0.5)
ax.plot(x, pred_hi, 'k--', label='pred_95%')
ax.plot(x, pred_lo, 'k--')
ax.fill_between(x, pred_hi, conf_hi, facecolor='lightcoral', alpha=0.2)
ax.fill_between(x, conf_lo, pred_lo, facecolor='lightcoral', alpha=0.2)
ax.legend(loc='best')
plt.show()

morikomorou’s blog

自分が学んだことなどの備忘録的なやつ