「入門はじめての多変量解析」をPythonで写経 Vol.5 ～ ２章「はじめての重回帰分析」④重回帰分析の統計的推定
２章「はじめての重回帰分析」
書籍の著者 石村貞夫 先生、石村光資郎 先生
書籍「入門はじめての多変量解析」２章「はじめての重回帰分析」の Python写経活動記録 です。
多変量解析の入門を Python と一緒に学ぶ写経シリーズです。
この記事は、重回帰分析の統計的推定 を取り扱います。
具体的には 予測値の区間推定、偏回帰係数の区間推定 に取り組みます。
ChatGPT 活用型学習で進めてまいります！
では書籍を開いて多変量解析の旅に出かけましょう🚀
はじめに
このブログシリーズは、書籍「入門はじめての多変量解析」（東京図書、「テキスト」と呼びます）の Python 写経を通じて得た「多変量解析の楽しさ」をご紹介します。
書籍の紹介と引用表記はリンク先の記事に掲載しています。
２章 はじめての重回帰分析
この記事は２章の以下のSectionを取り扱います。
2.7 重回帰分析の推定
記事に用いるデータは、テキストに掲載されたデータそのものを引用しています。
データ件数の少ないものはコード上でデータを登録し、データ件数の多いものはCSVファイル化してデータを読み込みしています。
この記事で用いるライブラリをインポートします。
### インポート
# 数値計算
import numpy as np
import pandas as pd
# 統計
import scipy.stats as stats
import statsmodels.formula.api as smf
from statsmodels.stats.outliers_influence import summary_table # 予測値の95%CI
# ユーティリティ
import io
# 描画
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo' # または import japanize_matplotlib
分析の準備
テキスト p.35 表 2.1.1 「温度・圧力・配向度のデータ」を引用いたします。
■ データの読み込み
pandas データフレームにデータを設定します。
### 温度・圧力・配向度のデータ p.35 表2.1.1
# データの登録
data1 = pd.DataFrame(
{'配向度': [45, 38, 41, 34, 59, 47, 35, 43, 54, 52],
'温度': [17.5, 17.0, 18.5, 16.0, 19.0, 19.5, 16.0, 18.0, 19.0, 19.5],
'圧力': [30, 25, 20, 30, 45, 35, 25, 35, 35, 40],
}, index=range(1, 11))
data1.index.name = 'サンプルNo.'
# 結果の表示
data1【実行結果】
データの個数（標本サイズ）は 10 です。
目的変数は「配向度」、説明変数の候補は「温度」「圧力」です。
■ 共通設定
説明変数名と目的変数名を設定します。
データフレームの列を特定する際に活躍します。
## 共通設定
# 説明変数
VARS = ['温度', '圧力']
# 目的変数
TARGET = '配向度'【実行結果】なし
■ Python ライブラリで重回帰分析結果を取得
ライブラリ statsmodels で、説明変数に「温度」「圧力」を用いて、重回帰分析を実施しておきます。
🖲️statsmodels
## 重回帰分析の実行 statsmodels利用
result1_sm = smf.ols(formula='配向度 ~ 温度 + 圧力', data=data1).fit()
result1_sm.summary()【実行結果】
重回帰分析の結果を格納した result1_sm は今回記事で頻出します！
重回帰分析のサマリーは以下のようになりました。
この記事は主に後段の「Intercept、温度、圧力」あたりを利用します。
統計的推定 p.68
統計的推定の概要について、総務省統計局「なるほど統計学園」に掲載の文章をお借りします。
推定は、母集団の特性値（平均や分散など）を標本のデータから統計学的に推測することで、推定には点推定と区間推定があります。
点推定で推定するのは1つの値で、区間推定ではある区間（幅）をもって値を推定します。
テキストは重回帰分析の統計的推定のうち、「予測値の区間推定」を紹介しています。
また、統計ソフトや Python のライブラリの重回帰分析結果に表示される「偏回帰係数の区間推定」もメジャーな推定でしょう。
予測値の区間推定 p.68～
目的変数の予測値 $${Y}$$ の区間推定です。
$${Y = b_0 + b_1 x_2 + b_2 x_2 + \cdots + b_p x_p}$$ で予測する「１点」の両側に「信頼区間」と呼ばれる「幅」を付与する感じです。
予測値の信頼区間を「予測区間」と呼ぶ場合もあります。
テキスト p.68 の公式をお借りします。
📊 予測値の区間推定の公式 p.68
説明変数 $${(x_1, x_2, \cdots, x_p)}$$ がある値 $${(x_{10}, x_{20}, \cdots, x_{p0})}$$ をとるとき、予測値の $${100(1-\alpha)\%}$$ 信頼区間は
$$
y_0 - t_{N-p-1} \left(\cfrac{\alpha}{2}\right) \sqrt{\left(1 + \cfrac{1}{N} + \cfrac{D_0^2}{N - 1}\right) V_E} \\
\\
\leq 予測値 \leq \\
\\
y_0 + t_{N-p-1} \left(\cfrac{\alpha}{2}\right) \sqrt{\left(1 + \cfrac{1}{N} + \cfrac{D_0^2}{N - 1}\right) V_E}
$$
【変数の意味】
$${y_0}$$：説明変数 $${(x_{10}, x_{20}, \cdots, x_{p0})}$$ に対応する予測値 $${Y}$$（点推定値）
$${D_0^2 = \sum_{j=1}^p \sum_{i=1}^p (x_{j0} - \bar{x}_j)(x_{i0} - \bar{x}_i) s^{ji}}$$
$${s^{ji}}$$：分散共分散行列の逆行列の $${j,i}$$ 成分
$${V_E}$$：誤差変動の不偏分散 $${S_E/(N-p-1)}$$ ※ $${S_E}$$ は残差平方和
🔢 予測値の信頼区間の計算
公式をPython関数化して、配向度データの予測値の信頼区間を計算しましょう。
予測値の信頼区間算出関数の定義から。
### 予測値の区間推定の関数定義 p.68
def pred_conf_interval(x_pred, X, result, alpha=0.05):
## 予測値y0の算出
y0 = result.predict(x_pred).values[0]
## t分布の上側α/2％点の算出
# 標本サイズn、説明変数の数pの取得
n, p = int(result.nobs), int(result.df_model)
# 自由度n-p-1のt分布の上側α/2％点の取得
t = stats.t.isf(q=alpha/2, df=n-p-1)
## s^jiの算出
# 説明変数間の分散共分散行列の算出
cov_mtx = X.cov(ddof=1).values
# s^jiの算出
s = np.linalg.inv(cov_mtx)
## D0^2の算出
# 説明変数の平均値の算出
x_bar = X.mean(axis=0).values
# D0^2の算出
D02 = sum(
[(x_pred.values[j] - x_bar[j]) * (x_pred.values[i] - x_bar[i]) * s[j, i]
for j in range(p) for i in range(p)])
## VEの算出
VE = result.mse_resid
## 信頼区間の掛け算部分の算出
multiplier = t * ((1 + 1/n + D02/(n-1)) * VE)**(1/2)
## 戻り値: 信頼区間
return y0 - multiplier, y0 + multiplierテストを兼ねて、テキスト p.68～69 のデータ No.1「温度=17.5、圧力=30」 の予測値の 95% 信頼区間を計算しましょう。
引数は次の３つです。
・x_pred = 予測に使う説明変数の値
・X = 重回帰分析で用いた説明変数
・result = statsmodels の 回帰分析結果
・alpha = 100 (1 - α) % の α の値。デフォルト 0.05（95% 信頼区間）
### 信頼区間の計算 p.68～69
# 設定
idx = 1 # 予測するデータのindex
# 予測に使用する説明変数データの取得
x_vals = data1.loc[idx, VARS]
# 信頼区間の計算 (95%信頼区間下端, 95%信頼区間上端)
pred_conf_interval(x_vals, data1[VARS], result1_sm)【実行結果】
$${33.15 \leq 予測値 \leq 50.85}$$ になりました。
説明変数が１つの単回帰分析の場合、信頼区間の可視化は簡単です！
横軸に説明変数、縦軸に目的変数をとり、縦方向に信頼区間の帯を描画します。
しかし重回帰分析の場合、説明変数が２つ以上あるので、平らなグラフ（ｘｙ平面）で表現できません。
そこで。。。
🛸ちょっと寄り道：実測値・予測値プロット 🛸
「目的変数だけ」に絞った可視化で信頼区間を活用しましょう！
実測値・予測値プロットは、横軸に目的変数の実測値、縦軸に目的変数の予測値をとって、散布図と４５度線を描画するチャートです。
「実測値＝予測値」を示す４５度線を目印にして、予測の適合具合い・外れ具合いを確認できます。
予測値の点に加えてエラーバーで信頼区間を示すことで、「予測値の信頼区間の幅を考慮した」予測の状況を確認できるようになります。
描画用の共通関数を作成しておきます。
# 実測値・予測値プロット（予測値の95%信頼区間エラーバー付き）の描画関数
def plot_conf_int_scatter(y_obs, y_pred, ci_lower, ci_upper):
# 描画領域の設定
plt.figure(figsize=(6, 4))
# 95%信頼区間のerrorbar付き散布図の描画
plt.errorbar(
y_obs,
y_pred,
yerr=[y_pred - ci_lower, ci_upper - y_pred],
fmt='o',
ms=8,
color='tab:blue',
ecolor='tab:blue',
elinewidth=1, # エラーバー線の太さ
capsize=3,
# alpha=0.7,
label='予測値 ± 95% CI'
)
# 観測値＝予測値の４５度線の描画
min_val = min(y_obs.min(), y_pred.min())
max_val = max(y_obs.max(), y_pred.max())
plt.plot([min_val, max_val], [min_val, max_val], color='tab:red', ls='--',
label='実測値 = 予測値')
# 修飾
plt.xlabel('実測値', fontsize=12)
plt.ylabel('予測値', fontsize=12)
plt.title('実測値 vs 予測値 ＋ 信頼区間')
plt.legend()
plt.show()さきほど作成した予測値の信頼区間算出関数を使って、プロットしましょう。
### 予測値と信頼区間の可視化
## data1の信頼区間算出
# 設定
df = data1[VARS] # 分析するデータフレーム
N = len(df) # 標本サイズ
# 信頼区間の算出
conf_intervals = np.zeros((N, 2))
for i in range(N):
x_vals = df.loc[i+1, VARS]
conf_intervals[i, :] = pred_conf_interval(x_vals, df, result1_sm)
## 描画関数引数用の設定
# 実測値
y_obs = data1[TARGET]
# 予測値
y_pred = result1_sm.fittedvalues
# 95%信頼区間の下端・上端
ci_lower, ci_upper = conf_intervals.T
## 描画
plot_conf_int_scatter(y_obs, y_pred, ci_lower, ci_upper)【実行結果】
すべての予測値の信頼区間は「実測値＝予測値」の４５度線（赤点線）を含んでいます。
実は statsmodels には回帰分析結果を用いて予測値の信頼区間を算出できる summary_table() があります。
summary_table() で予測値の信頼区間を取得して、プロットしましょう。
### 予測値と信頼区間の可視化 statsmodels利用
## 予測値と信頼区間を算出
# statsmodelsのツールで予測値と信頼区間を算出
table, vals, cols = summary_table(result1_sm)
# インデックス,実測値,予測値を取得
index, y_obs, y_pred = vals.T[:3]
# 予測値の信頼区間を取得
ci_lower, ci_upper = vals.T[6:8]
## 描画
plot_conf_int_scatter(y_obs, y_pred, ci_lower, ci_upper)【実行結果】
statsmodels を使えば、回帰分析で用いるさまざまな統計量を簡単に取得することができます！
🛸ちょっと寄り道：statsmodels の summary_table() を活用！ 🛸
ChatGPTに相談して調べたことをまとめます！
statsmodels 「予測値の区間推定」のデータは
statsmodels.stats.outliers_influence.summary_table
を用いて取得しています。
summary_table() と呼びましょう。
summary_table() で得られるデータを概観してみましょう。
# statsmodelの予測結果
df = pd.read_html(
io.StringIO(table.as_html()), header=[0, 1], index_col=0
)[0]
df【実行結果】
予測値の 95% 信頼区間を描画する際には「Predict ci」の「95% low」「95% upp」を使っています。
ChatGPTに各列の内容を整理してもらいました。
予測値に「平均予測値」と「個別予測値（新しい観測値予測）」の２つの概念が出てきました。
交通整理しましょう。
【 平均予測値と個別予測値の深堀り 】
ここまで用いてきた予測値の信頼区間は「個別予測値」の信頼区間 Predict ci です。
もう一方の「平均予測値」は重回帰式で計算される予測値のことであり、平均予測値の信頼区間は Mean ci です。
重回帰式が誤差 $${\varepsilon_i}$$ を含まないことに留意します。
２つの信頼区間の違いは…
Mean ci：重回帰式そのものの不確かさ（平均を予測する精度）
⇒ 「直線の高さ」のぶれ幅を考慮Predict ci：重回帰式＋個々のばらつき（新規観測の値そのものの範囲）
⇒ 「直線の高さのぶれ」と「データそのもののばらつき（誤差 $${\varepsilon_i}$$）」の両方を考慮した幅
ちなみに…
・Predicted Value は平均予測値そのもの
・Std Error Mean Predict は平均予測値の標準誤差
です。
➡️ 平均予測値に関連する項目
・平均予測値、標準誤差、平均予測値の 95% 信頼区間
➡️ 個別予測値に関連する項目
・個別予測値の 95% 信頼区間
これらを可視化して比べてみましょう。
ざっくり比較なので、チャートの適否には目を瞑っておきます。
### 予測値・標準誤差・平均予測値の95%信頼区間・個別予測値の95%信頼区間の違いを可視化
# 描画領域の設定
plt.figure(figsize=(8, 5))
# 実測値の散布図（青い点）の描画
plt.scatter(df.index, df.iloc[:, 0], label='実測値', zorder=10)
# 平均予測値の線の描画
plt.plot(df.index, df.iloc[:, 1], color='tab:red', lw=2,
label='重回帰式による予測値', zorder=8)
# 平均予測値の標準誤差の塗りつぶし
plt.fill_between(df.index,
df.iloc[:, 1] + df.iloc[:, 2], df.iloc[:, 1] - df.iloc[:, 2],
color='lightpink',label='平均予測値の標準誤差', zorder=6)
# 平均予測値の95%信頼区間の塗りつぶし
plt.fill_between(df.index, df.iloc[:, 3], df.iloc[:, 4],
color='mistyrose', label='平均予測値の95%信頼区間', zorder=4)
# 個別予測値の95%信頼区間の塗りつぶし
plt.fill_between(df.index, df.iloc[:, 5], df.iloc[:, 6],
color='gray', alpha=0.1, label='個別予測値の95%信頼区間')
# 修飾
plt.title('予測値・標準誤差・平均予測値の95%信頼区間・個別予測値の95%信頼区間の違い')
plt.xlabel('データ番号', fontsize=12)
plt.ylabel('目的変数：配向度', fontsize=12)
plt.legend()
plt.show()【実行結果】
誤差や区間の幅感、イメージできましたか？
偏回帰係数の区間推定
テキストで取り扱われていませんが、偏回帰係数の区間推定も見ておきましょう！
ひとまず statsmodels の回帰分析結果を用いて、偏回帰係数の信頼区間を体感します。
🔢 偏回帰係数の信頼区間の計算 by statsmodels
### statsmodelsの回帰分析結果から取り出し
result1_sm.summary().tables[1]【実行結果】
右端 [ 0.025 0.975 ] が 95% 信頼区間の下端 2.5% と上端 97.5% です。
■ 偏回帰係数の信頼区間の公式
ChatGPTに教えてもらった公式です。
なお変数の記号がテキストと異なっていますのでご注意下さい。
まず標準誤差を計算、次に信頼区間を計算します。
※標準誤差の計算は前回記事と同じです。
📊 $${j}$$ 番目の偏回帰係数の推定量 $${\hat{\beta}_j}$$の標準誤差
$$
\begin{align*}
\text{SE}(\hat{\beta_j}) &= \sqrt{s^2 \left[(X^{\top}X)^{-1} \right]_{jj}}\\
\\
s^2 &= \cfrac{\sum_{i=1}^N (y_i - \hat{y}_i)^2}{N-p-1}\\
\end{align*}
$$
【変数・記号の説明】
$$
\begin{array}{l:l}
変数・記号 & 説明 \\
\hline
\\
X & 説明変数（定数項を含む）\\
\\
[(X^{\top}X)^{-1}]_{jj} & [\ ]の行列の\ j\ 行\ j\ 列成分 \\
\\
s^2 & 誤差分散\ \sigma^2\ の推定量 \\
& テキストでは誤差変動の不偏分散\ V_E \\
\\
y_i & 目的変数 \\
& テキストでは\ Y \\
\\
\hat{y}_i & 目的変数の予測値 \\
\\
\sum_{i=1}^N (y_i - \hat{y}_i)^2 & 残差平方和 \\
\\
N & 標本サイズ \\
\\
p & 説明変数の数（定数項を除く） \\
\end{array}
$$
📊 $${j}$$ 番目の偏回帰係数 $${\hat{\beta}_j}$$ の $${100(1-\alpha)\%}$$ 信頼区間
$$
\hat{\beta}_j \pm t_{\nu, 1 - \tfrac{\alpha}{2}} \cdot \text{SE}(\hat{\beta}_j)
$$
$${t_{\nu, 1 - \tfrac{\alpha}{2}}}$$ は自由度 $${\nu=N-p-1}$$ の $${t}$$ 分布の上側 $${\alpha/2}$$ 点です。
🔢 偏回帰係数の信頼区間の計算 by 公式
公式を使って信頼区間をステップ・バイ・ステップで計算しましょう。
配向度データから説明変数と目的変数を取り出して、説明変数に定数項を追加しています。
### 公式に則った計算
## 設定と準備
# 説明変数
X = data1[VARS].values
# 目的変数
y = data1[TARGET].values
# 標本サイズN、説明変数の数p
N, p = X.shape
# 説明変数の最初の列に定数項を追加
X_const = np.column_stack([np.ones(N), X])偏回帰係数を算出します。
（せっかくなので偏回帰係数の計算もします）
## 偏回帰係数の推定
beta_hat = np.linalg.inv(X_const.T @ X_const) @ X_const.T @ y
beta_hat【実行結果】
切片、温度、圧力の偏回帰係数の推定値です。
偏回帰係数の標準誤差を算出します。
公式を使っていきます！
## 偏回帰係数の標準誤差の推定
# 誤差分散の推定値
sigma2_hat = sum(result1_sm.resid**2) / (N - p - 1)
# 説明変数の行列積の逆行列
XX_inv = np.linalg.inv(X_const.T @ X_const)
# 偏回帰係数の標準誤差の推定
se_beta_hat = np.sqrt(np.diag(sigma2_hat * XX_inv))
# 結果の表示
print(se_beta_hat)【実行結果】
切片、温度、圧力の標準誤差です。
statsmodels の回帰分析結果と合っています！
いよいよ信頼区間の計算です。
## 偏回帰係数の信頼区間の推定
## 設定と準備
# 有意水準
alpha = 0.05
# 自由度N-p-1のt分布のα/2%点
t_ppt = stats.t.isf(q=alpha/2, df=N-p-1)
## 信頼区間の算出
# 下端
ci_lower = beta_hat - se_beta_hat * t_ppt
# 上端
ci_upper = beta_hat + se_beta_hat * t_ppt
# 結果の表示
print(ci_lower)
print(ci_upper)【実行結果】
上段が 95% 信頼区間の下端、下段が上端です。
左から切片、温度、圧力の順です。
statsmodels の回帰分析結果と合っています！
せっかくなので前回記事で実践した「偏回帰係数の検定統計量 $${t}$$ 値と $${p}$$ 値を」計算して、表形式に取りまとめましょう。
## t値、p値の算出
# t値の算出
t_value = beta_hat / se_beta_hat
print('t値: ', t_value)
# t値のp値の算出
p_value = stats.t.sf(abs(t_value), df=N-p-1)*2 # 両側
print('p値: ', p_value)【実行結果】
## 重回帰分析のサマリー表の作成
pd.DataFrame(
{'係数': beta_hat, '標準誤差': se_beta_hat, 't値': t_value, 'p値': p_value,
f'{alpha/2:.1%}': ci_lower, f'{1-alpha/2:.1%}': ci_upper},
index=['切片'] + VARS
).round(3)【実行結果】
statsmodels の出力を再現できました！（やったね！）
🛸ちょっと寄り道：ChatGPTおすすめのチャート２選 🛸
前回記事でChatGPTに訊いたおすすめチャートのうち
1. フォレストプロット（Coefficient Plot）
2. ブートストラップ分布のヒストグラム
に取り組みます！
📈 フォレストプロット
説明変数ごとに偏回帰係数の推定値と信頼区間を描画したプロットです。
### フォレストプロット（偏回帰係数の推定値と95% CI）の描画
## 設定
# 偏回帰係数
params = result1_sm.params[1:]
# 偏回帰係数の95%信頼区間
conf = result1_sm.conf_int().loc[VARS]
## 描画設定
# 信頼区間のエラーバーの値の算出
errors = np.array([params - conf[0], conf[1] - params])
# 説明変数のy軸上の位置の算出
y_pos = np.arange(len(VARS))
## 描画
# 描画領域の設定
plt.figure()
# 偏回帰係数の点と95%信頼区間バーの描画
plt.errorbar(params, y_pos, xerr=errors, fmt='o')
# 有意性を判断する偏回帰係数＝０の垂直線の描画
plt.axvline(0, color='tab:red', ls='--')
# 修飾
plt.xlim(round(conf[0].min() - 1), round(conf[1].max() + 1))
plt.ylim(-1, len(VARS))
plt.yticks(y_pos, VARS)
plt.xlabel('偏回帰係数', fontsize=12)
plt.title('偏回帰係数の推定値と95%信頼区間')
plt.grid(lw=0.5, alpha=0.5, axis='x')
plt.show()【実行結果】
点が 偏回帰係数の推定値、横棒が信頼区間（幅）を示しています。
【ChatGPTによる解説】
フォレスト・プロットは、変数ごとの推定係数を横軸に並べ、95％信頼区間を棒（エラーバー）で表示した縦長の棒グラフです。
一目で「どの係数がゼロをまたぐか」「どれが有意か」 がわかります。
変数名を縦軸に、推定値を横軸に配置することで、複数の係数をすっきり比較可能。
信頼区間が０を含むときに有意ではない、０を含まないときに有意である、と判断できます。
偏回帰係数の推定値を「ぱっと」理解できて、しかも、有意性の検定と同様の情報も得られます！
📈 偏回帰係数のブートストラップ分布のヒストグラム
タイトルだけでは難解ですよね。。。
ChatGPTによる解説で、ざっくりと、ブートストラップのリサンプルとは？、ヒストグラムで読み取れることとは？ を見ておきましょう！
ブートストラップ分布ヒストグラムの原理
リサンプリング
もとのデータ（10サンプル）を「重さつきくじ引き」で何度も（たとえば500回）引き直し、同じサイズのデータセットを作ります。
これを「ブートストラップサンプル」と呼びます。
再推定
各ブートストラップサンプルで重回帰モデルをあてはめ、興味のある偏回帰係数を毎回計算します。
こうして「500個の偏回帰係数」が手に入ります。
ヒストグラム化
その500個の係数を棒グラフ（ヒストグラム）にまとめると、
「このモデルでデータがちょっとだけ変わっても、係数はだいたいこんな範囲に収まるよ」という推定値のばらつきを視覚化できます。
読み取り方
山の高さ（頻度）
山が高いところほど、「その値がもっとも”よく”出る」ということ。
山の幅（広がり）
山が狭い → ばらつきが小さく、推定値に自信が持ちやすい
山が広い → 推定値が不安定で、読み取りに注意が必要
山の中心（位置）
山の中心がゼロから大きく離れている → 「変数の効果はゼロでない可能性が高い」
山がゼロをまたいでいる → 「効果がゼロという帰無仮説も否定しにくい」
応用ポイント
95％信頼区間：ヒストグラムの下端2.5％点と上端97.5％点を区切って線を引くと、教科書の信頼区間と同じ意味で使えます。
分布の形：歪みがあるときは「サンプル数が少ない」か「外れ値の影響」が考えられるため、追加の検討材料になります。
このように、ブートストラップ分布のヒストグラムは「データの揺らぎに対して係数がどうぶれるか」を直感的につかむ強力なツールです😊
500 個のブートストラップ「偏回帰係数」標本を作成して、500個の偏回帰係数を計算し、ヒストグラムを描画します！
### ブートストラップによる偏回帰係数の分布ヒストグラム
## 設定と準備
# ブートストラップ標本の個数
B = 500
# 乱数生成器
rng = np.random.default_rng(seed=5) # seed値を変更するとサンプルが変わります
# ブートストラップ標本による重回帰分析の偏回帰係数を格納する配列の初期化
boot_coef = np.zeros((B, 2))
## ブートストラップ標本で重回帰分析の実行
for i in range(B):
# ブートストラップ標本の取得 ※replace=Trueで復元抽出
idx = rng.choice(a=len(data1), size=len(data1), replace=True)
data1_b = data1.iloc[idx]
# ブートストラップ標本で重回帰分析の実行
result1_b = smf.ols(formula='配向度 ~ 温度 + 圧力', data=data1_b).fit()
# 偏回帰係数の格納
boot_coef[i, :] = result1_b.params[VARS]
## 描画処理
# 描画領域の設定
fig, axes = plt.subplots(1, 2, figsize=(8, 3.5), tight_layout=True)
# 説明変数ごとにヒストグラム描画を繰り返し処理
for j, (var, ax) in enumerate(zip(VARS, axes.flat)):
# 偏回帰係数のヒストグラムの描画
ax.hist(boot_coef[:, j], bins=20, edgecolor='white', alpha=0.7)
# 重回帰分析で推定した偏回帰係数の垂直線の描画
ax.axvline(result1_sm.params.iloc[j+1], color='tab:red', ls='--', lw=2,
label='偏回帰係数')
# 重回帰分析で推定した偏回帰係数の95%信頼区間の垂直線（オレンジ）の描画
ax.axvline(result1_sm.conf_int().loc[var].values[0], color='tab:orange',
label='95%信頼区間')
ax.axvline(result1_sm.conf_int().loc[var].values[1], color='tab:orange')
# ブートストラップ標本から推定した偏回帰係数の95%信頼区間の垂直線（グレイ）の描画
ax.axvline(np.quantile(boot_coef[:, j], 0.025), color='gray',
label='B標本の95%信頼区間')
ax.axvline(np.quantile(boot_coef[:, j], 0.975), color='gray')
# 修飾
ax.set_xlabel(f'{var}の偏回帰係数', fontsize=12)
ax.set_ylabel('頻度', fontsize=12)
ax.set_title(f'{var}のブートストラップ分布')
ax.legend()
plt.show()【実行結果】
ブートストラップ標本から得た偏回帰係数の分布が青いヒストグラムで表現されています！
歪みはなさそうです（たぶん）。
そして、以下の線を添えました。
・重回帰分析結果の偏回帰係数（赤い点線）
・重回帰分析結果の95% 信頼区間（オレンジ線）
・ブートストラップ分布から得た 95% 信頼区間（グレイ線）
ブートストラップ分布から得た 95% 信頼区間は、重回帰分析結果の 95% 信頼区間とまあまあ近い感じがします。
やっぱり可視化はいいものですね！
記事の最後はChatGPTが締めくくります。
今回は夜空に想いを寄せて。
📘 ChatGPTのひとこと：
夜空に瞬く星々を、望遠鏡を少しずつ動かしながら見つめるように、
今回は「予測値の信頼区間」と「偏回帰係数の信頼区間」という“光の帯”で
データの星並びをそっと浮かび上がらせました。
信頼区間付きプロットは、まるで星座を描くように、
予測と係数の揺らぎをやさしく伝えてくれます😊
次回は書籍の分析例を手に、
同じ星空をもっと鮮明に映し出すレンズを選ぶように──
重回帰モデルの精度をさらに高める新章へと進みます。
そしてまた一緒に、静かな夜の学びを深めましょう✨
夜風に吹かれながら星図を広げて、
データの小さな光をひとつずつ確かめる時間を楽しみにしています。
今回の写経は以上です。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！
１．のんびり統計
統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。
２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で
書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！
３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で
書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！
４．楽しい写経 ベイズ・Python等
ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀
５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で
書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。
６．データサイエンスっぽいことを綴る
統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。
７．Python機械学習プログラミング実践記
書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。
いいなと思ったら応援しよう！
応援ありがとうございます。これからもがんばって記事を作成します！
