「数値シミュレーションで読み解く統計のしくみ」をPythonで写経 ～ 第７章7.1「回帰分析と確率モデル」
第７章「回帰分析とシミュレーション」
書籍の著者 小杉考司 先生、紀ノ定保礼 先生、清水裕士 先生
この記事は、テキスト「数値シミュレーションで読み解く統計のしくみ」第７章「回帰分析とシミュレーション」の7.1節「回帰分析と確率モデル」の Python写経活動 を取り扱います。
第７章は回帰分析を掘り下げます。
今回はパラメータリカバリを学びます。
R の基本的なプログラム記法はざっくり Python の記法と近いです。
コードの文字の細部をなぞって、R と Python の両方に接近してみましょう！
ではテキストを開いてシミュレーションの旅に出発です🚀
はじめに
テキスト「数値シミュレーションで読み解く統計のしくみ」のご紹介
このシリーズは書籍「数値シミュレーションで読み解く統計のしくみ〜Rでためしてわかる心理統計」（技術評論社、「テキスト」と呼びます）の Python写経です。
テキストは、2023年9月に発売され、副題「Rでためしてわかる心理統計」のとおり、統計処理に定評のある R の具体的なコードを用いて、心理統計の理解に役立つ数値シミュレーションを実践する素晴らしい書籍です。
引用表記
この記事は、出典に記載の書籍に掲載された文章及びコードを引用し、適宜、掲載文章とコードを改変して書いています。
【出典】
「数値シミュレーションで読み解く統計のしくみ〜Rでためしてわかる心理統計」初版第１刷、著者 小杉考司・紀ノ定保礼・清水裕士、技術評論社
記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！
7.1 回帰分析と確率モデル
この記事は Jupyter Notebook 形式（拡張子 .ipynb）でPythonコードを書きます。
概ね確率分布の特性値算出には scipy.stats を用い、乱数生成には numpy.random.generator を用いるようにしています。
主に利用するライブラリをインポートします。
### インポート
# 数値計算
import numpy as np
# 確率・統計
import scipy.stats as stats
import statsmodels.formula.api as smf # 回帰分析等の統計処理
# 描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'
7.1.1 確率モデルとしての回帰分析
座学の節です。
１つの説明変数$${x}$$と目的変数$${y}$$からなる線形モデルを定式化します。
切片$${\beta_0}$$、傾き$${\beta_1}$$、残差$${e}$$です。
残差が偶然誤差ならば、平均$${0}$$、標準偏差$${\sigma}$$の正規分布に従うと仮定できるそうです。
$$
y_i = \beta_0 + \beta_1 x_i + e_i \\
e_i \sim \text{Normal}\ (0,\ \sigma) \\
$$
上の式では$${y_i}$$が確率的に振る舞うので、次の式で表現できるそうです。
この式が回帰分析の確率モデルです。
$$
Y_i \sim \text{Normal} \ (\beta_0 + \beta_1 x_i,\ \sigma)
$$
7.1.2 シミュレーションとパラメータリカバリ
パラメータリカバリについて、テキストの説明をお借りします。
母数を係数などのモデルのパラメータ推定によって正しく復元できるかどうかのチェックを、パラメータリカバリと呼びます。
この項でのシミュレーションは、次のような感じでパラメータリカバリを実践します。
回帰分析の確率モデルに従って、モデル化したデータ生成メカニズムからデータを生成します
生成したデータについて回帰分析を行います
得られたパラメータ推定量をデータ生成メカニズムのパラメータ（真値）と比較して、パラメータを正しく復元できているかチェックします。
データ生成メカニズムに与えるパラメータは以下のとおりです。
・切片$${\beta_0 = 1}$$
・傾き$${\beta_1 = 0.5}$$
・残差の標準偏差$${\sigma=2}$$
statsmodels で回帰分析を実行します。
■ サンプルサイズ$${500}$$のパラメータリカバリの実践
データを生成して回帰分析を実行します。
### 258ページ 単回帰モデルのパラメータリカバリ n=100のケース
## 設定と準備
# サンプルサイズ
n = 500
# 推定したい値を設定（任意の値）
beta0 = 1 # 切片
beta1 = 0.5 # 傾き
sigma = 2 # 残差の標準偏差
rng = np.random.default_rng(seed=123)
## データの生成
# 説明変数の生成
x = rng.uniform(size=n, low=-1, high=1)
# 残差の生成
e = rng.normal(size=n, loc=0, scale=sigma)
# 目的変数の生成
y = beta0 + beta1 * x + e
## 統計モデルによる検証 statsmodelsのolsで回帰分析
model = smf.ols('y ~ x', data=dict(x=x, y=y)).fit()
print(model.summary())【実行結果】
切片 intercept の推定値は$${1.0719}$$、傾き x の推定値は$${0.5850}$$です。
データ生成メカニズムの切片$${1}$$、傾き$${0.5}$$に近い値となっています。
テキストによると「小数第１桁まで合っていますから、ほぼリカバリに成功したといってもいいかもしれません」です。
残差の標準偏差を確認しておきます。
# 残差の標準偏差
model.resid.std(ddof=model.nobs - model.df_resid) # ddof=2【実行結果】
■ サンプルサイズを$${20}$$に変更する場合の影響確認
サンプルサイズが小さい場合、パラメータリカバリにどのような影響があるか、確認します。
### 259ページ n = 20のケース
## 設定と準備
# サンプルサイズ
n = 20
## データの生成
rng = np.random.default_rng(seed=123)
# 説明変数の生成
x = rng.uniform(size=n, low=-1, high=1)
# 残差の生成
e = rng.normal(size=n, loc=0, scale=sigma)
# 目的変数の生成
y = beta0 + beta1 * x + e
## 統計モデルによる検証 statsmodelsのolsで回帰分析
model = smf.ols('y ~ x', data=dict(x=x, y=y)).fit()
print(model.summary())【実行結果】
切片 intercept の推定値は$${1.5290}$$、傾き x の推定値は$${0.4142}$$です。
データ生成メカニズムの切片$${1}$$、傾き$${0.5}$$からずれてしまっています。
テキストによると「サンプルサイズを変えると推定精度に影響があることは明らかです」です。
残差の標準偏差を確認しておきます。
# 残差の標準偏差
model.resid.std(ddof=model.nobs - model.df_resid) # ddof=2【実行結果】
■ サンプルサイズを$${500}$$に戻し、残差の標準偏差を$${10}$$に変更する場合の影響確認
### 260ページ n = 500, sigma = 10 のケース
# n=100のケースよりも回帰係数の推定値の精度は低い. n=500程度で母数を復元できない
## 設定と準備
# サンプルサイズ
n = 500
# 残差の標準偏差
sigma = 10
## データの生成
rng = np.random.default_rng(seed=123)
# 説明変数の生成
x = rng.uniform(size=n, low=-1, high=1)
# 残差の生成
e = rng.normal(size=n, loc=0, scale=sigma)
# 目的変数の生成
y = beta0 + beta1 * x + e
## 統計モデルによる検証 statsmodelsのolsで回帰分析
model = smf.ols('y ~ x', data=dict(x=x, y=y)).fit()
print(model.summary())【実行結果】
切片 intercept の推定値は$${1.3596}$$、傾き x の推定値は$${0.9251}$$です。
こちらもデータ生成メカニズムの切片$${1}$$、傾き$${0.5}$$からずれてしまっています。
テキストによると「サンプルサイズを大きくすれば、設定値に近づけられるかもしれませんが、少なくともサンプルサイズが$${500}$$程度では母数を復元できるわけではありません」です。
残差の標準偏差を確認しておきます。
# 残差の標準偏差
model.resid.std(ddof=model.nobs - model.df_resid) # ddof=2【実行結果】
テキストのまとめをまとめます。
「サンプルが確率的にゆらぐから」「データを大量に取得して、単純なモデルに当てはめれば正しい答えが出せるというものではない」とし、実際のデータを取得しなくともシミュレーションを活用することで「数百件のデータを生成することも、回帰分析にもほとんど時間は」かからないという、シミュレーション実施によるチェックを促しています。
今回の写経は以上です。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
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の教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。