「数学の言葉で世界を見たら」をPython で見たら… Vol.1 統計誤差・系統誤差

第１話「不確実な情報から判断する」

書籍の著者　大栗博司 先生

書籍「数学の言葉で世界を見たら」第１話「不確実な情報から判断する」の Python写経活動記録 です。 

標本調査結果に含まれる 統計誤差（※）と 系統誤差 を学べる好事例だと感じて、Pythonで意味合いを深堀りしました。
※標本誤差、偶然誤差とも呼ばれます。

書籍の冒頭に書かれたこのテーマに魅せられて、Python 写経の道に踏み入る決意をした次第です。

では書籍を開いて数学の旅に出発です🚀

はじめに

このブログシリーズは、書籍「数学の言葉で世界を見たら」（幻冬舎）で学んだ「数学の楽しさ」を「Python 写経の形式」でご紹介いたします。

【引用表記】
この記事は、出典に記載の書籍に掲載された文章とデータを引用し、適宜、掲載文章・データを改変して書いています。
【出典】
「数学の言葉で世界を見たら」 第3刷、著者 大栗博司、幻冬舎

序 米大統領選の予想が外れた理由

学びポイント

2024年11月の大統領選は「大接戦」予想が８割を占める状況下で、トランプ氏の圧勝となりました。
先生は「大接戦との世論調査結果が８割も揃うのは統計学的に異常」という点に気づいていたそうです。

それは、集めた情報＝標本に偏りがある「系統誤差」が起きている、ということ。
系統誤差が生じている場合、調査の正確性・信頼性を大きく損なうことがあるそうです。

先生は具体的な数値を挙げて「８割」の起こりにくさを説明しています。

• 統計誤差の観点では、大接戦＝支持率が五分五分の場合に「支持率の差が $${2.5\%}$$ 以下」となる確率は およそ「６割」になるはず。

• およそ「８割」の調査会社が「支持率の差が $${2.5%}$$ 以下」と発表している。６割で起きるはずの現象が８割で起きる確率は 100 億分の１であり、調査結果には系統誤差が含まれているのではないか。

「８割」の調査会社は $${250}$$ 社中の約 $${200}$$ 社です。
では、統計誤差の「６割」、６割で起きるはずの現象が８割で起きる確率「100 億分の１」はどのように計算されたのでしょう？

書籍には計算過程の詳細が紹介されていません…
とても気になります…

そこで…
生成 AI に計算過程を推測してもらい、Python で計算を追いかけます！

この記事で用いるライブラリをインポートします。

## インポート

# 数値計算
import numpy as np

# 統計処理
import scipy.stats as stats

# 可視化
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'  # または import japanize_matplotlib

支持率が50%のときに支持率の差が2.5%以内になる確率

統計誤差の観点で得られる確率「６割」の計算に進みます。

次の条件・手順で「支持率が $${50\%}$$ のときに支持率の差が $${2.5\%}$$ 以内になる確率」を求めます。

※数学面はたぶん雑です。不確かな記述を含むかもしれません。

1️⃣ 調査の概要

2️⃣ 標準誤差 $${\text{se}}$$

3️⃣ 支持率の推定量が近似的に従う正規分布
標本サイズに着目し、支持率の推定量 $${\hat{p}}$$ が、平均 $${p}$$、標準偏差 $${\text{se}}$$ の正規分布に近似的に従うと仮定します。

4️⃣ 正規分布における支持率 $${0.5 \pm 0.025/2}$$ の確率 $${P}$$
上記 3️⃣ の正規分布（※）の累積分布関数 $${\text{cdf}()}$$ を用いて：

※Python に計算を任せるので、数表利用を前提とする標準化は不要

Python で一気に解きます！

### 実際の支持率が五分五分のときに支持率の差が2.5%以内になる確率を計算

## 1. 設定と準備
# 聴取する人数（標本サイズ） 1000人
N = 1000
# 候補Aの真の支持率（母比率） 50%
p = 0.5
# 候補AとBの支持率の差 2.5%
diff = 0.025

## 2. 標準誤差の計算
se = np.sqrt(p*(1 - p)/N)
print(f'標準誤差は {se:.4f}')

## 3. 標本比率が近似的に従う正規分布の設定 scipy.stats利用
dist = stats.norm(loc=p, scale=se)

## 4. 確率の計算
# 2.5%以内の差の下限・上限の設定 ※48.75～51.25%
lower, upper = p - diff/2, p + diff/2
# 支持率の差が2.5%以内に収まる確率の計算（正規分布近似）
prob = 1 - 2 * dist.cdf(x=lower)
print(f'確率は　　 {prob:.2%}')

【実行結果】

確率は $${57.08\% \approx 60\%}$$ となりました。

この確率を可視化しましょう。

## 確率の可視化

# 近似する正規分布の確率密度関数曲線の描画
x_line = np.linspace(0.42, 0.58, 1001)
plt.plot(x_line, dist.pdf(x_line))
# 差が2.5%以内の確率に相当する面積を塗りつぶし
x_prob = np.linspace(lower, upper, 101)
plt.fill_between(x_prob, 0, dist.pdf(x_prob), alpha=0.2,
                 label=f'差が{diff:.1%}以内の確率')
# 差が2.5%以内の確率に相当する部分の両端の垂直点線の描画
plt.vlines(lower, 0, dist.pdf(lower), lw=1, ls='--')
plt.vlines(upper, 0, dist.pdf(upper), lw=1, ls='--')
# y=0の水平点線の描画
# plt.axhline(0, color='black', lw=0.5, ls='--')
# 修飾
plt.title('実際の支持率が五分五分のときに\n'
          f'支持率の差が{diff:.1%}以内になる確率={prob:.2%}')
plt.xlabel('支持率', fontsize=12)
plt.ylabel('確率密度', fontsize=12)
plt.ylim(bottom=0)
plt.legend();

【実行結果】

青い線が正規分布 $${\text{Normal} (p, \text{se})}$$ の確率密度関数です。

青塗りの領域が「支持率が $${50\%}$$のときに支持率の差が $${2.5\%}$$ 以内になる確率」の『６割』を表しています。

６割で起きるはずの現象が８割で起きる確率

「100 億分の１の確率」を見つける旅に出発です！

しかしこちらの計算は難解です。
そんなときは…
ChatGPT に頼ります！

ChatGPT と相談した結果、次のような計算でこの確率を求めました。
※数式は ChatGPT 任せです。

ChatGPTによると...
「KLダイバージェンス」と「チェルノフの不等式」を利用します。

1️⃣ KL ダイバージェンス
成功率 $${q}$$ が $${p}$$ からずれるときの情報量（KLダイバージェンス）を

とします。

2️⃣ チェルノフの不等式
1️⃣ のとき、標本成功率 $${\hat{p}}$$ が $${q}$$ 以上になる確率に関するチェルノフの不等式は

となります。

3️⃣ KLダイバージェンスの計算
書籍のケースについて、６割：$${p=0.6}$$、８割：$${q=0.8}$$ で情報量を計算すると、

になります。

4️⃣ チェルノフの不等式に当てはめる
3️⃣ の計算結果と調査会社数 $${n=250}$$ をチェルノフの不等式に当てはめて「少なくとも８割の調査会社が "支持率の差は 2.5% 以内" と判定する確率」を求めると、

です。

【結論】
確率は 100億分の1.15 以下（$${1.15 \times 10^{-10}}$$ 以下）です。

書籍の「６割で起きるはずの現象が８割で起きる確率を計算してみると なんと 100億分の1 だった」になりました。

Python で計算します。
「KLダイバージェンス関数 D」と「チェルノフの不等式・右辺の関数 prob」を上述の数式どおりに定義して、確率を求めます。
無名関数 lambda() を使用します。

### ６割で起きるはずの現象が８割で起きる確率を計算

## 設定と準備
# KLダイバージェンス関数の定義
D = lambda q, p: q * np.log(q/p) + (1 - q) * np.log((1 - q)/(1 - p))
# チェルノフの不等式の右辺の関数定義
prob = lambda n, p, q: np.exp(-n * D(q, p))

## 確率の計算
prob(n=250, p=0.6, q=0.8)

【実行結果】
$${1.15 \text{e-}10 = 1.15 \times 10^{-10}}$$ は 100 億分の 1.15 です！

「６割で起きるはずの現象が６～８割で起きる確率」を可視化しましょう。

### ６割で起きるはずの現象が６割～８割で起きる確率の可視化

## 設定
# 標本サイズ
n = 250
# p=６割
p = 0.6
# q=６割～８割（101刻み）
qs = np.linspace(0.6, 0.8, 101)

## 確率計算 
prob_list = [prob(n=n, p=p, q=q) for q in qs]

## 描画
# 確率の曲線の描画
plt.plot(qs, prob_list)
# 修飾
plt.xlabel('確率 $q$', fontsize=12)
plt.title(f'{p*10:.0f} 割で起きるはずの現象が $q$ で起きる確率', fontsize=12)
plt.grid(lw=0.5, alpha=0.5);

【実行結果】
横軸の $${q}$$ が７割あたりで確率がほぼ $${0}$$ になっています。

書籍の数値の裏側を覗くことができて、とても楽しかったです🍀

ところで、大統領選予想に系統誤差が入り込んだ原因や、系統誤差の疑いを検知できるシグナルがあるのか、など気になることはまだまだありますよね…
ぜひ、書籍でご確認下さい！

おわり

シリーズの記事

次の記事

目次

ブログの紹介

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

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