見出し画像

「入門はじめての統計解析」をPythonで写経 Vol.12 ~ 4章「はじめての統計的検定」⑤カイ二乗検定:適合度検定、独立性の検定

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

書籍の著者 石村貞夫 先生


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

書籍の図・表・計算を淡々とPython化する写経シリーズです。
この記事は4章の統計的検定テーマのうち 適合度検定独立性の検定 に取り組みます。

ChatGPTの活用も継続してまいります!

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

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

はじめに


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

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

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


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

4.10 適合度検定
4.11 独立性の検定

記事に用いるデータは、テキストに掲載されたデータそのものを引用しています。
データ件数の少ないものはコード上でデータを登録し、データ件数の多いものは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による「やさしいイントロ」をご堪能ください。


「アンケートの選択肢、どれも同じくらい人気のはずなんだけど…ほんとにそう?」
「仕事の種類によって、平日の過ごし方って違うのかな?」

こんなふうに、数字をじーっと見るというより、
“分け方に注目して”データを比べてみたくなることってありますよね。今回のテーマは、そんなときにぴったりの カイ二乗検定です。

この検定には、2つの仲間がいます。

🔸 適合度検定
→ データの出方が「このくらいのはず」と思っていた割合と合っているかどうかを調べます。
たとえば、選択肢がほんとうに均等に選ばれているのか?など、
期待と現実のズレに注目します。

🔸 独立性の検定
→ 「この2つの分け方は、おたがいに関係なくバラバラなのかな?」をチェックします。
たとえば、職業によって平日の過ごし方に偏りがあるか?みたいな感じです。
組み合わせに“かたより”があるかどうかを見ていく検定です。

分け方に注目することで、データの見え方がちょっと変わってきます。
「なんとなくそうかも?」が、「やっぱりそうだった!」になるかもしれません。
さあ、今回は“分けて比べる”データの楽しさに触れてみましょう!


Section 4.10 適合度検定

適合度検定は特定の母集団分布を仮定しないノンパラメトリック検定です。

■ 適合度検定の前に軽い例でウォーミングアップ
サイコロを60回投げて出た目の回数(「度数」といいます)が次の表のようになりました。

$$
\begin{array}{l:rrrrrr:r}
サイコロの目 & 1 & 2 & 3 & 4 & 5 & 6 & 合計\\
\hline
 \\
出た目の回数(観測度数) &  12 & 9 & 13 & 11 & 7 & 8 & 60\\
 \\
理論的な回数(期待度数) & 10 & 10 & 10 & 10 & 10 & 10 & 60\\
\end{array}
$$

サイコロの出た目の回数のように実際のデータ・標本採取したデータの度数を「観測度数」と呼びます。

ところで、サイコロの目は1から6まで均等に出るのが理論的です。
「理論どおりだったらこうなるはずだ」の度数を「期待度数」と呼びます。
理論どおりとは理論・仮説・想定する確率分布(理論分布)が当てはまります。
サイコロの場合は離散一様分布(どの目も均等)です。

観測度数と期待度数を棒グラフで可視化します。
期待度数は全て同じ値(一様・均等)ですが、観測度数は凸凹があります。

ここまでをまとめます。
サイコロの6つの目をカテゴリ変数として扱います。

適合度検定は「カテゴリ変数」(質的変数)の各カテゴリの度数(個数)に関する検定です。
具体的には「観測度数」が理論分布と「適合しているかどうか」(観測値が理論値に合致しているか)に関する検定です。

この例の適合度検定の結果は「有意水準 5% で帰無仮説は棄却されません」です。

サイコロのキャラクター:「いらすとや」さんより

■ 仮説
仮説を次のように設定します。

帰無仮説 $${H_0}$$:
 $${n}$$ 個のカテゴリの比率は $${p_1, p_2, \cdots, p_n}$$ である
対立仮説 $${H_1}$$:
 $${n}$$ 個のカテゴリの比率は $${p_1, p_2, \cdots, p_n}$$ とは異なる

テキストより引用

「理論分布と一致する/異なる」と読み替えても良さそうです。
なお、棄却域を求める際は分布の上側(片側)を見るようです。

◆ ◆ ◆

■ 適合度検定:検定統計量の公式 p.197
適合度検定における検定統計量の公式をテキストよりお借りします。

適合度検定の検定統計量

$$
T(f_i, N) = \cfrac{(f_1 - NP_1)^2}{NP_1} + \cfrac{(f_2 - NP_2)^2}{NP_2} + \cdots + \cfrac{(f_N - NP_N)^2}{NP_N}
$$

テキストの数式を引用

は自由度 $${n-1}$$ の $${\chi^2}$$ 分布に従います。
$${N}$$ は標本サイズ、$${n}$$ はカテゴリの個数、カテゴリ$${i}$$ に関して $${f_i}$$ は観測度数、$${p_i}$$ は比率($${\sum_{i=1}^n p_i = 1}$$)、$${Np_i}$$ は期待度数です。

式自体はややこしい印象かもですが、1つ1つの項は「観測度数から期待度数を引いた差の二乗を期待度数で割ってます」なのです!

公式に則って「適合度検定関数」を定義します。
確率計算には scipy.stats を利用します。

### 適合度検定関数 p.197 ※片側検定(上側)のみ

# 適合度検定のχ²検定統計量の算出関数
def chi2_stat_goodness_fit(x_list, ratio):
    # 期待度数Npの算出
    Nps = [sum(x_list) * r / sum(ratio) for r in ratio]
    # 期待度数<5を含む場合、アラートを出す
    text = '期待度数 < 5' if min(Nps) < 5 else ''
    # 戻り値:χ²検定統計量, 期待度数Nps, アラートtext
    return sum([(f - Np)**2 / Np for f, Np in zip(x_list, Nps)]), Nps, text

# 適合度検定関数 ※確率計算はscipy.stats利用
def goodness_fit_chi2test(x_list, ratio, alpha=0.05):
    # 標本サイズ
    n = len(x_list)
    # χ²検定統計量、期待度数、アラートの算出
    chi2_val, expected, text = chi2_stat_goodness_fit(x_list, ratio)
    # カイ二乗分布の設定
    chi2_dist = stats.chi2(df=n - 1)
    # 棄却限界値の算出(上側のみ)
    c_val = chi2_dist.ppf(q=1 - alpha)
    # p値の算出
    p_val = chi2_dist.sf(x=chi2_val)
    # 戻り値
    return {'chi2_value': chi2_val, 'c_value': c_val, 'alpha': alpha,
            'p_value': p_val, 'expected': expected, 'alert': text}

テストデータ:観測度数 $${[57, 33, 46, 14]}$$、理論的な比$${[4, 2, 3, 1]}$$ について、帰無仮説 $${H_0}$$ 「$${n}$$ 個のカテゴリの比率は $${p_1, \ldots, p_N}$$ である」、有意水準 5%($${\alpha=0.05}$$)、片側検定(greater)で適合度検定を実行します。

# テスト
obs, ratio = [57, 33, 46, 14], [4, 2, 3, 1]
goodness_fit_chi2test(obs, ratio)

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

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

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

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

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

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

  • expected:期待度数

  • alart:期待度数$${\leq 5}$$ のカテゴリを含むときアラートを出力

scipy.stats で検算します。

# scipy.statsで答え合わせ
expected = [sum(obs) * r / sum(ratio) for r in ratio]
stats.chisquare(obs, expected)

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

◆ ◆ ◆

■ 適合度検定(片側検定)例題 p.198
p.198 例)のアンケートで取得した血液型の観測度数と理論比率の適合度について、適合度検定で調べます。

仮説は次のとおりです。
血液型の理論分布を $${p_A=0.4, p_B=0.2, p_O=0.3, p_{AB}=0.1}$$ とすると
・帰無仮説 $${H_0}$$ 「観測した血液型の比は理論分布と一致する」
・対立仮説 $${H_1}$$ 「観測した血液型の比理論分布と一致しない」

① 自作関数で検定実行
観測度数(人数)と理論分布(比)のデータを設定します。

### 適合度検定 血液型アンケート p.198

# データ
data4 = pd.DataFrame({'人数': [57, 33, 46, 14], '比': [4, 2, 3, 1]},
                     index=['A型', 'B型', 'O型', 'AB型'])
data4.index.name = '血液型'
data4.T

【実行結果】

期待度数を算出します。

# 期待度数を算出

# data4をコピー
data4_ratio = data4.copy()
# 期待度数を算出
data4_ratio['期待度数'] = data4['人数'].sum() * data4['比'] / data4['比'].sum()
# 列名の絞り込み
data4_ratio = data4_ratio[['人数', '期待度数']]
# 結果の表示
data4_ratio.round(3).T

【実行結果】
観測度数(人数)と期待度数はよく似ています。

観測度数と期待度数を可視化して、分布の一致/ズレを直感的に確認します。
pandas データフレームの plot メソッドで描画します。

# 人数(観測度数)と期待度数の可視化
data4_ratio.plot.bar(rot=0, alpha=0.7);

【実行結果】
各血液型の観測度数(人数)と期待度数はよく似ています。

では有意水準 $${5\%}$$、片側検定で適合度検定を実行します。

### 適合度検定

# 関数利用
result = goodness_fit_chi2test(data4['人数'], data4['比'])
result

【実行結果】
検定統計量 $${0.5389}$$ は、棄却限界値$${7.8147}$$ より小さい(テキストによると棄却域に入らない)ので、有意水準 $${5\%}$$ で帰無仮説を棄却できません。
つまり、観測した血液型の比は 理論分布と異なるとは言えない です。

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

### 図示 p.199

## 設定と準備
# 統計関連の値
N = len(data4)                    # 標本サイズ
alpha = result['alpha']           # 有意水準
upper = 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_upper = np.linspace(upper, x_max, 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_upper, 0, chi2_dist.pdf(x_val_upper),
                 color=color, alpha=0.3, label=f'有意水準{alpha}の棄却域')
# 検定統計量の垂直線の描画
plt.axvline(chi2_val, color='tab:red', ls='--', label=f'検定統計量={chi2_val:.3f}')
# 修飾
plt.xticks([0, upper])
plt.legend();

【実行結果】

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

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

③ Python ライブラリで検定実行
ここでは、scipy.stats を利用します。
いずれも引数に2つのサンプルデータを渡しています。

🖲️scipy.stats

# scipy.stat利用
data4['期待度数'] = [data4['人数'].sum() * r / data4['比'].sum()
                    for r in data4['比']]
stats.chisquare(data4['人数'], data4['期待度数'])

【実行結果】
検定統計量と $${p}$$ 値を得ました。

理解度チェック 適合度検定 p.200

サヤエンドウの4種類の「色・形状」の比が理論比かどうかを、適合度検定で確かめます。

有意水準 $${5\%}$$、片側検定で、コードを淡々と書きます。
観測度数と理論比を設定します。

### 適合度検定 サヤエンドウ豆の種類の比 p.200

data5 = pd.DataFrame({'観測度数': [105, 38, 41, 16], '比': [9, 3, 3, 1]},
                     index=['黄色丸型', '黄色角型', '緑色丸型', '緑色角型'])
data5.index.name = 'サヤエンドウ豆'
data5.T

【実行結果】
「比」が理論比です。

観測度数と期待度数はよく似ています。

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

### 適合度検定

# 関数利用
result = goodness_fit_chi2test(data5['観測度数'], data5['比'])
result

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

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

# scipy.stat利用
data5['期待度数'] = [data5['観測度数'].sum() * r / data5['比'].sum()
                    for r in data5['比']]
stats.chisquare(data5['観測度数'], data5['期待度数'])

【実行結果】

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

### 図示 p.201

## 設定と準備
# 統計関連の値
N = len(data4)                    # 標本サイズ
alpha = result['alpha']           # 有意水準
upper = 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_upper = np.linspace(upper, x_max, 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_upper, 0, chi2_dist.pdf(x_val_upper),
                 color=color, alpha=0.3, label=f'有意水準{alpha}の棄却域')
# 検定統計量の垂直線の描画
plt.axvline(chi2_val, color='tab:red', ls='--', label=f'検定統計量={chi2_val:.3f}')
# 修飾
plt.xticks([0, upper])
plt.legend();

【実行結果】
上側(右側)の青い棄却域に検定統計量は乗っていません。
帰無仮説を棄却できないのです。

Section 4.11 独立性の検定

■ 独立性の検定の前に軽い例でウォーミングアップ
2つのカテゴリ変数が「独立があるか」または「関連性があるか」を調べるのが独立性の検定です。
2つのカテゴリ変数の「クロス集計表」の値が「観測度数」となり、一定の計算方法で「期待度数」を算出すると、あとは適合度検定とよく似た検定統計量を計算できます。

クロス集計表の例です。
訪日外国人に聞いた「食べてみたい日本食」の集計(仮想データ)です。
「観測度数」の表がアンケート結果の集計値、「期待度数」は観測度数から一定の方法で計算した結果です。
この内容で独立性の検定を行うと、有意水準 5% で有意になります。

まとめると…
独立性の検定は2つのカテゴリ変数(質的変数)のクロス集計表に関する検定です。
具体的には「観測度数」と「期待度数」から帰無仮説「2つのカテゴリ変数独立している(関連がない)」、対立仮説「独立ではない(関連がある)」を検定するものです。

◆ ◆ ◆

■ 仮説
仮説を次のように設定します。

帰無仮説 $${H_0}$$:
 2つの属性 $${A}$$ と $${B}$$ は独立である
対立仮説 $${H_1}$$:
 2つの属性 $${A}$$ と $${B}$$ の間には関連がある

テキストより引用

なお、棄却域を求める際は分布の上側(片側)を見るようです。

◆ ◆ ◆

■ 独立性の検定:検定統計量の公式 p.203
独立性の検定における検定統計量の公式をテキストよりお借りします。
カテゴリ $${A_m, B_n}$$ のクロス集計表を次のように書くとき:

$$
\begin{array}{c|cccc|c}
& B_1 & B_2 & \cdots & B_n & 合計 \\
\hline
A_1 & f_{11} & f_{12} & \cdots & f_{1n} & f_{1B} \\
A_2 & f_{21} & f_{22} & \cdots & f_{2n} & f_{2B} \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
A_m & f_{m1} & f_{m2} & \cdots & f_{mn} & f_{mB} \\
\hline
合計 & f_{A1} & f_{A2} & \cdots & f_{An} & N \\
\end{array}
$$

独立性の検定統計量

$$
T(f_{ij}, N) = \sum_{i=1}^m \sum_{j=1}^n \cfrac{(Nf_{ij} - f_{iB}f_{Aj})^2}{Nf_{iB}f_{Aj}}
$$

テキストの数式を引用

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

公式に則って「独立性の検定関数」を定義します。
確率計算には scipy.stats を利用します。

### 独立性の検定関数 p.202 ※片側検定(上側)のみ

# 独立性の検定のχ²検定統計量の算出関数
def chi2_stat_independence(crosstab):
    # crosstabを変数fに代入
    f = crosstab
    # 行数の算出
    m = len(f)
    # 列数の算出
    n = len(f[0])
    # 横計f_*Bの算出
    fB = [sum(x) for x in f]
    # 縦計f_A*の算出
    fA = [sum(x[j] for x in f) for j in range(n)]
    # 総合計Nの算出
    N = sum(fB)
    # 戻り値:χ²検定統計量, 行数m, 列数n
    return sum([(N*f[i][j] - fB[i]*fA[j])**2 / (N*fB[i]*fA[j])
                for i in range(m) for j in range(n)]), m, n

# 独立性の検定関数 ※確率計算はscipy.stats利用
def independence_chi2test(crosstab, alpha=0.05):
    # χ²検定統計量、行数m、列数nの算出
    chi2_val, m, n = chi2_stat_independence(crosstab)
    # 自由度の算出
    df = (m - 1) * (n - 1)
    # カイ二乗分布の設定
    chi2_dist = stats.chi2(df=df)
    # 棄却限界値の算出(上側のみ)
    c_val = chi2_dist.ppf(q=1 - alpha)
    # p値の算出
    p_val = chi2_dist.sf(x=chi2_val)
    # 戻り値
    return {'chi2_value': chi2_val, 'c_value': c_val, 'alpha': alpha,
            'p_value': p_val, 'df': df}

テストデータ $${[105, 38, 41, 16]}$$、$${[89, 24, 75, 12]}$$ について、帰無仮説 $${H_0}$$ 「2つのカテゴリは独立である」、有意水準 5%($${\alpha=0.05}$$)、上側検定(greater)で母平均の検定を実行します。

# テスト
crosstab = [[105, 38, 41, 16], [89, 24, 75, 12]]
independence_chi2test(crosstab)

【実行結果】
「有意水準 5% で有意となり、帰無仮説は棄却される」です。

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

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

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

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

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

  • df:自由度

scipy.stats で検算します。
前半には 出力結果の表示関数を書いています。

# scipy.statsで答え合わせ

# 戻り値をprintする関数の定義
def print_res(res):
    print('statistic:', res.statistic)
    print('pvalue   :', res.pvalue)
    print('dof      :', res.dof)
    print('expected_freq:')
    print(res.expected_freq)

# 独立性の検定の実行
res = stats.chi2_contingency(crosstab)
print_res(res)

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

◆ ◆ ◆

■ 独立性の検定(片側検定)例題 p.198
p.204 例)の国別の血液型の独立性について、独立性の検定で調べます。

仮説は次のとおりです。
・帰無仮説 $${H_0}$$ 「国と血液型は独立である」
・対立仮説 $${H_1}$$ 「国と血液型の間には関連がある」

① 自作関数で検定実行
観測度数(人数)と理論分布(比)のデータを設定します。

### 独立性の検定 血液型 p.204

# データ
data6 = pd.DataFrame({'日本人': [57, 33, 46, 14], 'フランス人': [89, 24, 75, 12]},
                     index=['A型', 'B型', 'O型', 'AB型']).T
data6.index.name = '血液型'
display(data6)

【実行結果】

観測度数の比率を算出して可視化し、国別の分布を直感的に確認します。

# 比率の作成
data6_ratio = data6.T.copy()
data6_ratio = data6_ratio / data6_ratio.sum()
display(data6_ratio.T.round(3))
# 可視化
data6_ratio.plot.bar(rot=0, alpha=0.7);

【実行結果】
両国の分布は割と似ている印象です。
果たして、独立である(関連がない)のでしょうか?

では有意水準 $${5\%}$$、片側検定で適合度検定を実行します。

### 独立性の検定

# 関数利用
result = independence_chi2test(data6.values)
result

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

つまり、国と血液型の間には 関連があると言える、です。

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

### 図示 p.205

## 設定と準備
# 統計関連の値
alpha = result['alpha']           # 有意水準
upper = result['c_value']         # 棄却限界値
chi2_val = result['chi2_value']   # 検定統計量
df = result['df']                 # 自由度
chi2_dist = stats.chi2(df=df)     # カイ二乗分布
# グラフ描画用設定
color = 'tab:blue'                            # 基本の色
x_min, x_max = 0, 25                          # 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))
# カイ二乗分布の確率密度関数の描画
plt.plot(x_val, chi2_dist.pdf(x_val), color=color, 
         label=f'自由度{df}の$\chi^2$分布')
# 上側棄却域の塗りつぶし
plt.fill_between(x_val_upper, 0, chi2_dist.pdf(x_val_upper),
                 color=color, alpha=0.3, label=f'有意水準{alpha}の棄却域')
# 検定統計量の垂直線の描画
plt.axvline(chi2_val, color='tab:red', ls='--', label=f'検定統計量={chi2_val:.3f}')
# 修飾
plt.xticks([0, upper])
plt.legend();

【実行結果】

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

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

③ Python ライブラリで検定実行
ここでは、scipy.stats を利用します。
いずれも引数に2つのサンプルデータを渡しています。

🖲️scipy.stats

# scipy.stats利用
res = stats.chi2_contingency(data6)
print_res(res)

【実行結果】
検定統計量、$${p}$$ 値、自由度、期待度数を得ました。

理解度チェック 独立性の検定 p.206

武道の段位と技の独立性/適合性を独立性の検定で確かめます。

有意水準 $${5\%}$$、片側検定で、コードを淡々と書きます。
観測度数を設定します。

### 独立性の検定 段位と打突部位 p.206

# データ
data7 = pd.DataFrame({'七段・八段': [17, 19, 1, 2], '三段・四段': [21, 12, 9, 5]},
                     index=['面', '小手', '胴', '突き']).T
data7.index.name = '段位'
data7

【実行結果】

観測度数の比率を算出して可視化します。

# 比率の作成
data7_ratio = data7.T.copy()
data7_ratio = data7_ratio / data7_ratio.sum()
display(data7_ratio.T.round(3))
# 可視化
data7_ratio.plot.bar(rot=0, alpha=0.7);

【実行結果】
小手・胴・突きに段位の特徴がありそうです。

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

### 独立性の検定

# 関数利用
result = independence_chi2test(data7.values)
result

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

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

# scipy.stats利用
res = stats.chi2_contingency(data7)
print_res(res)

【実行結果】

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

### 図示 p.207

## 設定と準備
# 統計関連の値
alpha = result['alpha']           # 有意水準
upper = result['c_value']         # 棄却限界値
chi2_val = result['chi2_value']   # 検定統計量
df = result['df']                 # 自由度
chi2_dist = stats.chi2(df=df)     # カイ二乗分布
# グラフ描画用設定
color = 'tab:blue'                            # 基本の色
x_min, x_max = 0, 25                          # 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))
# カイ二乗分布の確率密度関数の描画
plt.plot(x_val, chi2_dist.pdf(x_val), color=color, 
         label=f'自由度{df}の$\chi^2$分布')
# 上側棄却域の塗りつぶし
plt.fill_between(x_val_upper, 0, chi2_dist.pdf(x_val_upper),
                 color=color, alpha=0.3, label=f'有意水準{alpha}の棄却域')
# 検定統計量の垂直線の描画
plt.axvline(chi2_val, color='tab:red', ls='--', label=f'検定統計量={chi2_val:.3f}')
# 修飾
plt.xticks([0, upper])
plt.legend();

【実行結果】
上側(右側)の青い棄却域に検定統計量が位置しています。
帰無仮説は棄却されるのです。


記事の最後を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 応援ありがとうございます。これからもがんばって記事を作成します!

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