「入門はじめての統計解析」をPythonで写経 Vol.7 ～ ３章「はじめての統計的推定」③サンプルサイズ、最尤法

３章「はじめての統計的推定」

書籍の著者　石村貞夫 先生

この記事は、書籍「入門はじめての統計解析」３章「はじめての統計的推定」の Python写経活動 を取り扱います。

書籍の図・表・計算を淡々とPython化する写経シリーズです。
この記事は３章の統計的推定テーマのうち、サンプルサイズの計算、最尤法によるパラメータ推定 に取り組みます。
ChatGPTの活用も継続してまいります！

では書籍を開いて統計解析の旅に出発です🚀

はじめに

このブログシリーズは、書籍「入門はじめての統計解析」（東京図書、「テキスト」と呼びます）の Python 写経を通じて得た「統計解析の楽しさ」をご紹介します。

書籍の紹介と引用表記はリンク先の記事に掲載しています。

３章 はじめての統計的推定

この記事は３章の以下のSectionを取り扱います。

記事に用いるデータは、テキストに掲載されたデータそのものを引用しています。
データ件数の少ないものはコード上でデータを登録し、データ件数の多いものはCSVファイル化してデータを読み込みしています。

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

### インポート

# 数値計算
import math                      # python標準ライブラリ
import numpy as np
import pandas as pd
from scipy.special import gamma  # ガンマ関数

# 統計
import scipy.stats as stats

# ユーティリティ
import collections               # カウンター

# 描画
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'

イントロダクション

ChatGPTが今回テーマの概要をざっくり教えてくれました。

✏️ サンプルサイズの計算

「アンケートって何人に聞けば“信頼できる結果”になるの？」
「商品テストって、何人くらいのモニターでやればいいの？」

実はこうした“調べる人数”を決めるときに、活躍するのがサンプルサイズ計算です。
統計的に信頼できる結論を出すには、「どれくらいのデータが必要か？」を考えることがとても大切なのです。

✏️ 最尤法によるパラメータ推定

「このデータが生まれたのって、どんな理由がいちばんありそう？」

最尤法（さいゆうほう）は、そんな問いに答える方法です。
観測したデータにいちばんよく合う“確率のしくみ”を見つけて、
その中でどんな数字（パラメータ）がいちばんもっともらしいかを探します。

「最尤法」は難解な概念ですが、頑張って取り組んでいきましょう！
ではスタートです！

Section 3.6 サンプルサイズを決める？！

サンプルサイズは「データの個数 $${N}$$」のことです。
アンケートなどの調査をする前に「必要なデータの個数」を予め算出したい、というのがサンプルサイズを決めたい動機です。
「データを集めるコストを抑えて」しかも「統計的に信頼できる分析が可能」になるサンプルサイズを考えます。

テキストは「信頼区間の誤差」に注目してサンプルサイズを決定する公式から始まります。

◆ ◆ ◆

■ 母比率の区間推定におけるサンプルサイズ $${N}$$ の公式 p.128
利用率、視聴率、支持率、アンケートの「はい」の割合などの調査を思い浮かべます。
「母比率の区間推定」が活用できそうですね！

母比率の区間推定の場合の公式をテキストよりお借りします。
母比率の $${100(1-\alpha)\%}$$ 信頼区間の誤差を $${E}$$ 以内に抑えるサンプルサイズ $${N}$$は、次のように計算できます。
$${95\%}$$ 信頼区間の場合、$${\alpha=0.05}$$ です。

$${p}$$ は利用率・支持率などの母比率の想定値です。
$${z (\tfrac{\alpha}{2})}$$ は 標準正規分布の $${100\cdot\tfrac {\alpha}{2}}$$ %点です。
$${95\%}$$信頼区間の場合、$${\alpha=0.05}$$なので、$${z (\tfrac{\alpha}{2})}$$ は $${2.5\%}$$ 点です。 

母比率 $${p}$$ が予測できない場合は以下の式を用います。

安全を見て $${p(1-p)}$$ が最大になる $${p=0.5}$$ を想定して、$${0.5 \times (1-0.5) = 1/2 \times 1/2 = 1/4}$$ を用いています。

計算されるサンプルサイズ $${N}$$ は、信頼区間の幅が

になるように定まります。

ちなみにこの公式は「データが全部集まること」を前提としています。
アンケート調査の場合、配布したアンケートに答えてくれないケースもありそうです。
サンプルサイズをさらにアンケート回収率等で調整して「アンケート配布数」を定めることも考えられます。

◆

公式に則ってサンプルサイズを算出する関数を定義します。
標準正規分布の %点の算出には scipy.stats を利用します。

### 母比率の区間推定のサンプルサイズNの算出関数 p.128～130
# E: 誤差（比率）、alpha: 1-信頼係数
# p: 想定母比率（想定できないときはサンプルサイズ最大になるp=0.5を使用）

def calc_sample_size_pop_ratio(E, alpha, p=0.5):
    # 標準正規分布の上側100*(α/2)%点の取得 scipy.stats利用
    z = stats.norm.ppf(q=1 - alpha/2, loc=0, scale=1)
    return (z / E)**2 * p * (1 - p)

誤差 $${E=0.02}$$（誤差比率 0.02）、信頼係数 $${100(1-\alpha)\%}$$ の $${\alpha=0.05}$$ をテストデータにして、母比率のサンプルサイズを算出してみましょう。

# テスト
calc_sample_size_pop_ratio(E=0.02, alpha=0.05)

【実行結果】
計算結果の小数以下端数を切り上げて、最低サンプルサイズは $${2401}$$ です。

📢 ちなみに情報 📢
統計WEBさんの記事に「母比率の区間推定における必要なサンプルサイズの計算フォーム」があります。
ぜひ試してみましょう！

テストデータをインプットにして計算した結果の画面を引用いたします。

◆ ◆ ◆

■ 母比率のサンプルサイズの例題 p.130
p.130 例）の 想定母比率 $${p=0.34}$$、信頼区間 95%（$${\alpha=0.05}$$）、誤差 $${E=0.02}$$ に基づいて母比率のサンプルサイズを計算します。

最初に先ほどのサンプルサイズ算出関数を利用してサクッとサンプルサイズを計算しましょう。

### 関数利用 p.130
calc_sample_size_pop_ratio(p=0.34, E=0.02, alpha=0.05)

【実行結果】
$${2156}$$ 以上のサンプルサイズが必要です。

続いて、テキスト p.130 の計算手順に沿って１ステップづつ計算を進めて、サンプルサイズを計算しましょう。

(1) $${\alpha/2}$$ の算出

### テキストの計算を辿る p.130

# 誤差Eと想定母比率p
E, p = 0.02, 0.34

# α
alpha = 0.05
alpha / 2

【実行結果】

(2) 標準正規分布の上側$${100 \cdot \alpha/2 ％}$$点の算出

# 標準正規分布の95%点
z = stats.norm.ppf(q=1 - alpha/2, loc=0, scale=1)
z

【実行結果】

(3) サンプルサイズの算出

# サンプルサイズの算出
(z / E)**2 * p * (1 - p)

【実行結果】
$${2156}$$ 以上のサンプルサイズが必要です。

■ 母平均の区間推定におけるサンプルサイズ $${N}$$ の公式 p.130
続いて「母平均」の話に移ります！
母平均の区間推定の場合の公式をテキストよりお借りします。

母平均の $${100(1-\alpha)\%}$$ 信頼区間の誤差を $${E}$$ 以内に抑えるサンプルサイズ $${N}$$は、次のように計算できます。
母分散 $${\sigma^2}$$ の想定情報が必要です。

$${z (\tfrac{\alpha}{2})}$$ は 標準正規分布の $${100\cdot\tfrac {\alpha}{2}}$$ %点です。

計算されるサンプルサイズ $${N}$$ は、信頼区間の幅が

になるように定まります。

公式に則ってサンプルサイズを算出する関数を定義します。
標準正規分布の %点の算出には scipy.stats を利用します。

### 母平均の区間推定のサンプルサイズNの算出関数 p.130
# E: 誤差（個、cmなどのサンプルの単位）、alpha: 1-信頼係数, sigma2: 想定母分散 

def calc_sample_size_pop_mean(E, alpha, sigma2):
    # 標準正規分布の上側100*(α/2)%点の取得 scipy.stats利用
    z = stats.norm.isf(q=alpha/2, loc=0, scale=1)
    return (z / E * np.sqrt(sigma2))**2

誤差 $${E=1}$$（例：誤差 1mm）、信頼係数 $${100(1-\alpha)\%}$$ の $${\alpha=0.05}$$、母分散 $${\sigma^2=3}$$ をテストデータにして、母平均のサンプルサイズを算出してみましょう。

# テスト
calc_sample_size_pop_mean(E=1, alpha=0.05, sigma2=3)

【実行結果】
計算結果の小数以下端数を切り上げて、最低サンプルサイズは $${12}$$ です。

サンプルサイズの決定では $${N}$$ に注目しました。
次のテーマ「最尤法」では別の対象にフォーカスが移ります。

Section 3.7 最尤法とは？

■ パラメータ推定と最尤法
データがある確率分布に従っていると仮定するとき、「その分布の特徴（＝パラメータ）をデータから推定しよう！」とするのが、パラメータ推定です。
パラメータはたとえば…

• 二項分布 $${\text{Bin}(n, p)}$$ の場合

  • 推定したいパラメータは母比率 $${p}$$

• 正規分布 $${\mathcal{N}(\mu, \sigma^2)}$$ の場合

  • 推定したいパラメータは母平均 $${\mu}$$ や 母分散 $${\sigma^2}$$

パラメータは確率分布の「中心」や「ばらつき」などを決める重要な値 なのです。

最尤法（さいゆうほう）は、パラメータを推定するための代表的な方法のひとつです。
ちなみに、前回記事までに取り上げた「点推定」「区間推定」も母集団の特徴であるパラメータを推定する方法です。

◆ ◆ ◆

■ 最尤法の定義 p.137
テキストの定義をお借りします。

尤度関数（ゆうどかんすう）は「このデータが観測されたのは、このパラメータだったからかもしれない」といろんなパラメータを仮に当てはめて、それぞれのもっともらしさ（=尤度）を求める関数です。

具体的には、パラメータ $${\theta}$$ の値を仮に決めて、データ $${\{x_1, x_2, \ldots, x_N\}}$$ のひとつひとつが起こる確率（または確率密度）を掛け算（総乗）します：

この掛け算の結果＝尤度が大きいということは「そのパラメータのもとでこのデータが出る可能性が高い」ということ。
さまざまなパラメータ値を当てはめて尤度を計算し、いちばん尤度が大きくなるパラメータ値を選ぶのが、最尤法のゴールです。

◆ ◆ ◆

【閑話休題】
「最尤法のイントロ部分」（テキストからの引用部分を除く）はChatGPTと一緒に編み上げました。
ChatGPTの「直感的な」文章を生み出す力を記事原稿の質向上に役立てていきたいです！
では本題に戻ります。

◆ ◆ ◆

■ 母集団の分布のタイプが二項分布であるとわかっている場合 p.133
テキストの二項分布の尤度関数の数式をお借りします。

数式に則って尤度関数を【Python標準縛り】で定義します。

### 二項分布の尤度関数の定義

# 二項分布の尤度関数
def binom_likelihood_func(x_list, n, p):
    # 二項係数部分の掛け算
    binom_coef = math.prod(
        [math.factorial(n) / (math.factorial(x) * math.factorial(n - x))
         for x in x_list])
    # pの指数expo1と1-pの指数expo2のｎ算出
    expo1 = sum(x_list)
    expo2 = sum([n - x for x in x_list])
    # 尤度の計算
    llf = binom_coef * p**expo1 * (1 - p)**expo2
    # 戻り値：尤度
    return llf

テストデータでパラメータ $${n=8, p=0.32}$$ の場合の尤度関数を確認します。

# テスト
x_list = [4, 2, 3]
binom_likelihood_func(x_list, n=8, p=0.32)

【実行結果】
尤度関数は $${0.012}$$ です。

scipy.stats で検算します。
stats.binom.pmf() を用いて、指定したパラメータにおけるテストデータの二項分布確率を算出し、np.prod() で確率を掛け算しています。

# scipy.statsで答え合わせ
np.prod(stats.binom.pmf(x_list, n=8, p=0.32))

【実行結果】
合っていました！

■ 例 母集団が二項分布$${Bin(8,p)}$$のとき：母比率$${p}$$の最尤法 p.134
テキストの標本 $${\{4, 2, 3\}}$$、試行回数パラメータ $${n=8}$$ で、母比率パラメータ $${p}$$ を最尤法で推定します。
仮の $${p}$$ の値を設定しないと尤度関数を計算できません。
ここでは $${p}$$ の仮の値＝探索範囲を $${0.32}$$ から $${0.42}$$ まで $${0.01}$$ 間隔としました。

### 関数利用 p.134

## 設定
x_list = [4, 2, 3]  # 標本
n = 8               # 試行回数

## 初期値
# 探索するパラメータpの値
p_valus = np.arange(0.32, 0.43, 0.01)
# 結果を格納するデータフレーム
result_df = pd.DataFrame(columns=['パラメータp', '尤度関数L'])

## パラメータpの尤度関数の算出
for i, p in enumerate(p_valus):
    L = binom_likelihood_func(x_list, n=n, p=p)
    result_df.loc[i, :] = p, L

## 結果の表示:尤度最大の行をハイライト
(result_df
 .style
 .set_properties(**{'background-color': 'yellow'}, 
                 subset=pd.IndexSlice[result_df['尤度関数L'].argmax(), :])
 .format({'パラメータp': '{:.2f}', '尤度関数L': '{:.8f}'}))

【実行結果】
パラメータ $${p}$$ の探索値と対応する尤度関数 $${L}$$ の値です。
最も尤度関数の値が大きい $${p=0.38}$$ が最尤法の推定値（最尤推定量）です。

◆ ◆ ◆

■ 母集団の分布のタイプが正規分布であるとわかっている場合 p.135
テキストの正規分布の尤度関数の数式をお借りします。

数式に則って尤度関数を【Python標準縛り】で定義します。

### 正規分布の尤度関数の定義

# 正規分布の尤度関数
def norm_likelihood_func(list, mu, sigma):
    # 標本サイズの算出
    N = len(list)
    # eの指数部の計算
    expo = -1/2 * sum([((x - mu)/sigma)**2 for x in list])
    # 尤度の計算
    llf = (1 / (math.sqrt(2 * math.pi) * sigma))**N * math.e**expo
    # 戻り値：尤度
    return llf

テストデータでパラメータ $${\mu=4, \sigma=3}$$ の場合の尤度関数を確認します。
$${\sigma}$$ は標準偏差です。

# テスト
x_list = [5, 3, 4]
norm_likelihood_func(x_list, mu=4, sigma=3)

【実行結果】
尤度関数は $${0.0021}$$ です。

scipy.stats で検算します。
stats.norm.pdf() を用いて、指定したパラメータにおけるテストデータの正規分布の確率密度を算出し、np.prod() で確率を掛け算しています。

# scipy.statsで答え合わせ
np.prod(stats.norm.pdf(x_list, loc=4, scale=3))

【実行結果】
合っていました！

◆ ◆ ◆

■ 例 母集団が正規分布$${N(\mu,3^2)}$$のときのパラメータ：母平均$${\mu}$$の最尤法 p.136
テキストの標本 $${\{5, 3, 4\}}$$、分散パラメータ $${\sigma^2=3^2}$$ で、母平均パラメータ $${\mu}$$ を最尤法で推定します。
仮の $${\mu}$$ の値を設定しないと尤度関数を計算できません。
ここでは $${\mu}$$ の仮の値＝探索範囲を $${3.40}$$ から $${4.60}$$ まで $${0.02}$$ 間隔としました。

### 関数利用 p.136

## 設定
x_list = [5, 3, 4]  # 標本
sigma = 3           # 標準偏差

## 初期値
# 探索するパラメータμの値
mu_valus = np.arange(3.4, 4.7, 0.2)
# 結果を格納するデータフレーム
result_df2 = pd.DataFrame(columns=['パラメータμ', '尤度関数L'])

## パラメータμの尤度関数の算出
for i, mu in enumerate(mu_valus):
    L = norm_likelihood_func(x_list, mu=mu, sigma=sigma)
    result_df2.loc[i, :] = mu, L

## 結果の表示:尤度最大の行をハイライト
(result_df2
 .style
 .set_properties(**{'background-color': 'yellow'}, 
                 subset=pd.IndexSlice[result_df2['尤度関数L'].argmax(), :])
 .format({'パラメータμ': '{:.3f}', '尤度関数L': '{:.8f}'}))

【実行結果】
パラメータ $${\mu}$$ の探索値と対応する尤度関数 $${L}$$ の値です。
最も尤度関数の値が大きい $${\mu=4.000}$$ が最尤法の推定値（最尤推定量）です。

■ scipyの最適化関数で尤度最大化を計算する
ここまでの最尤法は、尤度を計算する「パラメータの探索範囲」を限定的に設定しました。
探索範囲を設定することはちょっと面倒ですね…

そこで、scipyの最適化関数 minimize() を利用してみましょう。
探索範囲を定めなくても尤度を最大化するパラメータを推定してくれます。
実際にはマイナスの尤度を最小化するような推定をします。

こちらのWeb記事を参考にいたしました。
ありがとうございます！

正規分布に従う母集団を想定します。
scipy の miniize() をインポートします。

# 追加インポート
from scipy.optimize import minimize

◆ ◆ ◆

１つ目は母分散既知で母平均 $${\mu}$$ を最尤法で推定します。
テキストの標本 $${\{5, 3, 4\}}$$、分散パラメータ $${\sigma^2=3^2}$$

#### 尤度最大化をscipyで最適化（最小化）してみる 【平均μの推定】
# 参考サイト: https://note.com/united_code/n/nd72d381f7617

# 尤度関数の定義
def likelihood_function_mu(params, list, sigma):
    # パラメータの初期値の取り出し
    mu = params
    # 標本サイズの算出
    N = len(list)
    # eの指数部の計算
    expo = -1/2 * sum([((x - mu)/sigma)**2 for x in list])
    # 尤度の計算
    llf = (1 / (math.sqrt(2 * math.pi) * sigma))**N * math.e**expo
    # 戻り値：尤度のマイナス（最小化するので）
    return -llf

# データの設定
x_list = [5, 3, 4]
sigma = 3           # 標準偏差

# パラメータの初期値の設定
x_bar = sum(x_list) / len(x_list)  # 標本平均
init_params = [x_bar]

# 尤度最大となるパラメータμ,σをscipyで探索
result = minimize(likelihood_function_mu, init_params, args=(x_list, sigma, ))

# パラメータ推定値の表示
print(f'平均μの推定値: {result.x[0]}')
print(f'尤度　　　　 : {-result.fun:.8f}')

【実行結果】
テキストの計算式で算出した母平均 $${\mu}$$ の推定値と尤度に一致しました。

（参考）テキストの計算式で算出した母平均 $${\mu}$$ の推定値と尤度

◆ ◆ ◆

２つ目は母分散未知とし、母平均 $${\mu}$$ と母分散 $${\sigma^2}$$（実際には母分散の標準偏差 $${\sigma}$$）を最尤法で推定します。

### 尤度最大化をscipyで最適化（最小化）してみる 【平均μと標準偏差σの推定】
# 参考サイト: https://note.com/united_code/n/nd72d381f7617

# 尤度関数の定義
def likelihood_function_mu_sigma(params, list):
    # パラメータの初期値の取り出し
    mu, sigma = params
    # 標本サイズの算出
    N = len(list)
    # eの指数部の計算
    expo = -1/2 * sum([((x - mu)/sigma)**2 for x in list])
    # 尤度の計算
    llf = (1 / (math.sqrt(2 * math.pi) * sigma))**N * math.e**expo
    # 戻り値：尤度のマイナス（最小化するので）
    return -llf

# データの設定
x_list = [5, 3, 4]

# パラメータの初期値の設定
N = len(x_list)                                                   # 標本サイズ
x_bar = sum(x_list) / N                                           # 標本平均
x_std = math.sqrt(sum([(x-x_bar)**2 for x in x_list]) / (N - 1))  # 標本標準偏差
init_params = [x_bar, x_std]

# 尤度最大となるパラメータμ,σをscipyで探索
result = minimize(likelihood_function_mu_sigma, init_params, args=(x_list,))

# パラメータ推定値の表示
mu_hat, sigma_hat = result.x
print(f'平均μの推定値　　: {mu_hat:.4f}')
print(f'標準偏差σの推定値: {sigma_hat:.4f}')
print(f'尤度　　　　　　 : {-result.fun:.8f}')

【実行結果】
標準偏差 $${\sigma}$$ は母分散既知ケースの$${3}$$ と比べて、こちらの推定値の方が小さいです。
３個のデータ点の標本分散も標準偏差も $${1}$$ ですし。
尤度は大きくなっています。

最尤法、面白かったですね！

記事の最後をChatGPTに締めくくってもらいましょう！

📘 ChatGPTのひとこと：

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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

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