「入門はじめての統計解析」をPythonで写経 Vol.9 ～ ４章「はじめての統計的検定」②２つの母平均の差の検定

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

書籍の著者　石村貞夫 先生

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

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

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

はじめに

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

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

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

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

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

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

### インポート

# 数値計算
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と一緒に練った導入文をどうぞ！

「A店とB店、どっちの接客がよかった？」「この商品の新バージョン、本当に旧バージョンより使いやすいの？」
私たちは日常の中で、2つのものを比べて判断することがたくさんあります。

統計にも、“比べて判断する”ための道具があります。
それが今回のテーマ、差の検定です。
今回は、2つのグループからデータ（＝標本）を集めて、違いがあるかどうかを見ていくというのがポイント。

「差があるのかな？」「たまたまの違いじゃない？」
そんな問いに、2つの標本をもとにしっかり答えていくのが、差の検定の役割なのです。

引き続きの「統計的検定」を頑張って取り組んでいきましょう。
まずはChatGPTによる「差の検定」の「やさしいイントロ」をご堪能ください。

前回の記事では「1つの集団」についての検定を学びました。
今回はステップアップして「2つの集団を比べる」場面に取り組みます。比較することで見えてくる“違い”に注目しながら、次の4つの検定についてやさしく見ていきましょう。

🍰 各検定のやさしい紹介

🔸2つの母平均の差の検定
→ 「A社とB社で平均年収に違いがあるのか？」など、2つの集団の“平均”に差があるかを調べる検定です。

🔸 対応のある2つの母平均の差の検定
→ 「同じ人がダイエット前と後で体重が変わったか？」など、“同じ対象を2回測った”ときの平均の差を調べます。
→ 比較の相手が“別の人”ではなく“同じ人のビフォーアフター”のようなときに使います。

✨ イントロの締め
どの検定も、「違いが本当にあるのか？ それとも偶然か？」を見極めるための道具です。
データから確かな比較を引き出す、そんな一歩を一緒に踏み出してみましょう。

Section 4.5 ２つの母平均の差の検定

２つの正規母集団の母平均の差を取り扱います。

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

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

◆ ◆ ◆

■ ２つの母平均の差の検定（等分散性を仮定）：検定統計量の公式 p.162
２つの母平均の差の検定（母分散未知・等分散性を仮定）における検定統計量の公式をテキストよりお借りします。

２つの標本 $${\{x_{11}, x_{12}, \ldots, x_{1N_1}\}}$$、$${\{x_{21}, x_{22}, \ldots, x_{2N_2}\}}$$ による２つの母平均の差 $${\mu_1 - \mu_2}$$ の検定統計量

は自由度 $${N_1+N_2-2}$$ の $${t}$$ 分布に従います。
グループ１、グループ２を示す添字 $${_{1, 2}}$$ を $${_i}$$ で置き換えると、$${\bar{x_i}}$$ は標本平均、$${N_i}$$ は標本サイズです。
$${s^2}$$ は「プールした分散」と呼ばれ、各グループの標本分散 $${s_i^2}$$ を用いて次の式で算出します。

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

### ２つの正規母集団の母平均の差の検定関数 p.162

# ２つの正規母集団の母平均の差のt検定統計量の算出関数
def t_stat_pop_mean_diff(x_list1, x_list2, mu01=0, mu02=0):
    # 各グループの標本サイズ
    N1, N2 = len(x_list1), len(x_list2)
    # 各グループの標本平均
    x_bar1, x_bar2 = sum(x_list1) / N1, sum(x_list2) / N2
    # 各グループの標本分散
    s21 = sum([(x1 - x_bar1)**2 for x1 in x_list1]) / (N1 - 1)
    s22 = sum([(x2 - x_bar2)**2 for x2 in x_list2]) / (N2 - 1)
    # s²
    s2 = ((N1 - 1)*s21 + (N2 - 1)*s22) / (N1 + N2 - 2)
    # 戻り値: t検定統計量
    return ((x_bar1 - x_bar2) - (mu01 - mu02)) / math.sqrt((1/N1 + 1/N2)*s2)

# ２つの正規母集団の母平均のt検定関数 ※確率計算はscipy.stats利用
def pop_mean_diff_ttest(
        x_list1, x_list2, mu01=0, mu02=0, alpha=0.05, alternative='two-sided'):
    # 各グループの標本サイズ
    N1, N2 = len(x_list1), len(x_list2)
    # t検定統計量の算出
    t_val = t_stat_pop_mean_diff(x_list1, x_list2, mu01, mu02)
    # 自由度(N1+N2-2)のt分布の設定
    t_dist = stats.t(df=N1 + N2 - 2)
    # 棄却限界値c_valとp値p_valの算出
    match alternative:
        case 'two-sided':
            c_val = t_dist.ppf(q=1 - alpha/2)
            c_val = -c_val, c_val
            p_val = t_dist.sf(x=abs(t_val)) * 2
        case 'less':
            c_val = t_dist.ppf(q=alpha)
            p_val = t_dist.cdf(x=t_val)
        case 'greater':
            c_val = t_dist.ppf(q=1 - alpha)
            p_val = t_dist.sf(x=t_val)
    return {'t_value': t_val, 'c_value': c_val, 'alpha': alpha, 'p_value': p_val}

テストデータ $${[1, 2, 3]}$$、$${[1.5, 1.8, 2.0]}$$ について、帰無仮説 $${H_0}$$ 「母平均の差 $${\mu_1 > \mu_2}$$」、有意水準 5%（$${\alpha=0.05}$$）、両側検定で母平均の検定を実行します。

# テスト
x_list1, x_list2 = [1, 2, 3], [1.5, 1.8, 2.0]
pop_mean_diff_ttest(x_list1, x_list2, alternative='greater')

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

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

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

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

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

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

scipy.stats で検算します。

### scipy.statsで答え合わせ
stats.ttest_ind(x_list1, x_list2, equal_var=True, alternative='greater')

【実行結果】
検定統計量と $${p}$$ 値です。
検算結果は一致しました。

◆ ◆ ◆

■ ２つの母平均の差の検定（両側検定）例題 p.164
p.164 例）の「利根川水系」と「信濃川水系」に生息するイワナの体長に差があるか、２つの母平均の差の検定で調べます。

仮説は次のとおりです。
・帰無仮説 $${H_0}$$ 「母平均の差 $${\mu_1=\mu_2}$$」
・対立仮説 $${H_1}$$ 「母平均の差 $${\mu_1 \neq \mu_2}$$」
　⇒両側検定：差の大小にかかわらず差の有無に興味がある

① 自作関数で検定実行
まず２つのデータを設定します。

### ２つの母平均の差の検定（両側検定） イワナの体長 p.164

# データ
tone = [165, 130, 182, 178, 194, 206, 160, 122, 212, 165, 247, 195]
shinano = [180, 180, 235, 270, 240, 285, 164, 152]

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

# 関数利用
result = pop_mean_diff_ttest(tone, shinano)
result

【実行結果】
検定統計量 $${T=-1.765}$$ は、棄却限界値（下側の閾値） $${-2.101}$$ より大きい（テキストによると棄却域に入らない）ので、有意水準 $${5\%}$$ で帰無仮説は棄却できません。

つまり、２つの河川のイワナの体長に差があるとは言えないです。

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

### 図示 p.166

## 設定と準備
# 統計関連の値
N1, N2 = len(tone), len(shinano)  # 標本サイズ
alpha = result['alpha']           # 有意水準
lower, upper = result['c_value']  # 棄却限界値
t_val = result['t_value']         # 検定統計量
t_dist = stats.t(df=N1 + N2 - 2)  # 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'自由度({N1}+{N2}-2)の$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+N_2-2=18}$$ の $${t}$$ 分布の確率密度関数です。
両端の青い領域が有意水準 $${5\%}$$ の棄却域です。
両側検定なので両端に棄却域が $${2.5\%}$$ 分づつあります。
赤い点線が検定統計量 $${T}$$ です。

$${T}$$ が両方の棄却域に入っていないので、帰無仮説を棄却できないのです。

③ Python ライブラリで検定実行
普段は定評のあるライブラリを活用しましょう。
ここでは、scipy.stats、pingouin、statsmodels を利用します。
いずれも引数に２つサンプルデータと等分散を仮定することを渡しています。

🖲️scipy.stats

# scipy.stats利用
stats.ttest_ind(tone, shinano, equal_var=True)

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

🖲️pingouin

# pingouin利用
pg.ttest(x=tone, y=shinano, correction=False)  # correction=True: ウェルチのt検定

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

💡 Tips を差し込みます
効果量は「差の大きさ」の参考値です。
差の検定における検定統計量や $${p}$$ 値でわかるのは「差の有無」。
「差の大きさ」は効果量を参考にするようです。
次のWeb記事が参考になります。

🖲️statsmodels

# statsmodels利用
ttest_ind(tone, shinano, usevar='pooled')

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

◆ ◆ ◆

■ ２つの母平均の差の検定（片側検定）例題 p.168
p.164 例）の冷水系である「信濃川水系」に生息するイワナが「利根川水系」のイワナよりも体長が大きいか、２つの母平均の差の検定で調べます。

仮説は次のとおりです。
・帰無仮説 $${H_0}$$ 「母平均の差 $${\mu_1=\mu_2}$$」
・対立仮説 $${H_1}$$ 「母平均の差 $${\mu_1 < \mu_2}$$」
　⇒片側検定（less）：グループ２の方が大きいかに興味がある

① 自作関数で検定実行
まず２つのデータを設定します。

### ２つの母平均の差の検定（両側検定） イワナの体長 p.168

# データ
tone = [165, 130, 182, 178, 194, 206, 160, 122, 212, 165, 247, 195]
shinano = [180, 180, 235, 270, 240, 285, 164, 152]

続いて、有意水準 $${5\%}$$、片側検定（less）で２つの母平均の差の検定を実行します。

# 関数利用
result = pop_mean_diff_ttest(tone, shinano, alternative='less')
result

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

つまり、信濃川水系のイワナの方が体長が大きいと言えます、です。

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

### 図示 p.166

## 設定と準備
# 統計関連の値
N1, N2 = len(tone), len(shinano)  # 標本サイズ
alpha = result['alpha']           # 有意水準
lower = result['c_value']         # 棄却限界値
t_val = result['t_value']         # 検定統計量
t_dist = stats.t(df=N1 + N2 - 2)  # 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軸の値

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

【実行結果】

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

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

③ Python ライブラリで検定実行
scipy.stats、pingouin、statsmodels を利用します。
いずれも引数に２つサンプルデータ、等分散を仮定すること、片側検定（less）を渡しています。

🖲️scipy.stats

# scipy.stats利用
stats.ttest_ind(tone, shinano, equal_var=True, alternative='less')

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

🖲️pingouin

# pingouin利用
pg.ttest(x=tone, y=shinano, correction=False, alternative='less')

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

🖲️statsmodels

# statsmodels利用
ttest_ind(tone, shinano, usevar='pooled', alternative='smaller')

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

■ ２つの母平均の差の検定（等分散性を仮定しない場合）：検定統計量の公式 p.171
ここからは「２つのグループのデータの母分散が異なる／母分散が等しいという仮定を置かない」場合を取り扱います。

２つの母平均の差の検定（母分散未知・等分散性を仮定しない）、通称：ウェルチの検定における検定統計量の公式をテキストよりお借りします。

２つの標本 $${\{x_{11}, x_{12}, \ldots, x_{1N_1}\}}$$、$${\{x_{21}, x_{22}, \ldots, x_{2N_2}\}}$$ による２つの母平均の差 $${\mu_1 - \mu_2}$$ の検定統計量

は自由度 $${m}$$ の $${t}$$ 分布に従います。
グループ１、グループ２を示す添字 $${_{1, 2}}$$ を $${_i}$$ で置き換えると、$${\bar{x_i}}$$ は標本平均、$${s^2}$$ は標本分散、$${N_i}$$ は標本サイズです。

自由度 $${m}$$ は次の式で近似します。
$${m}$$ が整数値でないときは最も近い整数値を $${m}$$ とします。

公式に則って「ウェルチの検定関数」を定義します。
確率計算には scipy.stats を利用します。
なおこの関数の自由度 $${m}$$ は整数値にせず、小数値のままとします。

### ２つの正規母集団の母平均の差の検定関数 p.162

# ２つの正規母集団の母平均の差のウェルチのt検定統計量の算出関数
def t_stat_welch_pop_mean_diff(list1, list2, mu01=0, mu02=0):
    # 各グループの標本サイズ
    N1, N2 = len(list1), len(list2)
    # 各グループの標本平均
    x_bar1, x_bar2 = sum(list1) / N1, sum(list2) / N2
    # 各グループの標本分散
    s21 = sum([(x1 - x_bar1)**2 for x1 in list1]) / (N1 - 1)
    s22 = sum([(x2 - x_bar2)**2 for x2 in list2]) / (N2 - 1)
    # 自由度 ※整数値に丸めていない（scipy.stats等と同様にfloat型の自由度にする）
    m = (s21/N1 + s22/N2)**2 / \
        (s21**2/(N1**2 * (N1 - 1)) + s22**2/(N2**2 * (N2 - 1)))
    # 戻り値: ウェルチのt検定統計量, 自由度m
    return ((x_bar1 - x_bar2) - (mu01 - mu02)) / (s21/N1 + s22/N2)**(1/2), m

# ２つの正規母集団の母平均のウェルチのt検定関数 ※確率計算はscipy.stats利用
def pop_mean_diff_welch_ttest(
        list1, list2, mu01=0, mu02=0, alpha=0.05, alternative='two-sided'):
    # t検定統計量と自由度の算出
    t_val, m = t_stat_welch_pop_mean_diff(list1, list2, mu01, mu02)
    # 自由度(m)のt分布の設定
    t_dist = stats.t(df=m)
    # 棄却限界値c_valとp値p_valの算出
    match alternative:
        case 'two-sided':
            c_val = t_dist.ppf(q=1 - alpha/2)
            c_val = -c_val, c_val
            p_val = t_dist.sf(x=abs(t_val)) * 2
        case 'less':
            c_val = t_dist.ppf(q=alpha)
            p_val = t_dist.cdf(x=t_val)
        case 'greater':
            c_val = t_dist.ppf(q=1 - alpha)
            p_val = t_dist.sf(x=t_val)
    return {'t_value': t_val, 'c_value': c_val, 'alpha': alpha,
            'p_value': p_val, 'df': m}

テストデータ $${[1, 2, 3]}$$、$${[1.5, 1.8, 2.0]}$$ について、帰無仮説 $${H_0}$$ 「母平均の差 $${\mu_1 > \mu_2}$$」、有意水準 5%（$${\alpha=0.05}$$）、両側検定で母平均の検定を実行します。

# テスト
list1, list2 = [1, 2, 3], [1.5, 1.8, 2.0]
pop_mean_diff_welch_ttest(list1, list2, alternative='greater')

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

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

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

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

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

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

• df：自由度

scipy.stats で検算します。

### scipy.statsで答え合わせ ※equal_var=Falseでウェルチの検定になる
stats.ttest_ind(list1, list2, equal_var=False, alternative='greater')

【実行結果】
検定統計量、$${p}$$ 値、自由度です。
検算結果は一致しました。

◆ ◆ ◆

■ ウェルチの検定（両側検定）例題 p.164
p.172 例）の「利根川水系」と「信濃川水系」に生息するイワナの体長に差があるか、ウェルチの検定で調べます。

仮説は次のとおりです。
・帰無仮説 $${H_0}$$ 「母平均の差 $${\mu_1=\mu_2}$$」
・対立仮説 $${H_1}$$ 「母平均の差 $${\mu_1 \neq \mu_2}$$」
　⇒両側検定：差の大小にかかわらず差の有無に興味がある

① 自作関数で検定実行
まず２つのデータを設定します。

### ２つの母平均の差のウェルチの検定（両側検定） イワナの体長 p.172

# データ
tone = [165, 130, 182, 178, 194, 206, 160, 122, 212, 165, 247, 195]
shinano = [180, 180, 235, 270, 240, 285, 164, 152]

続いて、有意水準 $${5\%}$$、両側検定でウェルチの検定を実行します。

# 関数利用
result = pop_mean_diff_welch_ttest(tone, shinano)
result

【実行結果】
検定統計量 $${T=-1.636}$$ は、棄却限界値（下側の閾値） $${-2.192}$$ より大きい（テキストによると棄却域に入らない）ので、有意水準 $${5\%}$$ で帰無仮説は棄却できません。

つまり、２つの河川のイワナの体長に差があるとは言えないです。

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

### 図示 p.173

## 設定と準備
# 統計関連の値
N1, N2 = len(tone), len(shinano)  # 標本サイズ
alpha = result['alpha']           # 有意水準
lower, upper = result['c_value']  # 棄却限界値
t_val = result['t_value']         # 検定統計量
df = result['df']                 # 自由度
t_dist = stats.t(df=df)           # 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'自由度{df:.4f}の$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, stats.t.pdf(x_val_upper, df=N - 1),
                 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();

【実行結果】

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

$${T}$$ が両方の棄却域に入っていないので、帰無仮説を棄却できないのです。

③ Python ライブラリで検定実行
scipy.stats、pingouin、statsmodels を利用します。
いずれも引数に２つサンプルデータと等分散を仮定しないことを渡しています。

🖲️scipy.stats

# scipy.stats利用 ※equal_var=False=Trueでウェルチの検定になる
stats.ttest_ind(tone, shinano, equal_var=False)

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

🖲️pingouin

# pingouin利用 ※correction=Trueでウェルチの検定になる
pg.ttest(x=tone, y=shinano, correction=True)

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

🖲️statsmodels

# statsmodels利用 ※usevar='unequal'でウェルチの検定になる
ttest_ind(tone, shinano, usevar='unequal')

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

理解度チェック ２つの母平均の差の検定 p.174

とある検査値の男女差に関する母平均の差の検定です。
等分散性を仮定し、有意水準 $${5\%}$$、両側検定で、コードで淡々と書きます。

標本データを設定します。
sample1 がグループ１：女性、sample2 がグループ２：男性です。

### ２つの母平均の差の検定（両側検定） 総コレステロール値 p.174

# データ
sample1 = [292, 351, 284, 278, 322, 295, 282, 317,
           305, 296, 267, 272, 343, 298, 275]
sample2 = [265, 272, 248, 276, 284, 258, 289, 307,
           284, 273, 301, 268, 293, 284, 318]

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

# 関数利用
result = pop_mean_diff_ttest(sample1, sample2)
result

【実行結果】
帰無仮説は棄却されます。

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

# scipy.stats利用
stats.ttest_ind(sample1, sample2, equal_var=True)

【実行結果】

pingouin で検定を実行します。

# pingouin利用
pg.ttest(x=sample1, y=sample2, correction=False)

【実行結果】

statsmodels で検定を実行します。

# statsmodels利用
ttest_ind(sample1, sample2, usevar='pooled')

【実行結果】

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

### 図示 p.175

## 設定と準備
# 統計関連の値
N1, N2 = len(sample1), len(sample2)  # 標本サイズ
alpha = result['alpha']              # 有意水準
lower, upper = result['c_value']     # 棄却限界値
t_val = result['t_value']            # 検定統計量
t_dist = stats.t(df=N1 + N2 - 2)     # 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'自由度({N1}+{N2}-2)の$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();

【実行結果】

ちょっと寄り道：Brunner-Munzel検定

母平均の差の検定をネットで調べていたときに出会ったのが Brunner-Munzel検定 です。
ノンパラメトリック検定ですので、母集団が正規分布である仮定は不要、等分散性の仮定も不要です。
実は５章でノンパラメトリック検定を実践するのですが、Brunner-Munzel検定を含んでいません。
「旬のうちに」お届けいたします（注：旬は私起点です）。

📣 特徴 📣

• 2標本に対して正規分布の仮定も、等分散の仮定のどちらも必要ないノンパラメトリック検定

• データの中央値に有意差があるか検出する

• 帰無仮説：各群から値を1つずつ取ったとき、両群で大きな値が得られる確率が等しい

• $${p}$$ 値が有意水準以下の場合に、対立仮説「2標本の間に差がある」が支持される

scipy 公式サイトのコードを辿ります。
両側検定です。

# Brunner-Munzel検定 scipy.stats公式のコードをなぞってみる
x1 = [1,2,1,1,1,1,1,1,1,1,2,4,1,1]
x2 = [3,3,4,3,1,2,3,1,1,5,4]
w, p_value = stats.brunnermunzel(x1, x2)
w, p_value

【実行結果】
統計量と $${p}$$ 値が出力されました。
有意水準 5% で帰無仮説は棄却される流れです。

ウェルチの検定と比べてみます。

# 参考：ウェルチの検定 scipy.stats利用 ※equal_var=False=Trueでウェルチの検定になる
stats.ttest_ind(x1, x2, equal_var=False)

【実行結果】
こちらの有意水準 5% で帰無仮説が棄却されます。

データのヒストグラムを見てみますと、正規分布ではなさそうです。
このデータの場合は、ウェルチの検定などの「正規分布を仮定」する検定を使うのは、避けたほうがよさそうです。

# ２つデータのヒストグラム
plt.hist(x1, alpha=0.5, label='$x_1$')
plt.hist(x2, alpha=0.5, label='$x_2$')
plt.legend();

【実行結果】

興味のある方は、次のWEBサイトで情報をGETしましょう。

Section 4.6 対応のある２つの母平均の差の検定

■ 前説
対応のある２つの母平均の差の検定は「同じ人」「同じ対象」から取得した２つの標本の母平均の差を対象とします。
２つの標本は「対応関係」があり、母平均の差は「１つの母集団」と考えられるので、「１標本の母平均の検定」（Section 4.2）を用います。

■ 母集団が従う確率分布
テキストは、対応のある２つの母平均の差の検定に関して、母集団が従う確率分布を明記していません。
ただし、対応のある２つの母平均の差の検定は「１標本の母平均の検定」（Section 4.2）を用いますので、おそらく母集団が正規分布に従うことを想定していると思われます。

■ 仮説
覆したい母平均の差 $${\mu_1 - \mu_2}$$ の仮説を次のように設定します。

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

◆ ◆ ◆

■ 対応のある２つの母平均の差の検定（片側検定）例題 p.176
p.176 例）の同一被験者に対する２種類の体温計の測定値を考えます。
１つめの体温計の方が高い体温を表示する傾向があるとのことです。
対応のある２つの母平均の差の検定で調べます。

仮説は次のとおりです。
・帰無仮説 $${H_0}$$ 「母平均の差 $${\mu_1 - \mu_2 = 0}$$」
・対立仮説 $${H_1}$$ 「母平均の差 $${\mu_1 - \mu_2 > 0}$$」
　⇒片側検定（greater）：１つ目の体温計の方が大きな値を表示するか

① 自作関数で検定実行
２つのデータを設定し、差を計算します。

### 対応のある２つの母平均の差の検定 ２つの体温計 p.176

# データ
sample1 = [37.1, 36.2, 36.6, 37.4, 36.8, 36.7, 36.9, 37.4, 36.6, 36.7]
sample2 = [36.8, 36.6, 36.5, 37.0, 36.0, 36.5, 36.6, 37.1, 36.4, 36.7]

# 2つの標本の差
diff = [round(x - y, 1) for x, y in zip(sample1, sample2)]
diff

【実行結果】
２つのデータの差は次のようになります。

続いて検定を実行します。

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

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

つまり、１つ目の体温計の方が温度高めに表示されると言えます。

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

### 図示 p.178

## 設定と準備
# 統計関連の値
N = len(diff)                     # 標本サイズ
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([lower, 0, upper])
plt.legend();

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

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

理解度チェック 対応のある２つの母平均の差の検定 p.179

ある人のダイエット前後の体重に関する母平均の差の検定で、ダイエットと体重増減の関係を調べます。
有意水準 $${5\%}$$、両側検定で、コードで淡々と書きます。

標本データを設定し、差を計算します。
sample1 がダイエット前、sample2 がダイエット後です。

### 対応のある２つの標本の母平均の差の検定 リンゴダイエット前後の体重 p.179

# データ
sample1 = [53.0, 50.2, 59.4, 61.9, 58.5, 56.4, 53.4]
sample2 = [51.2, 48.7, 53.5, 56.1, 52.4, 52.9, 53.3]

# 2つの標本の差
diff = [round(x - y, 1) for x, y in zip(sample1, sample2)]
diff

【実行結果の】
差を表示しました。
差が正なので、ダイエット後の体重のほうが小さいようです。

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

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

【実行結果】
帰無仮説は棄却されます。

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

# scipy.stats利用
stats.ttest_1samp(diff, popmean=0)

【実行結果】

pingouin で検定を実行します。

# pingouin利用
pg.ttest(x=diff, y=0)

【実行結果】

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

### 図示

## 設定と準備
# 統計関連の値
N = len(diff)                     # 標本サイズ
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();

【実行結果】

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

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