「数学の言葉で世界を見たら」をPython で見たら… Vol.2 サイコロと確率

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

書籍の著者　大栗博司 先生

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

第１節「まずはサイコロを振ってみる」を Python で振ってみます。
乱数生成と可視化の実践です。

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

はじめに

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

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

１ まずはサイコロを振ってみる

学びポイント

Python で疑似的にサイコロを振って、確率を体感します。
可視化の楽しさも体感します。

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

## インポート
import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'  # または import japanize_matplotlib

どの目も同じくらい出やすいサイコロを1000回振る

書籍 p.17 の【方法Ｂ】「サイコロを実際に転がして確率を計算する」をやってみます。
１から６の値が同じ確率で出現する「乱数」を 1000 個発生させます。
「どの目も同じ確率で出現するサイコロを 1000 回振るのと同じ」だと思い込みます。

乱数生成には numpy.random の generator を用います。
次のおまじないを唱えると…

「rng」でさまざまな乱数を生成できるようになります。
今回は rng.choice() で１～６の整数を等確率で乱数発生させます。

### どの目も同じくらい出やすいサイコロを1000回振る試行

## 設定
# 試行回数
N = 1000
# 乱数生成器
rng = np.random.default_rng()

## サイコロを振って可視化
# サイコロを振る
samples = rng.choice(a=[1, 2, 3, 4, 5, 6], size=N)
# 出た目の数をカウント
uniques, counts = np.unique(samples, return_counts=True)
# 棒グラフの描画
p = plt.bar(uniques, height=counts/N, alpha=0.7)
# 棒の上に出た目の比率を表示
plt.bar_label(p, labels=counts/N)
plt.xlabel('サイコロの目', fontsize=12)
plt.ylabel('確率', fontsize=12);

【実行結果】
等確率に条件設定しましたが、0.159 ～ 0.177 の範囲で出目の確率がばらついています。
これが乱数！確率の実現値なのかもです！

実行の都度、各目の確率値が変わります。
ぜひコードを動かしてみて下さい！
ちなみに毎回同じ乱数値を発生させたい場合には、「乱数シード」を設定します。
次のコードは「42」という乱数シードを設定しています。

癖のあるサイコロの確率の実験

書籍 p.16 には「癖のあるサイコロ」があって、1000 回振ったときに１の目が 496 回出たという例があります。

１の目の出現確率が 0.496、その他の目の出現確率が均等、という条件でサイコロ・シミュレーションを行います。

### １の目の出る確率が 0.496、それ以外の目は等確率のサイコロを1000回振る試行

## 設定
# 試行回数
N = 1000
# 乱数生成器
rng = np.random.default_rng()

## サイコロを振って可視化
# サイコロを振る
samples = rng.choice(a=[1, 2, 3, 4, 5, 6], p=[0.496] + [(1-0.496)/5]*5, size=N)
# 出た目の数をカウント
uniques, counts = np.unique(samples, return_counts=True)
# 棒グラフの描画
p = plt.bar(uniques, height=counts/N, alpha=0.7)
# 棒の上に出た目の比率を表示
plt.bar_label(p, labels=counts/N)
plt.xlabel('サイコロの目', fontsize=12)
plt.ylabel('確率', fontsize=12);

【実行結果】
１の目がぴったり 0.496 にはならないようです。
他の目はおよそ 0.1 の確率ですが、バラツキがあります。

実行の都度、各目の確率値が変わります。
ぜひコードを動かしてみて下さい！

どの目も同じくらい出やすいサイコロを２個、1000回振る

書籍 p.17 の「今度は２つのサイコロを振ったときのことを考えてみよう」の実験です。
２つのサイコロは影響しあっていないので、出目は「独立」に起きます。
１のゾロ目が出るのは $${1/6 \times 1/6 = 1/36 \approx 0.0278}$$ です。

### どの目も同じくらい出やすいサイコロを２個、1000回振る試行

## 設定
# 試行回数
N = 1000
# 乱数生成器
rng = np.random.default_rng()

## サイコロを振って可視化
# サイコロを２個、1000回振る
samples1 = rng.choice(a=[1, 2, 3, 4, 5, 6], size=N)
samples2 = rng.choice(a=[1, 2, 3, 4, 5, 6], size=N)
# ２つの出た目の確率を算出 ※pandasのcrosstabでクロス集計表を作成
cross_tab = pd.crosstab(
    samples1, samples2, rownames=['サイコロ1'], colnames=['サイコロ2']
) / N
# ヒートマップの描画
fig, ax = plt.subplots(figsize=(7, 6))
sns.heatmap(cross_tab, annot=True, fmt='.4f', cmap='Blues', ax=ax)
plt.title(f'出目が均等な２つのサイコロの出目確率 1/36={1/36:.4f}');

【実行結果】
確率値の２次元ヒートマップです。
結構バラツキがあるように見えます。
「THE 確率」です。

実行の都度、各目の確率値が変わります。
ぜひコードを動かしてみて下さい！

可視化は楽しいですね～♬

おわり

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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

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