「回帰分析から学ぶ計量経済学」をPythonで写経 ～ 第６章「機械学習への道」③ロジスティック回帰・ベイズの定理・主成分分析

第６章「機械学習への道」

書籍の著者　山澤成康 先生

この記事は、書籍「回帰分析から学ぶ計量経済学」第６章「機械学習への道」の Python写経活動 を取り扱います。

今回は第６章の３つの課題に取り組みます。
テーマは ロジスティック回帰、ベイズの定理、主成分分析 です。
では書籍を開いて回帰分析の旅に出発です🚀

はじめに

書籍「回帰分析から学ぶ計量経済学」のご紹介

このシリーズは書籍「回帰分析から学ぶ計量経済学: Excelで読み解く経済のしくみ」（オーム社、「テキスト」と呼びます）の Python 写経です。

テキストは、2023年11月に発売され、副題「Excelで読み解く経済のしくみ」のとおり、主に Excel を用いて、計量経済学を平易に学べる素晴らしい書籍です。
テキストの「はじめに」に著者の先生が執筆の動機を書かれています。

私にとって計量経済学は高嶺の花ですが、このテキストでさまざまな回帰分析のアプローチを知ることができました。
また、書籍の Excel 処理を Python に置き換える「寄り道写経」の実践を通じて、回帰分析のお気持ちに少し近づけた感じがいたします。

回帰分析に慣れ親しむのに丁度良いレベル感と内容ですので、これはぜひともブログにしたい！と思って現在に至ります。
計量経済学の色を薄め、データ分析の色を濃いめに書いてまいります！

引用表記

この記事は、出典に記載の書籍に掲載された文章と配布データを引用し、適宜、掲載文章・配布データを改変して書いています。
【出典】
「回帰分析から学ぶ計量経済学: Excelで読み解く経済のしくみ」
第１版第１刷、著者 山澤成康、オーム社

記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！

第６章 機械学習への道

この記事は第６章の以下の節を取り扱います。

記事に用いるデータは、オーム社の書籍紹介サイトからダウンロードできる Excel ファイル内のデータをもとにしてCSVファイルを作成し、data フォルダに格納しています。

第６章で用いるライブラリをインポートします。

### インポート

## 数値計算
import numpy as np
import pandas as pd

## 統計
import scipy.stats as stats

## 機械学習
# irisデータセット
from sklearn.datasets import load_iris
# データの準備
from sklearn.model_selection import train_test_split  # 学習・テストデータの分割
# 交差検証
from sklearn.model_selection import StratifiedKFold   # 層化Kfold交差検証
from sklearn.model_selection import cross_validate    # CVによる学習実行
# 分類アルゴリズム
from sklearn.linear_model import LogisticRegression   # ロジスティック回帰
from sklearn.tree import DecisionTreeClassifier       # 決定木
from sklearn.neighbors import KNeighborsClassifier    # K近傍法
from sklearn.svm import SVC                           # サポートベクターマシン
from sklearn.naive_bayes import GaussianNB            # ナイーブベイズ
from sklearn.ensemble import RandomForestClassifier   # ランダムフォレスト
from sklearn.decomposition import PCA                 # 主成分分析
from sklearn.cluster import KMeans                    # k-means法
from sklearn.cluster import AgglomerativeClustering   # 階層分析
# 評価
from sklearn.metrics import accuracy_score            # 正解率
from sklearn.metrics import classification_report     # 分類評価レポート
from sklearn.metrics import confusion_matrix          # 混同行列
from sklearn.tree import plot_tree                    # 決定木の描画

## 描画
import matplotlib.pyplot as plt
import seaborn as sns
import japanize_matplotlib

課題1 ロジスティック回帰 p.205

テキストの「模擬試験の得点とA大学の合否結果の推定式」をお借りします。
$${Y_i}$$が合格確率であり、$${50\%}$$以上のとき合格です。
$${X_i}$$が模擬試験の得点です。

合否判定を関数化して、模擬試験の得点が$${70}$$点の学生の合否判定をしてみます。

### 合否判定シミュレーション

# 合否判定関数の定義
def get_result(x):
    prob = 1 / (1 + np.exp(-(-14 + 0.21 * x)))
    result = [1 if p >= 0.5 else 0 for p in prob]
    return prob, result
    
# 模試の得点x
x = 70

# 判定
prob, result = get_result(np.array([x]))
print(f"判定: {'合格' if result[0]==1 else '否'},  確率: {prob[0]:.1%}")

【実行結果】
合格です！

模擬試験の得点ごとの合格確率を可視化してみましょう。

### ロジスティック回帰による確率と合否の関係

## データの作成
# 得点データの作成
x_vals = np.arange(50, 89, 2)
# 合格確率と合否の算出
prob, result = get_result(x_vals)

## 描画処理
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 3))
# 確率の曲線の描画（ロジスティック曲線）
ax.plot(x_vals, prob)
# 合否の散布図の描画
sns.scatterplot(x=x_vals, y=result,
                hue=result, palette=['tab:red', 'tab:green'], alpha=0.6,
                ax=ax)
# 凡例表示
handles, labels = ax.get_legend_handles_labels()
ax.legend(handles=handles, labels=['否', '合格'], title='合否')
# 修飾
ax.set(xlabel='得点', ylabel='合格確率')
ax.grid(lw=0.5)

【実行結果】
模擬試験の得点が$${67}$$点付近で合否が分かれるようです。

課題2 ベイズによる分析 p.205

メールに「当選」の文字が含まれているとき、テキストのケース１～４のメールが迷惑メールである確率をベイズの定理を用いて推定します。

テキストのデータをデータフレーム化します。

### データの作成

# 各値の設定
case1 = [100, 40, 30, 10, 50, 30, 20]
case2 = [100, 40, 30, 10, 30, 30, 0]
case3 = [100, 30, 30, 0, 50, 30, 20]
case4 = [100, 100, 50, 50, 100, 50, 50]
index = ['メール総数', '迷惑メール数', '迷惑メール数_当選含む',
         '迷惑メール数_当選含まない', '当選単語数', '当選単語数_迷惑メール',
         '当選単語数_迷惑メールでない']

# データフレーム化
df2 = pd.DataFrame(dict(ケース1=case1, ケース2=case2, ケース3=case3, ケース4=case4),
                   index=index)
display(df2)

【実行結果】

ベイズの定理で迷惑メールの事後確率（「当選」が含まれるメールが迷惑メールである確率）を算出する関数を定義します。

### ベイズによる迷惑メール確率算出関数の定義

def mail_fillter(x):
    # 事前確率
    prior = x['迷惑メール数'] / x['メール総数']
    # 尤度
    likelihood = x['迷惑メール数_当選含む'] / x['迷惑メール数']
    # 全確率
    total_prob = x['当選単語数'] / x['メール総数']
    # 事後確率
    posterior = prior * likelihood / total_prob
    # 戻り値：各確率を格納した辞書
    return dict(
        事前確率=prior, 尤度=likelihood, 全確率=total_prob, 事後確率=posterior)

ベイズの定理による迷惑メール分析を実行します。

### ベイズによる迷惑メール分析の実行
result = df2.apply(mail_fillter, axis=0, result_type='expand')
display(result)

【実行結果】
こんな感じになりました。

課題3 働く力の作成 p.206

都道府県別の「働く力」指標を主成分分析で作成するという課題です。

テキストのデータを読み込みます。

### データの読み込み
# 所定内給与: 男性・1000円, 昼夜間人口比率: ％, 完全失業率: ％

df_work = pd.read_csv('./data/06_ex3_work.csv', index_col=0)
print('df_work.shape:', df_work.shape)
display(df_work.head())

【実行結果】
所定内給与（男子）、昼夜間人口比率、完全失業率の３変数を用います。

主成分分析を実行します。scikit-learn の PCA を利用します。

### 主成分分析の実行

# 設定
n_components = 3
pc_names = [f'PC{i}' for i in range(1, n_components+1)]

# データの標準化
df_work_ss = pd.DataFrame(stats.zscore(df_work, ddof=1),
                          columns=df_work.columns, index=df_work.index)

# pcaの実行
pca = PCA(n_components=n_components)
result = pca.fit_transform(df_work_ss)

# 固有値、固有ベクトル、主成分得点、主成分負荷量のデータフレーム化
eig_val_df = pd.DataFrame({'固有値': pca.explained_variance_,
                           '寄与率': pca.explained_variance_ratio_,
                           '累積寄与率': np.cumsum(pca.explained_variance_ratio_)},
                           index=pc_names)
eig_vec_df = pd.DataFrame(pca.components_.T,
                          columns=pc_names, index=df_work_ss.columns)
pc_score_df = pd.DataFrame(result, columns=pc_names, index=df_work_ss.index)
loadings_df = eig_vec_df * eig_val_df['固有値'].apply(np.sqrt).T

# 主成分分析の結果の表示
print('主成分分析')
print('n =', pca.n_samples_, '\n')
print('相関行列の固有値解析')
display(eig_val_df.round(2))
print('固有ベクトル')
display(eig_vec_df.round(2))
print('主成分得点')
display((pc_score_df[['PC1']]
         .sort_values(by='PC1', ascending=False)
         .round(3).head(10)))
print('主成分負荷量')
display(loadings_df.round(3))

【実行結果】

■ 固有値・寄与率に関して
第１主成分で 40%、第２主成分までで 73%の累積寄与率となっています。
固有値が１超となった主成分は第１主成分のみです。
⇒テキストは第１主成分のみを利用する方向で検討を進めます。

■ 固有ベクトル（第１主成分）に関して
所定内給与および昼夜人口比率が高く、完全失業率が低い都道府県ほど、働く力が高いと考えられます。

■ 主成分得点（第１主成分）の上位都道府県に関して
働く力が大きい都道府県トップ３は「東京都」「愛知県」「大阪府」と言えそうです。

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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の教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。