拙著「Pythonでデータ分析」の第５章「単回帰分析」の更新と補足説明
最小2乗法とベイズ分析
第５章の説明（最小2乗法）に基づくスクリプトの更新と、ベイズ分析の場合のスクリプトを用意した。
誤差の2乗和を最小にするという最小２乗基準を用いると、回帰直線のパラメータは、拙著（5.1.5）式で与えられる。この式を用いた単回帰分析のPythonスクリプトをリストLS1のように用意した。スクリプトファイルは圧縮ファイルc05files.zipにまとめたので、ダウンロード展開すれば自由に利用できる。
リストLS1のスクリプトでは、入力データはExcelファイルとして用意する（図LS1）。
図LS1
１行目に変数名を入力して、データ値は2行目以降に入力する。
第１欄は、データの識別名を入力する。第２列に従属変数（
）の値、第３列に独立変数（
）の値を入力する。
単回帰式
において、独立変数は、平均値が０になるように中心化した値
を独立変数として用いている。
このとき、式（5.1.2a）における
は、中心化された独立変数の値が０のときの回帰直線の値、すなわち中心化前の独立変数の平均値
ににおける回帰直線の値を表す。
リストLS1のスクリプトを実行すると、入力データファイル名の入力が求められる。
入力データファイル(*.xlsx) =
入力データのExcelファイル名を設定する。図LS1のデータを入力ファイルとするときは、そのファイル名「気温と緯度.xlsx」を次のように入力する。
入力データファイル(*.xlsx) = 気温と緯度.xlsx
入力ファイル名を入力すると、計算が始まり、計算結果が端末に出力される。
b = [21.1
-0.6804719]
r = -0.8688500627741321
R^2 = 0.7549004315826133
bの行に［b0, b1］の形で定数と傾き
推定値が表示されている。中心化された独立変数が０のときの値が12.1、すなわち緯度の平均値における単回帰直線の値は21.1℃である。
傾きは、−0.68であり、この値は独立変数の値が１異なるときの回帰直線の値の違いを表す（Gelmanら、2021）。独立変数と従属変数の因果関係の強さを表すと解釈するためには、背景となる因果関係の理論（causal model）が必要である（Pearlら、2018）。
独立変数と従属変数の相関係数はr = -0.87であり、決定係数はR^2 = 0.75である。
分析結果のグラフ表示は、図LS2のように表示される。
図LS2
図LS2のフォームを閉じると、プログラムの実行終了である。
プログラムの実行終了後、端末に以下の出力が表示されている。
FigGraph.png was saved.
Results.txt was saved.
図LS2のグラフと計算結果の数値がそれぞれファイルに保存されたことが示されている。
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University Press.
Pearl, J., & Mackenzie, D. (2018). The book of why. Penguin books.
リストLS1 単回帰分析（最小２乗法）のPythonスクリプト（simple_regression.py）
import pandas as pd
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
import japanize_matplotlib
"""
Excelファイルからの読み込み
"""
in_file_name = input('入力データファイル(*.xlsx) = ')
xlsx = pd.ExcelFile(in_file_name)
data_xlsx = pd.read_excel(xlsx)
f = open('Results.txt', 'w')
f.write('\n入力データファイル = ' + in_file_name)
keys = data_xlsx.keys()
case_name = keys[0]
case_id = [ v for v in
data_xlsx[case_name]]
print(case_id)
y_name = keys[1]
y = [v for v in data_xlsx[y_name]]
print(y)
x_name = keys[2]
xv = [v for v in data_xlsx[x_name]]
xv = np.array(xv)
xv_c = xv - np.mean(xv)
f.write('\nデータ...')
f.write(f'\n{case_name:<15s}{y_name:<15s}{x_name:<15s}{"（中心化値）":<15s}')
for i in range(len(case_id)):
f.write(f'\n{case_id[i]:<15s}{y[i]:<15.2f}{xv[i]:<15.2f}{xv_c[i]:<15.2f}')
X = np.stack([np.ones(len(xv_c)),
xv_c]).T
print(X)
b = np.linalg.inv(X.T @ X) @ X.T @
y
# 式（5.1.5）
print('b =', b)
f.write(f'\n\n定数 = {b[0]:.3f} 傾き = {b[1]:.3f}')
r = ss.pearsonr(xv_c, y)[0]
print('r =', r)
R2 = r**2
# 決定係数
print('R^2 =', R2)
f.write(f'\n\n決定係数 = {R2:.5f} 相関係数 = {r:.5f}\n')
plt.plot(xv_c, y, 'o')
xvc_min = np.min(xv_c)
xvc_max = np.max(xv_c)
plt.plot([xvc_min, xvc_max], [b[0] +
b[1]*xvc_min, b[0] + b[1]*xvc_max],
label = f'Y = {b[0]:.2f} + {b[1]:.2}*X')
for i in range(len(case_id)):
plt.text(xv_c[i], y[i],
case_id[i])
plt.xlabel(x_name + '（中心化値）')
plt.ylabel(y_name)
plt.title(f'$R^2 = {R2:.3f}$')
plt.legend()
plt.savefig('FigGraph.png')
plt.show()
print('FigGraph.png was saved.')
f.close()
print('Results.txt was saved.')
このサイトに用意した。