「入門はじめての統計解析」をPythonで写経 Vol.5 ～ ３章「はじめての統計的推定」①点推定、不偏推定量、母平均の区間推定

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

書籍の著者　石村貞夫 先生

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

書籍の図・表・計算を淡々とPython化する写経シリーズです。
３章は「これぞ統計入門！」なテーマである 点推定 と 区間推定 に取り組みます。

◆ ◆ ◆

そしてこの記事からある実験を開始いたします。
それは…

ChatGPT をアシスタントに迎えて、生成AI時代の 漫才「学びサポート」を試行してまいります！

◆ ◆ ◆

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

はじめに

このブログシリーズは、書籍「入門はじめての統計解析」（東京図書、「テキスト」と呼びます）の 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に点推定、区間推定について訊いてみます。

【ChatGPTの回答】

素敵な活動ですね！
それでは、統計入門者向けに、「点推定」「区間推定」「母平均の区間推定」をやさしい言葉で順番に説明します。

🔸点推定（てんすいてい）とは？
たとえば「全国の中学生の平均身長」を知りたいとします。
でも、全国のすべての中学生を測るのはむずかしいですよね。
だから、一部の中学生（サンプル）を選んで調べて、その結果から全体を1つの数字で予想するのが「点推定」です。

📌 かんたんに言うと
→「これくらいが全体の平均だと思う！」とひとつの数字で言い切る方法です。
例：
調べた50人の中学生の平均身長が 160cm だったら、
「全国の中学生の平均身長は 160cm と推定します」
これが点推定です。

🔸区間推定（くかんすいてい）とは？
でも、たまたま選んだ50人の身長がたまたま高め・低めだったら？
そんな心配もありますよね。
そこで「本当の平均は、このあたりにあるはず！」と、数字に幅を持たせて予想するのが「区間推定」です。

📌 かんたんに言うと
→「本当の値は、この範囲の中にあると思うよ！」という幅のある予想です。
例：
「全国の中学生の平均身長は 158cm～162cm の間にあると思う」
といった感じです。

🔸母平均の区間推定（ぼへいきんのくかんすいてい）とは？
区間推定は、いろいろな「本当の値」に使えますが、特に“母集団の平均”を範囲で推定することを「母平均の区間推定」と呼びます。

📌 かんたんに言うと
→「全体（母集団）の平均が、この範囲にあると思うよ！」という特別な区間推定です。
母平均とは、「全体の平均」のこと。
だから「母平均の区間推定」は、
「全体の平均はこの範囲にあると考えられます」と伝える推定方法です。

🧁まとめ表

必要があれば、図やイラスト、日常のたとえ話も加えてご紹介できます！
「のんびり統計」読者層に合わせて、もっとポップにしたり、専門性を加えたりもできますので、お気軽にリクエストしてくださいね。

本文の例えに「中学生の身長」を使っているにもかかわらす、まとめでは「ケーキの味」に変更する柔軟さ（気分屋さん？）が面白いですね！

会場が温まったところで本題に入ります！

Section 3.2 点推定と区間推定

テキストは「不偏統計量」からスタートします。
難しい概念がでてきました。
こういうときは・・・ChatGPTに質問します！

ChatGPTによる不偏統計量の解説です。

🔸不偏統計量（ふへんとうけいりょう）とは？
まず、「統計量（とうけいりょう）」という言葉からおさらいしましょう。

🔹統計量ってなに？
統計量とは、データから計算される数のことです。
たとえば、平均・分散・中央値・最大値など、データを使って出した数字はぜんぶ「統計量」と呼ばれます。

では、本題の 「不偏統計量」 とはなんでしょう？

✨不偏統計量：かんたんに言うと…

もっと具体的に言うと…

• 私たちは「サンプル（標本）」から、「全体（母集団）」のことを知りたいと思っています。

• サンプルから計算した平均や分散が、何度もサンプリングしたときに、母集団の値とズレが少ないものを「不偏」と呼びます。

🔹たとえば、分散の話（よくある例）
分散には次の2つの式があります：

1. 母分散のときの分母： $${n}$$

2. 標本分散のときの分母： $${n−1}$$

「なんで分母を $${n−1}$$ にするの？」という疑問がよく出てきますよね。
🔍 それは、「n-1で割った方が、本当の分散（母分散）に近くなるから」です。
→ これが 不偏性（ふへんせい） です。
✅ このとき、n-1で割った分散を「不偏分散」と呼びます。

🐣まとめ

💬 最後に一言
不偏統計量は、「いつも正解！」というわけではありませんが、何度も調べたときにズレが平均的に少なくなるという意味で「信頼できる統計量」なんです。

もし図やサンプルコード付きの解説、身近なたとえ話バージョンなども必要でしたら、お気軽にお知らせくださいね！

ChatGPTは概ね「質問の論点」「説明本文」「まとめ」「次ステップの提案」の流れで、私たちに情報提供してくれます。
「次ステップの提案」を続けたい場合には、どんどんChatGPTに問いかけてみましょう✨️
提案文章をコピペして「図やサンプルコード付きの解説をお願いします」みたいな感じで。

では「点推定」からはじめましょう！

① 母集団が正規分布 $${N(\mu, \sigma^2)}$$ の場合

標準正規分布に従う母集団からランダムに標本サンプリングして、平均・分散等を計算しましょう。

分散や標準偏差を計算できるライブラリに numpy や pandas があります。
不偏分散のように分母を $${n-1}$$ にする場合は、引数 ddof に 1 を設定してやります。
ddof は自由度の「差（delta）」的な意味合いであり、$${n-1}$$ の 1 が ddof に相当します。

### 標準正規分布乱数で仮想データを作成

# 乱数シードの固定
np.random.seed(0)
# 標準正規分布乱数を100生成
samples = stats.norm.rvs(loc=0, scale=1, size=100)
# 基本統計量の算出
print('mean:', samples.mean())
print('var :', samples.var(ddof=1))
print('std :', samples.std(ddof=1))
# データの可視化
plt.hist(samples, edgecolor='white', alpha=0.7)
plt.axvline(samples.mean(), color='tab:red', lw=3, label='平均')
plt.legend();

【実行結果】
標本の平均、不偏分散、不偏分散の標準偏差を計算し、標本のヒストグラムを描画しました。

「母平均」「母分散」「母標準偏差」に関する不偏推定量の公式をテキストからお借りし、Pythonで実装いたします。
不偏統計量を算出する関数は【Python標準縛り】で書いてます。
※ただし不偏標準偏差のガンマ関数は scipy に頼っています…

■ 母平均 $${\mu}$$ の不偏推定値の公式 p.109

$${N}$$ は標本サイズ（データの個数）、$${x_n}$$ は標本（値を持つ個々のデータ）、$${\bar{x}}$$ は標本平均です。
この母平均の不偏推定量は「どのような分布に従う母集団の場合」も同じ計算式になります。

不偏推定量の算出関数を作り、その後、サンプルデータの不偏推定量を計算する流れで進めます。

### 母集団が正規分布の不偏推定値の算出関数 p.109

# 母平均μの不偏推定値 ※どのような母集団でも共通に使える
def unbias_mean(list):
    return sum(list) / len(list)

# 母平均μの不偏推定値
unbias_mean(samples)

【実行結果】
母平均は標本平均です。

母比率 $${\mu=0}$$ の標準正規分布母集団よりサンプリングしたデータから得た母平均の不偏統計量は $${0.06}$$ となりました。
難しいことを抜きにすると、この値はデータの平均値です！

■ 母分散 $${\sigma^2}$$ の不偏推定値の公式（母平均 $${\mu}$$ が既知） p.109

### 母集団が正規分布の不偏推定値の算出関数 p.109

# 母分散σ²の不偏推定値（母平均μが既知）
def norm_unbias_variance_mu(list, mu):
    N = len(list)
    return 1/N * sum([(x - mu)**2 for x in list])

# 母分散σ²の不偏推定値（母平均μが既知）
norm_unbias_variance_mu(samples, 0)

【実行結果】
母平均 $${\mu=0}$$ にして計算しており、標本の不偏分散とは異なります。

■ 母分散 $${\sigma^2}$$ の不偏推定値の公式（母平均 $${\mu}$$ が未知） p.109

この母分散の不偏推定量は「どのような分布に従う母集団の場合」も同じ計算式になります。

### 母集団が正規分布の不偏推定値の算出関数 p.109

# 母分散σ²の不偏推定値（母平均μが未知） ※どのような母集団でも共通に使える
def unbias_variance(list):
    N = len(list)
    x_bar = sum(list) / len(list)
    return 1/(N-1) * sum([(x - x_bar)**2 for x in list])

# 母分散σ²の不偏推定値（母平均μが未知）
unbias_variance(samples)

【実行結果】
こちらは標本の不偏分散です。

母分散 $${\sigma^2=1}$$ の標準正規分布母集団よりサンプリングしたデータから得た母分散の不偏統計量は $${1.03}$$ となりました。

■ 母標準偏差の不偏推定値の公式 p.109
$${\Gamma}$$ はガンマ関数と呼ばれるものです。

### 母集団が正規分布の不偏推定値の算出関数 p.109

# 母標準偏差の不偏推定値
def norm_unbias_stddev(list):
    N = len(list)
    x_bar = sum(list) / N
    return gamma((N - 1)/2) / (gamma(N/2)*math.sqrt(2/N)) \
           * math.sqrt(sum([(x - x_bar)**2 for x in list])/N)

# 母標準偏差の不偏推定値
norm_unbias_stddev(samples)

【実行結果】

⚠️ 注意ポイント ⚠️
計算した「母標準偏差の不偏推定量」の値は「不偏分散の標準偏差」の値（$${1.0129 \cdots}$$）と一致しません。
std(ddof=1) で計算する「不偏分散の標準偏差」は「不偏標準偏差ではない」のです！

（参考）係数 $${c4}$$ を用いた不偏標準偏差
品質管理・QC分野で活用される s 管理図の係数「$${c4}$$」を用いて、不偏標準偏差を計算することができます。
次のWebサイトを参考にさせていただきました。
ありがとうございます！

### 母集団が正規分布の不偏推定値の算出関数 p.109

# 母標準偏差の不偏推定値2
# https://biolab.sakura.ne.jp/unbiased-standard-deviation.html
def norm_unbias_stddev2(list):
    N = len(list)
    x_bar = sum(list) / N
    U2 = (1 / (N - 1)) * sum([(x - x_bar)**2 for x in list])           # 不偏分散
    c4 = (math.sqrt(2 / (N - 1)) * gamma(N / 2)) / gamma((N - 1) / 2)  # c4係数
    return math.sqrt(U2) / c4

# 母標準偏差の不偏推定値2
norm_unbias_stddev2(samples)

【実行結果】
先ほど計算した「母標準偏差の不偏推定量」の値と一致しています。

② 母集団が二項分布 $${B(n, p)}$$ の場合

二項分布 $${Bin(n=1, p=0.4)}$$ に従う母集団からサンプリングしましょう。

### 二項分布乱数で仮想データを作成

# 乱数シードの固定
np.random.seed(0)
# 二項分布乱数 Bin(1, 0.4)を100生成
samples = stats.binom.rvs(n=1, p=0.4, size=100)
# 基本統計量の算出
print('mean:', samples.mean())
print('var :', samples.var(ddof=1))
print('std :', samples.std(ddof=1))
# データの可視化
plt.hist(samples, edgecolor='white', alpha=0.7)
plt.axvline(samples.mean(), color='tab:red', lw=3, label='平均')
plt.legend();

【実行結果】
標本の平均、不偏分散、不偏分散の標準偏差を計算し、標本のヒストグラムを描画しました。

サンプルデータの値は $${0}$$ か $${1}$$ です。
そして確認する母数は母比率 $${p}$$ です。
ということで、二値の比率を直感的に把握できる円グラフを描画しましょう。

# 円グラフの描画
results = collections.Counter(samples)
plt.pie(results.values(), explode=[0.1, 0], labels=results.keys(),
        autopct='%1.1f%%', startangle=90, colors=['lightblue', 'lightpink'],
        textprops=dict(fontsize=14));

【実行結果】
値 $${1}$$ の比率 $${0.38}$$ を覚えておきます。

「母比率」に関する不偏推定量の公式をテキストからお借りし、Pythonで実装いたします。

■ 母比率 $${p}$$ の不偏推定値の公式 p.109

$${N}$$ は標本サイズ（データの個数）、$${m}$$ はあるカテゴリに属する標本サイズ（データの個数）です。
サンプルデータの場合、カテゴリは $${1}$$ であり、$${m}$$ は 値 $${1}$$ のデータの個数です。

母比率の不偏推定量を計算しましょう。

# 母比率の不偏統計量 m/N

## 設定
# 母比率を算出するカテゴリ値（今回は1）
category = 1

## 計算
# 設定したカテゴリに属するデータの個数 m
m = collections.Counter(samples)[category]
# 標本サイズ N
N = len(samples)
# 母比率の不偏推定値 p_hat
p_hat = m / N

## 結果表示
print(p_hat)

【実行結果】

母比率 $${p=0.4}$$ の二項分布母集団よりサンプリングしたデータから得た母比率の不偏統計量は $${0.38}$$ となりました。

◆ ◆ ◆

③  母集団がポアソン分布 $${P(\lambda)}$$ の場合

ポアソン分布 $${Poisson(\lambda=2)}$$ に従う母集団からサンプリングしましょう。

### ポアソン分布乱数で仮想データを作成

# 乱数シードの固定
np.random.seed(0)
# ポアソン分布乱数 Poisson(2)を100生成
samples = stats.poisson.rvs(mu=2, size=100)
# 基本統計量の算出
print('mean:', samples.mean())
print('var :', samples.var(ddof=1))
print('std :', samples.std(ddof=1))
# データの可視化
plt.hist(samples, bins=8, edgecolor='white', alpha=0.7)
plt.axvline(samples.mean(), color='tab:red', lw=3, label='平均')
plt.legend();

【実行結果】
標本の平均、不偏分散、不偏分散の標準偏差を計算し、標本のヒストグラムを描画しました。

「母平均」に関する不偏推定量の公式をテキストからお借りし、Pythonで実装いたします。

■ 母平均 $${\lambda}$$ の不偏推定値の公式 p.109

正規分布母集団の母平均と同じ計算式です。
母平均の不偏推定量は標本平均 $${\bar{x}}$$ （エックス・バー）なのです。

 # 母平均λの不偏推定量
unbias_mean(samples)

【実行結果】
母平均は標本平均です。

Section 3.3 母平均の区間推定

「区間推定」に進みます！「幅」のある推定です！

① 正規母集団 $${N(\mu, \sigma^2)}$$ の母平均の区間推定

母分散 $${\sigma^2}$$ が未知の場合の母平均の区間推定の公式をテキストよりお借りします。
正規母集団から標本 $${\{x_1, x_2, \cdots, x_N\}}$$ をランダムに抽出したとき、母平均 $${\mu}$$の$${100(1-\alpha)\%}$$ 信頼区間は、次のように計算できます。

標本平均 $${\bar{x}}$$、標本分散 $${s^2}$$、データの個数 $${N}$$、自由度 $${N-1}$$ の $${t}$$ 分布の $${100 \cdot \frac{\alpha}{2}}$$ %点 $${t_{N-1}(\tfrac{\alpha}{2})}$$ です。

信頼区間の数式が難しそうですよね。
$${\alpha=0.05}$$ のとき、$${95\%}$$ 信頼区間となります。
$${\alpha=0.01}$$ のとき、$${99\%}$$ 信頼区間となります。

公式に則って信頼区間を算出する関数を定義します。
$${t}$$ 分布の %点の算出には scipy.stats を利用します。

### 母平均の区間推定（正規母集団、母分散未知） p.112

def interval_pop_mean(list, alpha):
    # 標本サイズ
    N = len(list)
    # 標本平均
    x_bar = sum(list) / N
    # 標本（不偏分散）
    s2 = (1 / (N - 1)) * sum([(x - x_bar)**2 for x in list])
    # 標準誤差
    sem = math.sqrt(s2 / N)
    # 自由度N-1のt分布の100*(1 - α/2)%点の取得 scipy.stats利用
    t_pp = stats.t.ppf(q=1 - alpha / 2, df=N - 1)
    # 戻り値：100(1-α)%信頼区間
    return x_bar - t_pp * sem, x_bar + t_pp * sem

# テスト
x = [1, 3, 6, 9]
interval_pop_mean(x, 0.05)

【実行結果】
標本 $${\{1, 3, 6, 9\}}$$ より母平均の 95% 信頼区間を算出しました。

scipy.stats で検算します。

### scipy.statsで答え合わせ

## 設定と準備
confidence = 0.95   # 信頼係数
N = len(x)          # 標本サイズ
df = N - 1          # 自由度
mean = sum(x) / N   # 平均
sem = np.std(x, ddof=1) / np.sqrt(N)  # 標準誤差

## 信頼区間の算出
stats.t.interval(confidence=confidence, df=df, loc=mean, scale=sem)

【実行結果】
検算結果は一致しました。

■ 母平均の区間推定の例題 p.116, 117
p.116 例１）の 標本 $${\{7.86, 7.89, 7.84, 7.90, 7.82\}}$$ から母平均の 99% 信頼区間を推定します。

### 溶液のpH p.116
x = [7.86, 7.89, 7.84, 7.90, 7.82]
interval_pop_mean(x, alpha=0.01)

【実行結果】
母平均の 99% 信頼区間は $${7.793 \leq \mu \leq 7.931}$$ です。

続いて、p.117 例２）の 標本サイズ $${60}$$、標本平均 $${157.77}$$、標本分散 $${26.419}$$ を用いて、95% 信頼区間を推定します。
この標本分散は不偏分散（分母が $${N-1}$$）です。
$${t}$$ 分布の %点の算出には scipy.stats を利用します。

### 60人の平均身長 p.117

# 設定と準備
N = 60                   # 標本サイズ
interval = 0.95          # 信頼係数
x_bar = 157.77           # 標本平均
s2 = 26.419              # 標本分散
sem = math.sqrt(s2 / N)  # 標準誤差
t_pp = stats.t.ppf(q=1 - (1 - interval)/2, df=N - 1)  # t_{59}(0.975)

# 95%信頼区間の算出
(x_bar - t_pp * sem, x_bar + t_pp * sem)

【実行結果】
$${t}$$分布を用いて推定した母平均の 95% 信頼区間は $${156.44 \leq \mu \leq 159.10}$$ です。

テキストは標本サイズが大きいとして、標準正規分布近似を用いて信頼区間を計算しています。

② 指数分布に従う母集団の母平均の区間推定

指数分布に従う母集団の母平均の区間推定の公式をテキスト p.118 よりお借りします。

指数分布に従う母集団の母平均 $${\mu}$$ の $${100(1-\alpha)\%}$$ 信頼区間は次のように計算できます。

$${N}$$ 個の標本を $${\{x_1, x_2, \cdots, x_N\}}$$、自由度 $${2N}$$ の $${\chi^2}$$ 分布の $${pp}$$ %点を $${\chi^2_{2N}(pp)}$$ としています。

公式に則って信頼区間を算出する関数を定義します。
$${\chi^2}$$ 分布の %点の算出には scipy.stats を利用します。

### 母平均の区間推定（母集団が指数分布に従う） p.118

def expo_interval_pop_mean(list, alpha):
    # 標本サイズ
    N = len(list)
    # 標本平均
    x_bar = sum(list) / N
    # 自由度2Nのカイ二乗分布の上側100*(α/2)%点の取得 scipy.stats利用
    chi2_pp1 = stats.chi2.isf(q=alpha / 2, df=2 * N)
    # 自由度2Nのカイ二乗分布の上側100*(1 - α/2)%点の取得 scipy.stats利用
    chi2_pp2 = stats.chi2.isf(q=1 - alpha / 2, df=2 * N)
    # 戻り値：100(1-α)%信頼区間
    return 2 * N * x_bar / chi2_pp1, 2 * N * x_bar / chi2_pp2

テストデータを作成してこの関数をテストします。
１つ先の検算で「ある程度のデータの個数」を必要としていますので、指数分布に従う 100 個の乱数を作って、テスト実行します。

### テストデータを作成してテスト実行

## テストデータの作成
# 乱数シードの固定
np.random.seed(0)
# 指数分布乱数のテストデータ（100個）生成
x = stats.expon.rvs(size=100)
print('テストデータ先頭5個:', x[:5])

## テストの実行
# αの設定
alpha = 0.05
# テスト実行
result = expo_interval_pop_mean(x, alpha=alpha)
print('95%信頼区間       :', result)

【実行結果】
テストデータから推定した母平均の 95% 信頼区間は $${0.76 \leq \mu \leq 1.13}$$ です。

では検算に進みたいのですが…

🚨 ここで問題発生！ 🚨
正規母集団のときの scipy.stats のような Python ライブラリが指数分布の場合には見つかりませんでした。
そこで、汎用的に信頼区間を推定する方法の一つ ブートストラップ法 を「緩やかな検算」に使います。

以下、趣味的ではありますが、ブートストラップ法による信頼区間の推定方法にお付き合い下さい。
（難解と感じた場合は、読み飛ばしていただいて大丈夫です）

ブートストラップ法は、テストデータから指定個数（通常、かなり大きな数）の「再標本」を「復元抽出法」で取得し、再標本から信頼区間を推定する方法です。
頭が「？？？」状態かもですが、非正規母集団に対する汎用的な信頼区間推定方法ですので、「便利な推定方法があるってこと」だけでも心にお留置きくださいね！
なお、「推定値はあくまで近似値」ですので、公式を使った推定値と完全一致するものでは無いことに留意しましょう。

まずブートストラップ法による信頼区間算出関数を作ります。
scipy.stats の bootstrap を利用します。

### ブートストラップ法による信頼区間の推定関数の作成 scipy利用

# 追加インポート
from scipy.stats import bootstrap

def ci_popmean_bootstrap(x, alpha):
    # ブートストラップによる平均の信頼区間（デフォルトは95%信頼区間）
    res = bootstrap(
        data=(x,),                     # データ
        statistic=np.mean,             # 統計量：平均
        confidence_level=1 - alpha,    # 信頼係数 1-α
        n_resamples=10000,             # ブートストラップ法の再標本の数
        method='BCa',                  # 推定方法
    )
    print(f'母平均の{(1-alpha)*100}%信頼区間（ブートストラップ）: '
          f'({res.confidence_interval.low:.8f}, '
          f'{res.confidence_interval.high:.8f})')

では、テストデータの信頼区間をブートストラップ法で推定しましょう。

# ブートストラップ方による母平均の95%信頼区間の推定
res = ci_popmean_bootstrap(x, alpha)

【実行結果】
テストデータから推定した母平均の 95% 信頼区間は $${0.76 \leq \mu \leq 1.14}$$ です。
ブートストラップ法の信頼区間とテキスト公式による信頼区間がほぼ同じになったので、検算ＯＫといたしましょう。

ちなみにブートストラップ法コードのベースはChatGPTに作ってもらいました！

③ ポアソン分布に従う母集団の母平均の区間推定

ポアソン分布に従う母集団の母平均 $${\lambda}$$ の区間推定の公式をテキスト p.118 よりお借りします。

ポアソン分布に従う母集団の母平均 $${\lambda}$$の $${100(1-\alpha)\%}$$ 信頼区間は次のように計算できます。

$${N}$$ 個の標本を $${\{x_1, x_2, \cdots, x_N\}}$$、$${d=2(x_1 + x_2 + \cdots + x_N + 1)}$$、自由度 $${d-2}$$ の $${\chi^2}$$ 分布の $${pp}$$ %点を $${\chi^2_{d-2}(pp)}$$ としています。

公式に則って信頼区間を算出する関数を定義します。
$${\chi^2}$$ 分布の %点の算出には scipy.stats を利用します。

### 母平均の区間推定（母集団がポアソン分布に従う） p.118

def poisson_interval_pop_mean(list, alpha):
    # 標本サイズ
    N = len(list)
    # 自由度に用いるd
    d = 2 * (sum(list) + 1)
    # 自由度d-2のカイ二乗分布の上側100*(1-α/2)%点の取得 scipy.stats利用
    chi2_pp1 = stats.chi2.isf(q=1 - alpha / 2, df=d - 2)
    # 自由度dのカイ二乗分布の上側100*(α/2)%点の取得 scipy.stats利用
    chi2_pp2 = stats.chi2.isf(q=alpha / 2, df=d)
    # 戻り値：100(1-α)%信頼区間
    return chi2_pp1 / (2 * N), chi2_pp2 / (2 * N)

テストデータを作成してこの関数をテストします。
１つ先の検算で「ある程度のデータの個数」を必要としていますので、ポアソン分布 $${Poisson(\lambda=2)}$$ に従う 100 個の乱数を作って、テスト実行します。

### テストデータを作成してテスト実行

## テストデータの作成
# 乱数シードの固定
np.random.seed(0)
# ポアソン分布乱数Poisson(λ=2)のテストデータ（100個）生成
x = stats.poisson.rvs(mu=2, size=100)
print('テストデータ先頭5個:', x[:5])

## テストの実行
# αの設定
alpha = 0.05
# テスト実行
result = poisson_interval_pop_mean(x, alpha=alpha)
print('95%信頼区間       :', result)

【実行結果】
テストデータから推定した母平均の 95% 信頼区間は $${1.77 \leq \mu \leq 2.34}$$ です。

こちらもブートストラップ法の推定値と比べて検証しましょう。

# ブートストラップ方による母平均の95%信頼区間の推定
res = ci_popmean_bootstrap(x, alpha)

【実行結果】
テストデータから推定した母平均の 95% 信頼区間は $${1.75 \leq \mu \leq 2.38}$$ です。
ブートストラップ法の信頼区間とテキスト公式による信頼区間がかなり近い値になったので、検算ＯＫといたしましょう。

理解度チェック 母平均の区間推定 p.119～

製品のボトル内容量に関する母平均の 99% 信頼区間を算出します。

■ 標本データを設定して、自作関数を使って信頼区間を算出

### 関数利用 p.119
# 標本データ
data = [807, 811, 801, 798, 798, 795, 803, 805, 804]
# 99% 信頼区間
interval_pop_mean(data, alpha=0.01)

【実行結果】
99% 信頼区間は $${796.85 \leq \mu \leq 808.04}$$ です。

■ テキストの計算手順に沿って実行
(1) $${\alpha / 2}$$ の算出

### テキストの計算を辿る p.119
# αの設定
alpha = 0.01
alpha / 2

【実行結果】

(2) 標本サイズ $${N}$$ の算出

# 標本サイズの算出
N = len(data)
N

【実行結果】

(3) 自由度 $${N-1}$$ の $${t}$$ 分布の上側 $${\alpha/2 \%}$$ 点の算出

# 自由度N-1のt分布の上側α/2%点の算出 scipy.stats利用
t_pp = stats.t.ppf(q=1 - alpha/2, df=N - 1)
t_pp

【実行結果】

(4) 標本平均 $${\bar{x}}$$ の算出

# 標本平均の算出
x_bar = sum(data) / N
x_bar

【実行結果】

(5) 標本分散 $${s^2}$$ の算出

# 標本分散の算出
s2 = sum((x - x_bar)**2 for x in data) / (N - 1)
s2

【実行結果】

(6) 信頼係数 99% の信頼区間の下端の算出

# 信頼係数99%の信頼区間の下端
x_bar - t_pp * math.sqrt(s2 / N)

【実行結果】

(7) 信頼係数 99% の信頼区間の上端の算出

# 信頼係数99%の信頼区間の上端
x_bar + t_pp * math.sqrt(s2 / N)

【実行結果】

99% 信頼区間は $${796.85 \leq \mu \leq 808.04}$$ です！

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

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