「入門はじめての統計解析」をPythonで写経 Vol.13 ～ ４章「はじめての統計的検定」⑥外れ値の検定、正規性の検定、歪度と尖度の検定

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

書籍の著者　石村貞夫 先生

この記事は、書籍「入門はじめての統計解析」４章「はじめての統計的検定」の 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による「やさしいイントロ」をご堪能ください。

「この数字、なんだかひとつだけ浮いてる気がする…」
「このデータ、きれいな山型に見えるけど、ほんとうに正規分布なのかな？」

データを集めて分析するとき、“値そのもの”だけでなく、“並び方やかたち”に目を向けたくなることがあります。
今回のテーマは、そんなときに役立つ 外れ値や分布のかたちをチェックする検定たちです。

✨ たとえば、こんな検定があります：

• グラブス検定・スミルノフ検定：目立って飛び出したデータが、本当に“外れ値”なのかを判断します

• コルモゴロフ・スミルノフ検定：観測されたデータが、正規分布のような理論分布と似ているかを確かめます

• 正規Q-Qプロット：数字の並び方が、正規分布らしいパターンに沿っているかを直感的に見られます

• 尖度・歪度の検定：分布の“とんがり具合”や“ゆがみ”が、どの程度なのかを測定します

これまで学んできた統計的検定では、「2つの差があるか？」「関係があるか？」などをデータから確かめてきました。
今回はその視点を少し広げて、「このデータ、そもそもちゃんと使っていいかな？」という“データの前提”を見直す検定を学んでいきます。
統計的検定の締めくくりとして、データの“性格”を見抜く目を育てるラストステージ。
これからの分析にもきっと役立つ、そんなチェックの方法を一緒に身につけていきましょう！

Section 4.12 外れ値の検定

■ ウォーミングアップ
テキストによると外れ値は「非常に飛び離れた値のデータ」です。
外れ値を可視化して直感的に捉えてみましょう。

# 例
url = 'https://vincentarelbundock.github.io/Rdatasets/csv/carData/Davis.csv'
ex_df = pd.read_csv(url)
ex_df.plot.scatter(x='weight', y='height', s=70, ec='white', alpha=0.9,
                   xlabel='体重', ylabel='身長', title='Davisデータセット');

【実行結果】
右下に１点、非常に飛び離れた「外れ値」が存在します。

箱ひげ図で外れ値の様子を確かめてみましょう。

# weightの箱ひげ図
plt.boxplot(ex_df['weight'], orientation='horizontal', tick_labels=['体重']);

【実行結果】
最右点が外れ値候補です。

外れ値の検定は、データ分析の前処理段階で「外れ値」を取り除くかどうかを検討するときに活用できます。

ところで…
Davisデータセットを見ると、いつも雷神の屏風を思い出します。
（脱線でした）

■ 外れ値の検定：検定統計量の公式 p.208
テキストは「グラブス・スミルノフの外れ値の検定」を紹介しています。
母集団が従う確率分布に「正規分布」を想定しているようです。
検定統計量の公式をテキストよりお借りします。

正規母集団から抽出した標本 $${\{ x_1, x_2, \ldots, x_N\}}$$ において、データ $${x_k}$$ だけが非常に飛び離れているとき、この外れ値の検定統計量 $${T(x_k)}$$ は

となり、棄却域はグラブス・スミルノフの数表より求めます（テキスト p.292 の数表 5 参照）。
$${\bar{x}}$$ は標本平均、$${s}$$ は標本分散 $${s^2}$$ の標準偏差です。

テキストの公式は「外れ値候補データに対する標準化の結果を絶対値にしている」です！

仮説は次のとおりです。
・帰無仮説 $${H_0}$$ 「$${x_k}$$は外れ値ではない」
・対立仮説 $${H_1}$$ 「$${x_k}$$は外れ値である」

検定の都度、「数表」をめくって棄却域を取得するのは大変だと感じ…
次のWebサイトを参考にして自作関数の制作に取り込みました。
ありがとうございます！

公式やWebサイトの情報に基づいて「グラブス・スミルノフ検定関数」を定義します。
確率計算には scipy.stats を利用します。

### グラブス・スミルノフの外れ値の検定関数 p.210 ★p値,t分布,棄却域は自信なし

# 参考サイト: https://ja.wikipedia.org/wiki/%E5%A4%96%E3%82%8C%E5%80%A4
#            https://lactivator.net/2020/09/16/abnormal_value/

# 上側100α％点の近似的な有意点の算出関数
def t_significant_point(N, alpha, alternative='two-sided'):
    # 両側検定のとき、αを1/2する
    coef = 2 if alternative=='two-sided' else 1
    # 有意点tに用いるt分布の％点の算出
    t = stats.t.isf(q=alpha/(coef * N), df=N - 2)
    # 戻り値: 有意点t
    return (N - 1) * t / (N*(N - 2) + N*t**2)**(1/2)

# スミルノフ・グラブス検定のT検定統計量の算出関数
def t_stats_smirnov_grubbs(x_list):
    # 標本サイズ
    N = len(x_list)
    # 標本平均
    x_bar = sum(x_list) / N
    # 標本分散
    s2 = sum([(x - x_bar)**2 for x in x_list]) / (N - 1)
    # T検定統計量の算出
    Ts = [abs(x - x_bar) / s2**(1/2) for x in x_list]
    # T検定統計量が最大となる値を外れ値outlierとする
    T_max = max(Ts)
    T_max_idx = Ts.index(T_max)
    outlier = x_list[T_max_idx]
    # 外れ値が最小値か最大値かを判定
    min_max = 'min' if outlier == min(x_list) else 'max'
    # 戻り値: 外れ値、外れ値のインデックス検定統計量
    return outlier, T_max_idx, T_max, min_max

# スミルノフ・グラブス検定関数
def smirnov_grubbs_test(x_list, alpha=0.05, alternative='two-sided'):
    # 標本サイズ
    N = len(x_list)
    # 有意点c_valの算出
    c_val = t_significant_point(N, alpha, alternative)
    # T検定統計量の算出
    outlier, idx, T, min_max = t_stats_smirnov_grubbs(x_list)
    # p値の算出
    coef = 2 if alternative=='two-sided' else 1
    p_val = stats.t.sf(x=T, df=N - 2) * coef
    # 戻り値: 外れ値、外れ値のインデックス、T検定統計量、有意点
    return {'outlier': outlier, 'outlier_index': idx, 'min_max': min_max,
            'T_stat': T, 'c_value': c_val, 'p_value': p_val, 'alpha': alpha,
            'alternative': alternative}

テストデータ  $${[3.4, 3.5, 3.3, 2.2, 3.3, 3.4, 3.6, 3.2]}$$ について、帰無仮説 $${H_0}$$ 「検定統計量が最も大きなデータ点は外れ値ではない」、有意水準 5%（$${\alpha=0.05}$$）、片側検定（less）で外れ値の検定を実行します。

# テスト
x_list = [3.4, 3.5, 3.3, 2.2, 3.3, 3.4, 3.6, 3.2]
smirnov_grubbs_test(x_list, alternative='less')

【実行結果】
「インデックス $${3}$$ の データ点（値 $${2.2}$$）は、有意水準 5% で有意であり、帰無仮説は棄却され、外れ値と言える」です。

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

• outlier：検定統計量 $${T}$$ 最大値のデータ点の値

• outlier_index：外れ値候補 outlier のインデックス

• min_max：外れ値候補が最小値の場合 min、最大値の場合 max

• T_stat：検定統計量 $${T}$$

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

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

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

• alternative：両側検定 two-sided、片側検定（下側） less、片側検定（上側） greater

◆ ◆ ◆

■ グラブス・スミルノフの外れ値の検定（片側検定）例題 p.210
p.210 例）の血中物質の標本に含まれる小さな値 $${x_4}$$ が外れ値かどうか、外れ値の検定で調べます。

仮説は次のとおりです。
・帰無仮説 $${H_0}$$ 「$${x_4=2.2}$$ は外れ値ではない」
・対立仮説 $${H_1}$$ 「$${x_4=2.2}$$ は外れ値である」

① 自作関数で検定実行
データを設定して箱ひげ図で可視化します。

### グラブス・スミルノフの外れ値の検定 ヘモグロビン量 p.210

# データ
data8 = [3.4, 3.5, 3.3, 2.2, 3.3, 3.4, 3.6, 3.2]

# 箱ひげ図
plt.boxplot(data8, orientation='horizontal');

【実行結果】
めっちゃ外れ値な感じがいたします。

有意水準 $${5\%}$$、片側検定で外れ値の検定を実行します。

# 関数利用
result = smirnov_grubbs_test(data8, alternative='less')
result

【実行結果】
インデックス $${3}$$（Python的に「４番目」）のデータ点 $${x_k=2.2}$$ の検定統計量 $${2.3724}$$ は、棄却限界値$${2.0317}$$ より大きい（テキストによると棄却域に入る）ので、有意水準 $${5\%}$$ で帰無仮説は棄却されます。

つまり、値 $${2.2}$$ のデータ点は 外れ値であると言える、です。

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

### 図示 p.210

## 設定と準備
# 統計関連の値
N = len(data8)                    # 標本サイズ
alpha = result['alpha']           # 有意水準
lower = -result['c_value']        # 棄却限界値_下端
coef = -1 if result['min_max']=='min' else 1   # 最小値⇒T_statsをマイナスにする
t_val = result['T_stat'] * coef   # 検定統計量
t_dist = stats.t(df=N - 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'自由度({N}-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.title(f"外れ値 {result['outlier']} ({result['min_max']})")
plt.xticks([lower, 0])
plt.legend();

【実行結果】
ベル型の曲線は自由度 $${n-2=6}$$ の $${t}$$ 分布の確率密度関数です。
右側の青い領域が有意水準 $${5\%}$$ の棄却域です。
赤い点線が検定統計量 $${T}$$ です。

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

理解度チェック グラブス・スミルノフの外れ値の検定 p.211

物質の融点側転データに含まれる「大きく飛び離れたデータ」が外れ値かどうか、外れ値の検定で確かめます。

有意水準 $${5\%}$$、片側検定で、コードを淡々と書きます。
データを設定して箱ひげ図で可視化します。

### グラブス・スミルノフの外れ値の検定 マンガンの融点 p.211

# データ
data9 = [1276, 1265, 1302, 1283, 1271, 1269, 1285]

# 箱ひげ図
plt.boxplot(data9, orientation='horizontal');

【実行結果】
上側（右側）に長いヒゲが見られますが、箱ひげ図における外れ値（○で示される点）には該当しないようです。

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

# 関数利用
result = smirnov_grubbs_test(data9, alternative='greater')
result

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

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

### 図示 p.210

## 設定と準備
# 統計関連の値
N = len(data9)                    # 標本サイズ
alpha = result['alpha']           # 有意水準
upper = result['c_value']         # 棄却限界値_上端
coef = -1 if result['min_max']=='min' else 1  # 最小値⇒T_statsをマイナスにする
t_val = result['T_stat'] * coef   # 検定統計量
t_dist = stats.t(df=N - 2)        # 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}-2)の$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.title(f"外れ値 {result['outlier']} ({result['min_max']})")
plt.xticks([0, upper])
plt.legend();

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

Section 4.13 正規性の検定

■ ウォーミングアップ
正規性の検定はデータが「正規分布」に従っているかどうかの検定です。
データ分析のいろんな段階で「正規性の検定」が活躍すると思います！

• 母平均の検定など、母集団分布に正規分布を想定した検定手法の利用前

  • 標本が正規分布に従っていることが検定手法採用の適否に

• 線形回帰分析で推定したパラメータを用いたモデルの評価

  • 残差が正規分布に従っているかどうかはモデルの適否に

テキストは「コルモゴロフ・スミルノフの検定」に加えて、正規性を可視化できる「正規Q-Qプロット」を紹介しています。

◆ ◆ ◆

■ コルモゴロフ・スミルノフの検定
テキストに数式の紹介は無く、統計用ソフト SPSS の出力結果を紹介しています。
仮説は次のとおりです。

• 帰無仮説：データの分布が仮定した確率分布に一致する

• 対立仮説：データの分布が仮定した確率分布に一致しない

実はこの検定は「さまざまな確率分布に対応」しています。
正規性の検定の場合には、仮定した確率分布を「正規分布」とし、検定を行います。

■ 正規Q-Qプロット
データが正規分布に従う場合、ある変換を行ったデータ点が直線的に並ぶことを可視化するチャートです。
変換内容は統計Webなどの参考サイトでぜひご確認ください！

■ 正規性の検定と可視化
仮想データを作成して正規性の検定と可視化を試してみましょう。

① 仮想データの作成
自由度 $${3}$$ の $${t}$$ 分布に従う乱数を 100 個生成します。

### コルモゴロフ・スミルノフの検定のお試し

# データの生成：標準正規分布乱数
rng = np.random.default_rng(seed=0)
data10 = rng.standard_t(df=3, size=100)

# ヒストグラムの描画
fig, ax = plt.subplots()
twinx = ax.twinx()
sns.histplot(x=data10, kde=True, ec='white', alpha=0.7, ax=ax)
twinx.plot(x_val:=np.linspace(-4.5, 4.5, 1001), stats.norm.pdf(x_val), 
           color='tab:red', lw=2, ls='--', label='標準正規分布')
twinx.legend();

【実行結果】
仮想データのヒストグラムとKDE、参考値としての標準正規分布の確率密度関数（赤い点線）です。
データは標準正規分布に近いような、そうで無いような、どちらでしょう？

② コルモゴロフ・スミルノフの検定の実行
scipy.stats の kstest() を実行します。
引数は、データ、標準正規分布の累積分布関数、です。

# コルモゴロフ・スミルノフの検定 scipy.stats利用
# 帰無仮説：データの分布が仮定した確率分布に一致する
# 対立仮説：データの分布が仮定した確率分布に一致しない
# ２つ目の引数 cdf にstats.norm.cdfを与えることで、標準正規分布を指定している

res = stats.kstest(data10, stats.norm.cdf)
print('検定統計量:', res.statistic)
print('p値　　　 :', res.pvalue)

【実行結果】

$${p}$$ 値 $${0.3350}$$ は有意水準 $${5\%}$$ で有意とは言えず、帰無仮説を棄却できません。
つまり、仮想データは標準正規分布と一致しないとは言えない、です。
仮想データは標準正規分布に従うだろう、としておきましょう。

③ シャピロ・ウィルク検定の実行
シャピロ・ウィルク検定も正規性の検定の一種です。
仮説は次のとおりです。

• 帰無仮説：データが正規分布に従う

• 対立仮説：データが正規分布に従わない

scipy.stats の shapiro() を実行しましょう。
引数はデータです。

# シャピロ・ウィルク検定 scipy.stats利用
# 帰無仮説：データが正規分布に従う、対立仮説：データが正規分布に従わない
res = stats.shapiro(data10)
print('検定統計量:', res.statistic)
print('p値　　　 :', res.pvalue)

【実行結果】

$${p}$$ 値は $${0.0239}$$ は有意水準 $${5\%}$$ で有意と言えるので、帰無仮説は棄却されます。
つまり、仮想データが正規分布に従わないと言える、です。
仮想データは標準正規分布に従わないだろう、としておきましょう。

コルモゴロフ・スミルノフの検定とシャピロ・ウィルク検定の結果が真逆になってしまいました。。。どうしよう。。。

そうだ、可視化してみよう！

③ 正規Q-Qプロットの描画

# 正規Q-Qプロット scipy.stats利用
stats.probplot(data10, dist='norm', plot=plt)
plt.show();

【実行結果】
データが正規分布に従う場合、青いデータ点は赤い直線上に乗ります。
この図では、データの両端において直線に乗らない状況となっております。
「仮想データは正規分布に従わないんでしょうね」ですかねぇ…

正規分布に従うデータの正規Q-Qプロットが気になりますね！
直線上にデータ点が並ぶ正規Q-Qプロットの例を見てみましょう。
100 個の標準正規分布乱数で確かめます。

# 標準正規分布乱数の正規Q-Qプロット

## 標準正規分布に従うデータの作成
# 乱数生成器の設定
rng = np.random.default_rng(seed=7)
# 標準正規分布乱数 100 個生成
std_norm_rvs = rng.normal(size=100)

## 可視化
# 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
# 標準正規分布乱数データのヒストグラム＋KDEの描画
sns.histplot(x=std_norm_rvs, kde=True, ec='white', alpha=0.7, ax=ax1);
# 標準正規分布乱数データの正規Q-Qプロットの描画
stats.probplot(std_norm_rvs, dist='norm', plot=ax2)
plt.show();

【実行結果】
青いデータ点は直線上に乗っかっています！

Section 4.14 歪度と尖度の検定

■ ウォーミングアップ
テキストの言葉をお借りしますと、歪度と尖度の検定は
「正規母集団のチェックに利用します」

正規分布の場合に
・歪度＝０
・尖度＝０（または３）
となることを利用した検定です。

歪度・尖度の計算はテキストの p.30 ～ 31 や次の記事をご覧ください。

◆ ◆ ◆

■ 歪度と尖度の検定
テキストは「信頼区間を活用した検定方法」を紹介しています。

• 信頼区間に０を含む場合、帰無仮説「歪度（または尖度）＝０」を棄却できない

• 棄却できない
  ⇒ データは正規分布の歪度（または尖度）だろう
  ⇒ データは正規分布に従うだろう

という流れです。

■ 歪度・尖度の信頼区間の算出
信頼区間の算出に必要な要素は次の３つです。

1. 歪度・尖度の値（統計量）

2. $${\alpha}$$（$${0.05}$$ など）の $${z}$$ 値

3. 歪度・尖度の標準誤差

信頼区間は「統計値 ± $${z}$$ 値 × 標準誤差」で求められます。

2. の $${z}$$ 値の計算方法は既に学んでおり、計算できます。
1. の計算は、テキストの p.30 ～ 31 の公式を使って計算できます。
しかし、3.「標準誤差」の計算は、テキストに数式・公式が掲載されておらず、どのように計算してよいのやら…

Webサイトを彷徨い（30分くらい）、尖度・歪度の標準誤差の計算方法を探しましたが、妥当そうなものは見つかりませんでした。

【参考】尖度・歪度の標準誤差の計算方法に言及した記事
１番目の記事のおかげで２番目の記事に到達できました。
ありがとうございます！
２番目の記事に近似値の計算式がありますが、かなり粗い感じです。

ここは…
ChatGPTに頼るしかありません！
ChatGPTが教えてくれた「歪度・尖度の標準誤差の公式」を掲載いたします。
ただし裏取りできていないので、公式の適否は分かりかねます。

■ 歪度と尖度の検定・信頼区間編の実装
仮想データを作成し、ChatGPTが書いてくれた歪度・尖度の信頼区間算出関数を使って信頼区間を求めます。
標準誤差の計算は上述の式で行っています。

① 仮想データの作成
平均 $${\mu=0}$$、分散 $${\sigma^2=1.5^2}$$ の正規分布に従う乱数を 100 個生成します。

### 歪度と尖度の検定 p.213

# データの生成：平均0, 標準偏差1.5の正規分布乱数
rng = np.random.default_rng(seed=5)
data11 = rng.normal(size=100, loc=0, scale=1.5)

# ヒストグラムの描画
fig, ax = plt.subplots()
twinx = ax.twinx()
sns.histplot(x=data11, kde=True, ec='white', alpha=0.7, ax=ax)
twinx.plot(x_val:=np.linspace(-4.5, 4.5, 1001),
           stats.norm.pdf(x_val, loc=0, scale=1.5), 
           color='tab:red', lw=2, ls='--', label='正規分布 $N(0, 1.5)$')
twinx.legend();

【実行結果】
仮想データのヒストグラムとKDE、参考値としての正規分布の確率密度関数（赤い点線）です。
仮想データはいい感じに崩れていて、正規分布に近いような、そうで無いような、どちらでしょう？

② 信頼区間を用いて歪度と尖度の検定を行う
ChatGPT秘伝の公式に基づいて関数化したコードです。
引数は、データ、$${\alpha}$$ です。
$${95\%}$$ 信頼区間の場合、$${\alpha=0.05}$$ です。

### 歪度skewnessと尖度kurtosisの信頼区間算出関数 by ChatGPT

def skew_kurtosis_confint(data, alpha=0.05):

    # 設定と準備
    n = len(data)                      # 標本サイズ
    z = stats.norm.ppf(1 - alpha / 2)  # 100x(1-α)%点のzスコア

    # 歪度と尖度の算出
    skewness = stats.skew(data, bias=False)
    kurt = stats.kurtosis(data, fisher=True, bias=False) # 正規分布=0

    # 歪度の標準誤差と信頼区間の算出
    se_skew = np.sqrt(6 * n * (n - 1) / ((n - 2) * (n + 1) * (n + 3)))
    skew_ci = (skewness - z * se_skew, skewness + z * se_skew)

    # 尖度の標準誤差と信頼区間の算出
    se_kurt = np.sqrt(24 * n * (n - 1)**2 / ((n - 3)*(n - 2)*(n + 3)*(n + 5)))
    kurt_ci = (kurt - z * se_kurt, kurt + z * se_kurt)

    return {
        'skewness': skewness,
        'skew_CI': skew_ci,
        'kurtosis': kurt,
        'kurtosis_CI': kurt_ci,
        'SE_skew': se_skew,
        'SE_kurt': se_kurt
    }

仮想データについて、歪度・尖度の $${95\%}$$ 信頼区間を求めましょう。

# 歪度・尖度の信頼区間の算出
# 帰無仮説：歪度・尖度は０である、対立仮説：歪度・尖度は０ではない
# 信頼区間に０が含まれる場合、帰無仮説を棄却できない

skew_kurtosis_confint(data11)

【実行結果】
歪度の信頼区間は skew_CI、尖度の信頼区間は kurtsis_CI です。

歪度の $${95\%}$$ 信頼区間は０を含んでいます。
尖度の $${95\%}$$ 信頼区間も同様に０を含んでいます。
それゆえ、有意水準 $${5\%}$$ で 帰無仮説を棄却できず、データの歪度と尖度が正規分布と異なるとは言えない、です。
データは正規分布だろう、という感じです。

◆ ◆ ◆

■ 歪度と尖度の検定・scipy.stats 編の実装
実は…今更…申し上げにくいのですが…
scipy.stats に歪度と尖度の検定関数が実装されています！
こちらは検定統計量と $$p$$ 値を計算してくれます。

歪度の検定を行います。

# 歪度の検定 statisticはzスコア
# 帰無仮説：標本の母集団の歪度が正規分布の歪度である

res = stats.skewtest(data11)
print('検定統計量:', res.statistic)
print('p値　　　 :', res.pvalue)

【結果】

$${p}$$ 値 $${=0.1666}$$ は有意水準 $${5\%}$$（$${0.05}$$）より大きいです。
有意水準 $${5\%}$$ で 帰無仮説を棄却できず、データの歪度が正規分布と異なるとは言えない、です。
歪度的にデータは正規分布だろう、という感じです。

続いて尖度の検定を行います。

# 尖度の検定 statisticはzスコア
# 帰無仮説：標本の母集団の尖度が正規分布の尖度である

res = stats.kurtosistest(data11)
print('検定統計量:', res.statistic)
print('p値　　　 :', res.pvalue)

【結果】

$${p}$$ 値 $${=0.8725}$$ は有意水準 $${5\%}$$（$${0.05}$$）より大きいです。
有意水準 $${5\%}$$ で 帰無仮説を棄却できず、データの尖度が正規分布と異なるとは言えない、です。
尖度的にデータは正規分布だろう、という感じです。

６回にわたる「４章 はじめての統計的検定」はこの記事が最後となります。
記事の最後を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の教科書です。
よかったらぜひ、お試しくださいませ。

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