見出し画像

「入門はじめての統計解析」をPythonで写経 Vol.8 ~ 4章「はじめての統計的検定」①母平均の検定、母分散の検定、母比率の検定

4章「はじめての統計的検定」

書籍の著者 石村貞夫 先生


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

書籍の図・表・計算を淡々とPython化する写経シリーズです。
この記事は4章の統計的検定テーマのうち、母平均・母分散・母比率の検定 に取り組みます。
ChatGPTの活用も継続してまいります!

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

いろいろな色を表すイラスト・桃(ピンク):「いらすとや」さんより

はじめに


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

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

4章 はじめての統計的検定


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

4.2 母平均の検定
4.3 母分散の検定
4.4 母比率の検定

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

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

### インポート

# 数値計算
import math                                                # python標準ライブラリ
import numpy as np
import pandas as pd

# 統計
import scipy.stats as stats
import pingouin as pg
from statsmodels.stats.proportion import proportions_ztest # 母比率のz検定
from statsmodels.stats.weightstats import ttest_ind        # 2標本のt検定

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

イントロダクション

統計的検定の学びを開始するにあたり、ChatGPTがエールを送ってくれました。


プレゼンの説得力がグッと上がるときって、「なんとなく」じゃなくて「データで根拠がある」とき。
統計的検定は、その“根拠”を与えてくれるツールです。
「本当に差があるのか?」「偶然じゃないか?」という迷いに、数字で答えをくれる。
データで迷わない判断力、ここから一緒に育てていきましょう!


「統計的検定」も難解な概念ですが、頑張って取り組んでいきましょう。
まずはChatGPTによる統計的検定の「やさしいイントロ」をご堪能ください。


🎯 統計的検定ってなに?──まずは“仮説を疑う”という考え方
「たぶんこうだろう」と思っていること(=仮説)に対して、
「そのままでいいの?」と問いかけ、データを使って判断する。
それが統計的検定の考え方です。

🍰 各検定のやさしい紹介

🔸 母平均の検定
「この商品の満足度は平均で80点と言っていいの?」と平均値が基準と違うかどうかを確かめたいとき。

🔸 母分散の検定
「品質のバラつきが基準を超えていないか?」と安定性・バラつきの大きさを評価したいとき。

🔸 母比率の検定
「このキャンペーンの反応率は前回より上がった?」と割合の変化を比べたいとき。

🧾 イントロの締めに
このように検定は「データを使って“そのままでいいの?”と問い直す」ための道具です。
本当のところはどうなのか? 一緒にじっくり考えていきましょう。


統計的検定の手続きは、ざっくり次のようになります。
詳細はぜひ📘テキスト Section 4.1 をお読みください!

1.【計画段階】調査・研究の前に決めること

  • 検定方法の選択
    → 比較したい対象や仮説の内容に合った検定を選ぶ
    (例:母平均の検定、母分散の検定、母比率の検定 など)
     ※検定方法が決まると検定統計量も決まる

  • 仮説の設定

    • 帰無仮説 $${H_0}$$:覆したい仮説(例:「差がない」など)

    • 対立仮説 $${H_1}$$:支持したい仮説(例:「差がある」など)

  • 有意水準の設定
    → 判断のしきい値(例:5%=0.05)があらかじめ決める

2.【実行段階】データを得て行うこと

  • 検定統計量の計算
    → 得られたデータをもとに、帰無仮説のもとでの検定統計量を計算する

  • 帰無仮説の判定
    → 検定統計量と棄却域、または、$${p}$$値と有意水準を比較して、「帰無仮説が棄却されるかどうか」を判断する

Section 4.2 母平均の検定

【母平均の検定のサマリー】

- 母集団が従う確率分布:正規分布
- パラメータ:母平均 $${\mu}$$
- 検定統計量の分布:自由度 $${N-1}$$ の$${t}$$ 分布に従う

【仮説】
覆したい母平均の仮説 $${\mu_0}$$ を次のように設定します。

帰無仮説 $${H_0}$$:$${\mu=\mu_0}$$

テキストより引用

支持したい仮説=対立仮説は、次の3パターンのいずれかを設定します。

① 両側検定:母平均が $${\mu_0}$$ と異なることを検定する場合
 ・対立仮説 $${H_1}$$:$${\mu \neq \mu_0}$$
② 片側検定:母平均が $${\mu_0}$$ より小さいことを検定する場合
 ・対立仮説 $${H_1}$$:$${\mu < \mu_0}$$
③ 片側検定:母平均が $${\mu_0}$$ より大きいことを検定する場合
 ・対立仮説 $${H_1}$$:$${\mu > \mu_0}$$

テキストより引用

◆ ◆ ◆

■ 正規母集団の母平均の検定:検定統計量の公式 p.146

母平均の検定における検定統計量の公式をテキストよりお借りします。

標本平均 $${\bar{x}}$$、標本分散 $${s^2}$$、標本サイズ $${N}$$、帰無仮説の母平均 $${\mu_0}$$ とすると、標本 $${\{x_1, x_2, \ldots, x_N\}}$$ による母平均 $${\mu}$$ の検定統計量

$$
T(\bar{x}, s^2, N) = \cfrac{\bar{x} - \mu_0}{\sqrt{\cfrac{s^2}{N}}}
$$

テキストの数式を引用

は自由度 $${N-1}$$ の $${t}$$ 分布に従います。

公式に則って「正規母集団の母平均の検定関数」を定義します。
確率計算には scipy.stats を利用します。

### 正規母集団の母平均の検定関数(母分散未知) p.146

# 正規母集団の母平均のt検定統計量の算出関数
def t_stat_pop_mean(x_list, mu0):
    N = len(x_list)
    x_bar = sum(x_list) / N
    s2 = sum([(x - x_bar)**2 for x in x_list]) / (N - 1)
    return (x_bar - mu0) / math.sqrt(s2 / N)

# 正規母集団の母平均のt検定関数 ※確率計算はscipy.stats利用
def pop_mean_ttest(x_list, mu0, alpha=0.05, alternative='two-sided'):
    N = len(x_list)
    t_val = t_stat_pop_mean(x_list, mu0)
    match alternative:
        case 'two-sided':
            c_val = stats.t.ppf(q=1 - alpha/2, df=N - 1)
            c_val = -c_val, c_val
            p_val = stats.t.sf(x=abs(t_val), df=N - 1) * 2
        case 'less':
            c_val = stats.t.ppf(q=alpha, df=N - 1)
            p_val = stats.t.cdf(x=t_val, df=N - 1)
        case 'greater':
            c_val = stats.t.ppf(q=1 - alpha, df=N - 1)
            p_val = stats.t.sf(x=t_val, df=N - 1)
    return {'t_value': t_val, 'c_value': c_val, 'alpha': alpha, 'p_value': p_val}

テストデータ $${[14.1, 14.0, 14.1, 14.0]}$$ について、帰無仮説 $${H_0}$$ 「母平均 $${\mu=\mu_0=14}$$」、有意水準 5%($${\alpha=0.05}$$)、両側検定で母平均の検定を実行します。

# テスト
x_list = [14.1, 14.0, 14.1, 14.0]
pop_mean_ttest(x_list, mu0=14, alpha=0.05)

【実行結果】
「有意水準 5% で有意とは言えず、帰無仮説を棄却できない」です。

出力内容を簡単に説明します。

  • t_value:検定統計量 $${T}$$ の実現値

  • c_value:棄却限界値
    帰無仮説を棄却できる検定統計量 $${T}$$ のしきい値
    $${T}$$ の実現値が棄却限界値の外のとき、帰無仮説を棄却できる

  • alpha:有意水準 $${\alpha}$$

  • p_value:$${p}$$ 値
    帰無仮説のもとで検定統計量が実現値になる確率
    $${p}$$ 値が有意水準以下(または未満)のとき、帰無仮説を棄却できる

◆ ◆ ◆

■ 母平均の検定(両側検定)例題 p.147
p.147 例)の最近の製造品から無作為抽出して測定した直径データが仕様どおりとなっているかを母平均の検定で調べます。

仮説は次のとおりです。
・帰無仮説 $${H_0}$$ 「母平均 $${\mu=15.4}$$」
・対立仮説 $${H_1}$$ 「母平均 $${\mu \neq 15.4}$$」
 ⇒両側検定:仕様より大きくても小さくてもNGだから

① 自作関数で検定実行
まずデータと帰無仮説の $${\mu_0=15.4}$$ を設定します。

### 母平均の検定(両側検定) 部品の直径 p.147

# データの設定
samples = [15.5, 15.7, 15.4, 15.4, 15.6, 15.4, 15.6, 15.5, 15.4]  # 標本
mu0 = 15.4  # 帰無仮説の母平均μ₀

続いて、有意水準 $${5\%}$$、両側検定で母平均の検定を実行します。

# 関数利用
result = pop_mean_ttest(samples, mu0=mu0)
result

【実行結果】
検定統計量 $${T=2.683}$$ は、棄却限界値(上側の閾値) $${2.306}$$ より大きい(テキストによると棄却域に入る)ので、有意水準 $${5\%}$$ で帰無仮説は棄却されます。

つまり、最近の製造品の直径は仕様と異なっています!

「NG」のマーク:「いらすとや」さんより

② 検定の可視化
$${t}$$ 分布を可視化して棄却される様子を見てみましょう。

### 図示 p.148

## 設定と準備
# 統計関連の値
N = len(samples)                  # 標本サイズ
alpha = result['alpha']           # 有意水準
lower, upper = result['c_value']  # 棄却限界値
t_val = result['t_value']         # 検定統計量
t_dist = stats.t(df=N - 1)        # t分布
# グラフ描画用設定
color = 'tab:blue'                            # 基本の色
x_min, x_max = -5, 5                          # x軸の最小値、最大値
x_val = np.linspace(x_min, x_max, 1001)       # 確率密度関数用のx軸の値
x_val_lower = np.linspace(x_min, lower, 101)  # 下側棄却域用のx軸の値
x_val_upper = np.linspace(upper, x_max, 101)  # 上側棄却域用のx軸の値

## 描画処理
# 描画領域の設定
plt.figure(figsize=(7, 3))
# t分布の確率密度関数の描画
plt.plot(x_val, t_dist.pdf(x_val), color=color, label=f'自由度({N}-1)の$t$分布')
# 下側棄却域の塗りつぶし
plt.fill_between(x_val_lower, 0, t_dist.pdf(x_val_lower),
                 color=color, alpha=0.3, label=f'有意水準{alpha}の棄却域')
# 上側棄却域の塗りつぶし
plt.fill_between(x_val_upper, 0, t_dist.pdf(x_val_upper),
                 color=color, alpha=0.3)
# 検定統計量の垂直線の描画
plt.axvline(t_val, color='tab:red', ls='--', label=f'検定統計量={t_val:.3f}')
# 修飾
plt.xticks([lower, 0, upper])
plt.legend();

【実行結果】

ベル型の曲線は自由度 $${N-1=8}$$ の $${t}$$ 分布の確率密度関数です。
両端の青い領域が有意水準 $${5\%}$$ の棄却域です。
両側検定なので両端に棄却域が $${2.5\%}$$ 分づつあります。
赤い点線が検定統計量 $${T}$$ です。

$${T}$$ が右側(上側)の棄却域に入っているので、帰無仮説が棄却されるのです。

③ Python ライブラリで検定実行
普段は定評のあるライブラリを活用しましょう。
ここでは、scipy.stats と pingouin を利用します。
いずれも引数にサンプルデータと帰無仮説の母平均 $${\mu_0}$$ を渡しています。

🖲️scipy.stats

# scipy.stats利用
stats.ttest_1samp(samples, popmean=mu0)

【実行結果】
検定統計量 $${T}$$、$${p}$$ 値、$${t}$$ 分布の自由度を得ました。

🖲️pingouin

# pingouin利用
pg.ttest(x=samples, y=mu0)

【実行結果】
主な出力は、T:検定統計量 $${T}$$、p-val:$${p}$$ 値、CI95%:95%信頼区間です。

◆ ◆ ◆

■ 母平均の検定(片側検定)例題 p.150
p.150 例)のジャガイモの収穫データです。
新しい肥料によって収穫量が以前よりも「増えている」(大きくなっている)かどうかを母平均の検定で調べます。

仮説は次のとおりです。
・帰無仮説 $${H_0}$$ 「母平均 $${\mu=41.4}$$」
・対立仮説 $${H_1}$$ 「母平均 $${\mu > 41.4}$$」
 ⇒片側検定(greater):以前よりも大きくなったかに興味がある

① 自作関数で検定実行
まずデータと帰無仮説の $${\mu_0=41.4}$$ を設定します。

### 母平均の検定(片側検定) ジャガイモの収穫量 p.150

# データの設定
samples = [42.9, 43.7, 41.2, 40.8, 41.3, 44.2]  # 標本
mu0 = 41.4  # 帰無仮説の母平均μ₀

続いて、有意水準 $${5\%}$$、片側検定(greater)で母平均の検定を実行します。

# 関数利用
result = pop_mean_ttest(samples, mu0=mu0, alternative='greater')
result

【実行結果】
検定統計量 $${T=1.615}$$ は、棄却限界値(上側の閾値) $${2.015}$$ より小さくて棄却域から外れるので、有意水準 $${5\%}$$ で帰無仮説は棄却できません。

つまり、肥料の効果があると言えません。

「NG」のマーク:「いらすとや」さんより

② 検定の可視化
$${t}$$ 分布を可視化して棄却できない様子を見てみましょう。

### 図示 p.151

## 設定と準備
# 統計関連の値
N = len(samples)                  # 標本サイズ
alpha = result['alpha']           # 有意水準
upper = result['c_value']         # 棄却限界値
t_val = result['t_value']         # 検定統計量
t_dist = stats.t(df=N - 1)        # t分布
# グラフ描画用設定
color = 'tab:blue'                            # 基本の色
x_min, x_max = -5, 5                          # x軸の最小値、最大値
x_val = np.linspace(x_min, x_max, 1001)       # 確率密度関数用のx軸の値
x_val_upper = np.linspace(upper, x_max, 101)  # 上側棄却域用のx軸の値

## 描画処理
# 描画領域の設定
plt.figure(figsize=(7, 3))
# t分布の確率密度関数の描画
plt.plot(x_val, t_dist.pdf(x_val), color=color, label=f'自由度({N}-1)の$t$分布')
# 上側棄却域の塗りつぶし
plt.fill_between(x_val_upper, 0, t_dist.pdf(x_val_upper),
                 color=color, alpha=0.3, label=f'有意水準{alpha}の棄却域')
# 検定統計量の垂直線の描画
plt.axvline(t_val, color='tab:red', ls='--', label=f'検定統計量={t_val:.3f}')
# 修飾
plt.xticks([0, upper])
plt.legend();

【実行結果】
上側の片側検定なので、青い領域の棄却域は上側に $${5\%}$$ あります。

ベル型の曲線は自由度 $${N-1=5}$$ の $${t}$$ 分布の確率密度関数です。
上側の青い領域が有意水準 $${5\%}$$ の棄却域です。
赤い点線が検定統計量 $${T}$$ です。

$${T}$$ が棄却域から外れているので、帰無仮説を棄却できないのです。

③ Python ライブラリで検定実行
scipy.stats と pingouin を利用します。
いずれも引数にサンプルデータ、帰無仮説の母平均 $${\mu_0}$$、片側検定(greater)を渡しています。

🖲️scipy.stats

# scipy.stats利用
stats.ttest_1samp(samples, popmean=mu0, alternative='greater')

【実行結果】
検定統計量 $${T}$$、$${p}$$ 値、$${t}$$ 分布の自由度を得ました。

🖲️pingouin

# pingouin利用
pg.ttest(x=samples, y=mu0, alternative='greater')

【実行結果】
主な出力は、T:検定統計量 $${T}$$、p-val:$${p}$$ 値、CI95%:95%信頼区間です。

理解度チェック 母平均の検定 p.152

商品の強度に関する母平均の検定です。
有意水準 $${5\%}$$、両側検定で、コードで淡々と解きます。

標本データ、帰無仮説の母平均 $${\mu_0}$$ を設定します。

### 母平均の検定(片側検定) ザイルの破断強度 p.152

# データの設定
samples = [4480, 4510, 4570, 4360, 4240, 4520,
           4260, 4650, 4380, 4130, 4530, 4290]  # 標本
mu0 = 4500  # 帰無仮説の母平均μ₀

自作関数で検定を実行します。

# 関数利用
result = pop_mean_ttest(samples, mu0=mu0)
result

【実行結果】
帰無仮説を棄却できませんでした。

scipy.stats で検定を実行します。

# scipy.stats利用
stats.ttest_1samp(samples, popmean=mu0)

【実行結果】

pingouin で検定を実行します。

# pingouin利用
pg.ttest(x=samples, y=mu0)

【実行結果】

$${t}$$ 分布を可視化して棄却できない様子を見てみましょう。

### 図示 p.153

## 設定と準備
# 統計関連の値
N = len(samples)                  # 標本サイズ
alpha = result['alpha']           # 有意水準
lower, upper = result['c_value']  # 棄却限界値
t_val = result['t_value']         # 検定統計量
t_dist = stats.t(df=N - 1)        # t分布
# グラフ描画用設定
color = 'tab:blue'                            # 基本の色
x_min, x_max = -5, 5                          # x軸の最小値、最大値
x_val = np.linspace(x_min, x_max, 1001)       # 確率密度関数用のx軸の値
x_val_lower = np.linspace(x_min, lower, 101)  # 下側棄却域用のx軸の値
x_val_upper = np.linspace(upper, x_max, 101)  # 上側棄却域用のx軸の値

## 描画処理
# 描画領域の設定
plt.figure(figsize=(7, 3))
# t分布の確率密度関数の描画
plt.plot(x_val, t_dist.pdf(x_val), color=color, label=f'自由度({N}-1)の$t$分布')
# 下側棄却域の塗りつぶし
plt.fill_between(x_val_lower, 0, t_dist.pdf(x_val_lower),
                 color=color, alpha=0.3, label=f'有意水準{alpha}の棄却域')
# 上側棄却域の塗りつぶし
plt.fill_between(x_val_upper, 0, t_dist.pdf(x_val_upper),
                 color=color, alpha=0.3)
# 検定統計量の垂直線の描画
plt.axvline(t_val, color='tab:red', ls='--', label=f'検定統計量={t_val:.3f}')
# 修飾
plt.xticks([lower, 0, upper])
plt.legend();

【実行結果】

Section 4.3 母分散の検定

【母分散の検定のサマリー】

- 母集団が従う確率分布:正規分布
- パラメータ:母分散 $${\sigma^2}$$
- 検定統計量の分布:自由度 $${N-1}$$ の$${\chi^2}$$ 分布に従う

【仮説】
覆したい母分散の仮説 $${\sigma_0^2}$$ を次のように設定します。

帰無仮説 $${H_0}$$:$${\sigma^2=\sigma_0^2}$$

テキストより引用

支持したい仮説=対立仮説は、次の3パターンのいずれかを設定します。

① 両側検定:母分散が $${\sigma^2_0}$$ と異なることを検定する場合
 ・対立仮説 $${H_1}$$:$${\sigma^2 \neq \sigma^2_0}$$
② 片側検定:母分散が $${\sigma^2_0}$$ より小さいことを検定する場合
 ・対立仮説 $${H_1}$$:$${\sigma^2 < \sigma^2_0}$$
③ 片側検定:母分散が $${\sigma^2_0}$$ より大きいことを検定する場合
 ・対立仮説 $${H_1}$$:$${\sigma^2 > \sigma^2_0}$$

テキストより引用

◆ ◆ ◆

■ 正規母集団の母分散の検定:検定統計量の公式 p.154
母分散の検定における検定統計量の公式をテキストよりお借りします。

標本分散 $${s^2}$$、標本サイズ $${N}$$、帰無仮説の母分散 $${\sigma^2_0}$$ とすると、標本 $${\{x_1, x_2, \ldots, x_N\}}$$ による母分散 $${\mu}$$ の検定統計量

$$
T(s^2, N) = \cfrac{(N - 1)s^2}{\sigma_0^2}
$$

テキストの数式を引用

は自由度 $${N-1}$$ の $${\chi^2}$$ 分布に従います。

公式に則って「正規母集団の母分散の検定関数」を定義します。
確率計算には scipy.stats を利用します。

### 正規母集団の母分散の検定関数 p.154 ※両側検定のp値が自信ない

# 正規母集団の母分散のχ²検定統計量の算出関数
def chi2_stat_pop_variance(x_list, sigma20):
    N = len(x_list)
    x_bar = sum(x_list) / N
    s2 = sum([(x - x_bar)**2 for x in x_list]) / (N - 1)
    return ((N - 1) * s2) / sigma20

# 正規母集団の母分散のχ²検定関数 ※確率計算はscipy.stats利用
def pop_variance_chi2test(x_list, sigma20, alpha=0.05, alternative='two-sided'):
    N = len(x_list)
    chi2_val = chi2_stat_pop_variance(x_list, sigma20)
    match alternative:
        case 'two-sided':
            c_val_lower = stats.chi2.ppf(q=alpha/2, df=N - 1)
            c_val_upper = stats.chi2.ppf(q=1 - alpha/2, df=N - 1)
            c_val = c_val_lower, c_val_upper
            half_val = stats.chi2.ppf(q=0.5, df=N - 1)
            if chi2_val <= half_val:
                p_val = stats.chi2.cdf(x=chi2_val, df=N - 1) * 2
            else:
                p_val = stats.chi2.sf(x=chi2_val, df=N - 1) * 2
        case 'less':
            c_val = stats.chi2.ppf(q=alpha, df=N - 1)
            p_val = stats.chi2.cdf(x=chi2_val, df=N - 1)
        case 'greater':
            c_val = stats.chi2.ppf(q=1 - alpha, df=N - 1)
            p_val = stats.chi2.sf(x=chi2_val, df=N - 1)
    return {'chi2_value': chi2_val, 'c_value': c_val, 'alpha': alpha,
            'p_value': p_val}

テストデータ $${[1, 2, 1, 2, 3]}$$ について、帰無仮説 $${H_0}$$ 「母分散 $${\sigma^2=\sigma^2_0=6}$$」、有意水準 5%($${\alpha=0.05}$$)、両側検定で母分散の検定を実行します。

# テスト
x_list = [1, 2, 1, 2, 3]
pop_variance_chi2test(x_list, sigma20=6, alpha=0.05)

【実行結果】
「有意水準 5% で有意となり、帰無仮説を棄却できます」です。

出力内容を簡単に説明します。

  • chi2_value:検定統計量 $${T}$$ の実現値

  • c_value:棄却限界値
    帰無仮説を棄却できる検定統計量 $${T}$$ のしきい値
    $${T}$$ の実現値が棄却限界値の外のとき、帰無仮説を棄却できる

  • alpha:有意水準 $${\alpha}$$

  • p_value:$${p}$$ 値
    帰無仮説のもとで検定統計量が実現値になる確率
    $${p}$$ 値が有意水準以下(または未満)のとき、帰無仮説を棄却できる

◆ ◆ ◆

■ 母分散の検定(片側検定)例題 p.155
p.155 例)の新しい製造方法によって製造品の直径のバラツキが改善されたか(バラツキが小さくなったか)を母分散の検定で調べます。

仮説は次のとおりです。
・帰無仮説 $${H_0}$$ 「母分散 $${\sigma^2 = 0.4^2}$$」
・対立仮説 $${H_1}$$ 「母分散 $${\sigma^2 < 0.4^2}$$」
 ⇒片側検定(less):バラツキが小さくなったかに興味がある

① 自作関数で検定実行
データと帰無仮説の $${\sigma^2_0=0.4^2}$$ を設定し、検定を実行します。

### 母分散の検定(片側検定)  リベットの直径 p.155

# 関数利用
samples = [35.2, 34.5, 34.9, 35.2, 34.8, 35.1, 34.9, 35.2, 34.9, 34.8]
result = pop_variance_chi2test(samples, sigma20=0.4**2, alpha=0.05,
                               alternative='less')
result

【実行結果】
検定統計量 $${T=2.906}$$ は、棄却限界値(下側の閾値) $${3.325}$$ より小さい(テキストによると棄却域に入る)ので、有意水準 $${5\%}$$ で帰無仮説は棄却されます。

つまり、新しい製造方法によりバラツキが小さくなりました!

「OK」のマーク:「いらすとや」さんより

② 検定の可視化
$${chi^2}$$ 分布を可視化して棄却される様子を見てみましょう。

### 図示 p.153

## 設定と準備
# 統計関連の値
N = len(samples)                  # 標本サイズ
alpha = result['alpha']           # 有意水準
lower = result['c_value']         # 棄却限界値
chi2_val = result['chi2_value']   # 検定統計量
chi2_dist = stats.chi2(df=N - 1)  # カイ二乗分布
# グラフ描画用設定
color = 'tab:blue'                            # 基本の色
x_min, x_max = 0, 25                          # x軸の最小値、最大値
x_val = np.linspace(x_min, x_max, 1001)       # 確率密度関数用のx軸の値
x_val_lower = np.linspace(x_min, lower, 101)  # 下側棄却域用のx軸の値

## 描画処理
# 描画領域の設定
plt.figure(figsize=(7, 3))
# カイ二乗分布の確率密度関数の描画
plt.plot(x_val, chi2_dist.pdf(x_val), color=color, 
         label=f'自由度({N}-1)の$\chi^2$分布')
# 下側棄却域の塗りつぶし
plt.fill_between(x_val_lower, 0, chi2_dist.pdf(x_val_lower),
                 color=color, alpha=0.3, label=f'有意水準{alpha}の棄却域')
# 検定統計量の垂直線の描画
plt.axvline(chi2_val, color='tab:red', ls='--', label=f'検定統計量={chi2_val:.3f}')
# 修飾
plt.xticks([lower, 0])
plt.legend();

【実行結果】

左に偏った曲線は自由度 $${N-1=9}$$ の $${\chi^2}$$ 分布の確率密度関数です。
下側の青い領域が有意水準 $${5\%}$$ の棄却域です。
赤い点線が検定統計量 $${T}$$ です。

$${T}$$ が棄却域に入っているので、帰無仮説が棄却されるのです。

Section 4.4 母比率の検定

【母比率の検定のサマリー】

- 母集団が従う確率分布:二項分布
- パラメータ:母比率 $${p}$$
- 検定統計量の分布:標準正規分布に従う

【仮説】
覆したい母比率の仮説 $${p_0}$$ を次のように設定します。

帰無仮説 $${H_0}$$:$${p=p_0}$$

テキストより引用

支持したい仮説=対立仮説は、次の3パターンのいずれかを設定します。

① 両側検定:母比率が $${p_0}$$ と異なることを検定する場合
 ・対立仮説 $${H_1}$$:$${p \neq p_0}$$
② 片側検定:母比率が $${p_0}$$ より小さいことを検定する場合
 ・対立仮説 $${H_1}$$:$${p < p_0}$$
③ 片側検定:母比率が $${p_0}$$ より大きいことを検定する場合
 ・対立仮説 $${H_1}$$:$${p > p_0}$$

テキストより引用

◆ ◆ ◆

■ 二項母集団の母比率の検定:検定統計量の公式 p.158
母比率の検定における検定統計量の公式をテキストよりお借りします。

標本比率 $${m/N}$$、標本サイズ $${N}$$、帰無仮説の母比率 $${p_0}$$ とすると、母比率 $${p}$$ の検定統計量

$$
T(s^2, N) = \cfrac{(N - 1)s^2}{\sigma_0^2}
$$

テキストの数式を引用

は標準正規分布に従います。

公式に則って「二項母集団の母比率の検定関数」を定義します。
確率計算には scipy.stats を利用します。

### 二項母集団の母比率の検定関数 p.158

# 二項母集団の母分散のz検定統計量の算出関数
def z_stat_pop_ratio(m, N, p0):
    p_hat = m / N
    return (p_hat - p0) / math.sqrt(p0 * (1 - p0) / N)

# 二項母集団の母比率のz検定関数 ※確率計算はscipy.stats利用
def pop_ratio_ztest(m, N, p0, alpha=0.05, alternative='two-sided'):
    z_val = z_stat_pop_ratio(m, N, p0)
    std_norm_dist = stats.norm(loc=0, scale=1)
    match alternative:
        case 'two-sided':
            c_val = std_norm_dist.ppf(q=1 - alpha/2)
            c_val = -c_val, c_val
            p_val = std_norm_dist.sf(x=abs(z_val)) * 2
        case 'less':
            c_val = std_norm_dist.ppf(q=alpha)
            p_val = std_norm_dist.cdf(x=z_val)
        case 'greater':
            c_val = std_norm_dist.ppf(q=1 - alpha)
            p_val = std_norm_dist.sf(x=z_val)
    return {'z_value': z_val, 'c_value': c_val, 'alpha': alpha, 'p_value': p_val}

標本サイズ $${N=346}$$、興味のあるカテゴリのデータの個数 $${m=29}$$ について、帰無仮説 $${H_0}$$ 「母比率 $${p=p_0=0.12}$$」、有意水準 5%($${\alpha=0.05}$$)、片側検定(less)で母比率の検定を実行します。

# テスト
pop_ratio_ztest(m=29, N=346, p0=0.12, alternative='less')

【実行結果】
「有意水準 5% で有意となり、帰無仮説を棄却できます」です。

出力内容を簡単に説明します。

  • z_value:検定統計量 $${T}$$ の実現値

  • c_value:棄却限界値
    帰無仮説を棄却できる検定統計量 $${T}$$ のしきい値
    $${T}$$ の実現値が棄却限界値の外のとき、帰無仮説を棄却できる

  • alpha:有意水準 $${\alpha}$$

  • p_value:$${p}$$ 値
    帰無仮説のもとで検定統計量が実現値になる確率
    $${p}$$ 値が有意水準以下(または未満)のとき、帰無仮説を棄却できる

💡 参考情報:Python ライブラリ活用時の留意点
母平均の検定に使えそうなライブラリを探したところ、statsmodels がヒットしましたので、試してみます。

### statsmodelsで答え合わせしようと思ったが、異なる値が算出される
z_val, p_val = proportions_ztest(count=29, nobs=346, value=0.12,
                                 alternative='smaller')
print('z_value:', z_val, ', p_value:', p_val)

【実行結果】
テキストの公式に基づく自作関数と異なる結果が出力されました。

いろいろ調べた結果、以下のような状態だと想定しました。

  • おそらく proportions_ztest は主として「2標本の母比率の差の検定」を想定

  • プールした標本比率を計算に取り入れている

  • この記事のような1標本の母比率の検定もプールした標本比率の計算を用いている。

上記の仮説を試してみました。

### 二項母集団の母比率の検定関数
# statsmodelsと同じロジック(プールした標本比率の公式を使っていると思われる)

# 二項母集団の母分散のz検定統計量の算出関数2
def z_stat_pop_ratio2(m, N, p0):
    p_hat = m / N
    return (p_hat - p0) / math.sqrt(p_hat * (1 - p_hat) / N) # ★相違点

# 二項母集団の母比率のz検定関数2 ※確率計算はscipy.stats利用
def pop_ratio_ztest2(m, N, p0, alpha=0.05, alternative='two-sided'):
    z_val = z_stat_pop_ratio2(m, N, p0)
    std_norm_dist = stats.norm(loc=0, scale=1)
    match alternative:
        case 'two-sided':
            c_val = std_norm_dist.ppf(q=1 - alpha/2)
            c_val = -c_val, c_val
            p_val = std_norm_dist.sf(x=abs(z_val)) * 2
        case 'less':
            c_val = std_norm_dist.ppf(q=alpha)
            p_val = std_norm_dist.cdf(x=z_val)
        case 'greater':
            c_val = std_norm_dist.ppf(q=1 - alpha)
            p_val = std_norm_dist.sf(x=z_val)
    return {'z_value': z_val, 'c_value': c_val, 'alpha': alpha, 'p_value': p_val}

# テスト
pop_ratio_ztest2(m=29, N=346, p0=0.12, alternative='less')

【実行結果】
仮説はあたっている可能性が高いです。

statsmodels で1標本の母比率の検定を試す際は、計算方法に相違があることを留意しましょう!

◆ ◆ ◆

■ 母比率の検定(片側検定)例題 p.159
p.159 例)の新指導法による留年率減少の効果について、母比率の検定で調べます。

仮説は次のとおりです。
・帰無仮説 $${H_0}$$ 「母比率 $${p = 0.12}$$」
・対立仮説 $${H_1}$$ 「母比率 $${p < 0.12}$$」
 ⇒片側検定(less):減少したかに興味がある

データと帰無仮説の $${\sigma^2_0=0.4^2}$$ を設定し、自作関数で検定を実行します。

### 母比率の検定(片側検定) 新指導法の効果 p.159

# 関数利用
result =pop_ratio_ztest(m=29, N=346, p0=0.12, alternative='less')
result

【実行結果】
検定統計量 $${T=-2.071}$$ は、棄却限界値(下側の閾値) $${-1.645}$$ より小さい(テキストによると棄却域に入る)ので、有意水準 $${5\%}$$ で帰無仮説は棄却されます。

つまり、新指導法により留年率が小さくなりました!

「OK」のマーク:「いらすとや」さんより

標準正規分布を可視化して棄却される様子を見てみましょう。

### 図示 p.160

## 設定と準備
# 統計関連の値
alpha = result['alpha']           # 有意水準
lower = result['c_value']         # 棄却限界値
z_val = result['z_value']         # 検定統計量
std_norm_dist = stats.norm(loc=0, scale=1)  # 標準正規分布
# グラフ描画用設定
color = 'tab:blue'                            # 基本の色
x_min, x_max = -3, 3                          # x軸の最小値、最大値
x_val = np.linspace(x_min, x_max, 1001)       # 確率密度関数用のx軸の値
x_val_lower = np.linspace(x_min, lower, 101)  # 下側棄却域用のx軸の値

## 描画処理
# 描画領域の設定
plt.figure(figsize=(7, 3))
# 標準正規分布の確率密度関数の描画
plt.plot(x_val, std_norm_dist.pdf(x_val), color=color, label=f'標準正規分布')
# 下側棄却域の塗りつぶし
plt.fill_between(x_val_lower, 0, std_norm_dist.pdf(x_val_lower),
                 color=color, alpha=0.3, label=f'有意水準{alpha}の棄却域')
# 検定統計量の垂直線の描画
plt.axvline(z_val, color='tab:red', ls='--', label=f'検定統計量={z_val:.3f}')
# 修飾
plt.xticks([lower, 0])
plt.legend();

【実行結果】

ベル型の曲線は標準正規分布の確率密度関数です。
下側の青い領域が有意水準 $${5\%}$$ の棄却域です。
赤い点線が検定統計量 $${T}$$ です。

$${T}$$ が棄却域に入っているので、帰無仮説が棄却されるのです。


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

📘 ChatGPTのひとこと:

“たまたまそう見えただけかも?”という迷いに、数字で向き合うのが統計的検定。
検定の手順をひとつずつたどることで、データに対するまなざしが少しずつ磨かれていきます。
今日のひと学びが、あなたの判断にもっと確かな自信をもたらしますように。

今回の写経は以上です。


シリーズの記事

次の記事

前の記事

目次

ブログの紹介


note で7つのシリーズ記事を書いています。
ぜひ覗いていってくださいね!

1.のんびり統計

統計検定2級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定2級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

2.実験!たのしいベイズモデリング1&2をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング2」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは R言語+Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう!

3.実験!岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験!岩波データサイエンスvol.1」の4人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう!

4.楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

5.RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん!
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

6.データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

7.Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

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

いいなと思ったら応援しよう!

ネイピア DS 応援ありがとうございます。これからもがんばって記事を作成します!

この記事が参加している募集