見出し画像

RとPythonと表計算 #01 二項分布

本シリーズ記事は、R言語、Python、表計算ソフトを使った統計解析の手法を、主要な確率分布をテーマに沿ってまとめていくものです。
高等学校の数学や情報の授業におけるデータ分析・活用を想定しながら書いてまいります。

私はデータ分析に表計算ソフトやPythonを利用していたのですが、このあたりで一度、Rも勉強してみようと思い立ったのがきっかけです。
特に、私のように「Pythonと比較しながらRも学習してみよう」と思っている方にお読みいただけますと嬉しく思います。

初回は二項分布について整理していきたいと思っています。

Pythonは下記のライブラリの読込はすでに行っているものとします。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats

概要

一定の確率で起こる2つの事象のうち一方が起こる試行を複数回行ったときにできる確率の分布を二項分布(binomial distribution)といいます。

2つの事象のうち一方が起こる試行のことをベルヌーイ試行といい、一方の事象が起こる確率を「成功確率」と呼ぶことにします。

成功確率が$${ p }$$であるようなベルヌーイ試行を$${ n }$$回行ったときの成功回数$${ X }$$とすると$${ X }$$は確率変数であり、このとき$${ X }$$は二項分布$${ B(n,p) }$$に従うといいます。

確率質量関数

確率変数$${ X }$$が二項分布$${ B(n,p) }$$に従うときの確率は次の式で与えられます。

$$
P(X=r) = _nC_rp^r (1-p)^{n-r} \quad (r=0,1,2, \cdots , n)
$$

($${ p }$$ : 成功確率,$${ r }$$: 成功回数,$${ n }$$: 試行回数)

まず、次の確率をR・Python・表計算ソフトで求めます。

[例] サイコロを30回投げたとき、1の目が5回出る確率

次に、$${ r=0,1,2 \cdots , 30 }$$として、$${ r }$$を横軸に、「サイコロを30回投げたときに1の目が出る回数がちょうど$${ r }$$回以下出る確率」を縦軸にとった棒グラフを描画することにします。

R

dbinom(成功回数, 試行回数, 成功確率)

dbinom(5,30,1/6)

実行結果:0.192108131051634

x <- 0:30
n <- 30
p <- 1/6
barplot(dbinom(x,n,p), ylim=c(0, 0.20), names.arg = x ,col="blue")

実行結果:

図1-R: 確率質量関数(R)

Python

stats.binom.pmf(成功回数, 試行回数, 成功確率)

stats.binom.pmf(5,30,1/6)

実行結果:0.19210813105163368

x = np.arange(0,31,1)
n = 30
p = 1/6
pm = stats.binom.pmf(x, n, p)
plt.bar(x, pm, color = "blue")
plt.grid()
plt.show()

実行結果:

図1-P: 確率質量関数(Python)

表計算ソフト

=BINOM.DIST(成功回数, 試行回数, 成功確率, false)

図1-S1: 確率質量関数(スプレッドシート関数)
  • グラフの種類:縦棒グラフ

  • データ範囲:A1:B32

図1-S2: 確率質量関数(スプレッドシートのグラフ)

累積分布関数

確率変数$${ X }$$が二項分布$${ B(n,p) }$$に従うときとします。
このとき、$${ X }$$の累積分布関数$${ F(r) }$$は次の式で与えられます。

$$
F(r) = P(X\leqq r) = \sum_{k=0}^r  _nC_k  p^k(1-p)^{n-k}
$$

これは、成功確率が$${ p }$$であるようなベルヌーイ試行を$${ n }$$回行ったとき、成功回数が$${ r }$$回以下であるような確率を表しています。

次の確率をR・Python・表計算ソフトで求めます。

[例] サイコロを30回投げたとき、1の目が出る回数が5回以下である確率

また、$${ r=0,1,2 \cdots , 30 }$$として、$${ r }$$を横軸に、「サイコロを30回投げたときに1の目が出る回数が$${ k }$$回以下出る確率」を縦軸にとった棒グラフを描画することにします。

R

pbinom(成功回数, 試行回数, 成功確率)

pbinom(5,30,1/6)

実行結果:0.616447014690065

x <- 0:30
n <- 30
p <- 1/6
barplot(pbinom(x,n,p),names.arg = x, col="green")

実行結果:

図2-R: 累積分布関数(R)

Python

stats.binom.cdf(成功回数, 試行回数, 成功確率)

stats.binom.cdf(5,30,1/6)

実行結果:0.6164470146900646

x = np.arange(0,31,1)
n = 30
p = 1/6
cd = stats.binom.cdf(x, n, p)
plt.bar(x, cd, color = "green")
plt.grid()
plt.show()
図2-P: 累積分布関数(Python)

表計算ソフト

=BINOM.DIST(成功回数, 試行回数, 成功確率, true)

図2-S1: 累積分布関数(スプレッドシート関数)
  • グラフの種類:縦棒グラフ

  • データ範囲:A1:A32, C1: C32

図2-S2: 累積分布関数(スプレッドシートのグラフ)

分位点

累積分布関数の逆関数であるパーセントポイント関数の値のことです。
つまり、確率変数$${ X }$$が二項分布$${ B(n,p) }$$に従うとき、$${ X }$$の$${ 100\alpha }$$%点とは、$${ P(X \leqq r ) \geqq \alpha }$$となる最小の自然数$${ r }$$のことを言います。

次の例を考えます。これは有意水準5%で片側検定を行うときの棄却域を求める問題です。

[例] サイコロを30回投げたときに1の目が$${ r }$$回以上出る確率が5%未満になるような自然数$${ r }$$の最小値を求めなさい。(上側5%点

これは、$${ P(X \leqq r ) \geqq 0.95 }$$となる最小の自然数$${ r }$$を求める問題に他なりません。

R

qbinom(累積確率, 試行回数, 成功確率)

qbinom(0.95, 30, 1/6)

実行結果:9

これは$${ P(X \leqq r ) \geqq 0.95 }$$となる最小の自然数$${ r }$$の値が9であることを意味しています。

Python

stats.binom.ppf(累積確率, 試行回数, 成功確率)

stats.binom.ppf(0.95,30,1/6)

実行結果:9.0

これは$${ P(X \leqq r ) \geqq 0.95 }$$となる最小の自然数$${ r }$$の値が9であることを意味しています。

表計算ソフト

=BINOM.INV(試行回数, 成功確率, 累積確率)
下では、累積確率は「1 - 有意水準」で求めています。

図3-S: 分位数(スプレッドシート関数)

実際に、上図の列Cを見れば、8回以下の確率が約0.9494、9回以下の確率が0.9803であることが分かります。

期待値と分散

確率変数$${ X }$$が二項分布$${ B(n,p) }$$に従うとき、$${ X }$$の期待値$${ E(X) }$$と分散$${ V(X) }$$は次の式で与えられます。

$$
E(X) = np ,\quad V(X) = np(1-p)
$$

Python

Pythonでは、二項分布に従う確率変数を定義して、その確率変数の期待値や分散を次のように求めることができます。

X = stats.binom(30,1/6)    # Xは二項分布B(30,1/6)に従う
print(f"E(X) = {X.mean()}, V(X) = {X.var()} ")

実行結果:E(X) = 5.0, V(X) = 4.166666666666667



次回は「ポアソン分布」についてまとめます。
最後までお読みいただきありがとうございました。