統計検定準一級の標本抽出のところに層化抽出法があります。一見優しい様で、なかなか奥の深い内容です。まとめてみました。これは「データの分析方法を学ぶ会」の勉強会で悩んだ問題の１つです。定期的に勉強会を開いています。https://study-data-analysis.connpass.com/

コードはGitthubに掲載しました。
https://github.com/innovation1005/qiita_innovation1005

また、日曜の朝9:30と火曜の午後7:00から勉強会を開いています。
https://study-data-analysis.connpass.com/
参考文献は統計学実践ワークブックです。今回の投稿は２１章に関するものです。

層化抽出法

層化抽出法（Stratified Sampling）は、統計学において標本を抽出する方法の一つです。

母集団がその特性値に関していくつかの異なる特性値に分割できるとき、母集団をその特性値に応じて分割し、各層からある大きさの標本を個別に抽出した方が、単純無作為抽出によるよりも標本平均の推定の精度がよくなることが期待できます。特に各層の間で特性値の違いが顕著である場合に特に有用です。

層化抽出法の一般的なプロセスは以下のようになります：

層の定義: 母集団を相互に排他的な層または部分母集団に分割します。これらの層は、分析の目的に応じて選ばれ、関連する特性に基づいていることが一般的です（例：年齢、性別、収入など）。

各層からの抽出: 各層から独立に標本を抽出します。この抽出は、単純無作為抽出、非復元抽出、系統抽出など、さまざまな方法で行うことができます。

結合: 各層から抽出された標本を結合し、全体の標本を構築します。

層化抽出法は、実施が複雑であり、各層から適切に標本を抽出する必要があるため、設計と実施が困難であるとされることがあります。しかし、特定の状況下で単純無作為抽出よりも優れた推定を提供することができる強力な手法です。

無限母集団と有限母集団

無限母集団 (Infinite Population)：理論的に無限個の観測可能な要素からなる集団です。
有限母集団 (Finite Population)：有限個の観測可能な要素からなる集団です。

復元抽出と非復元抽出

復元抽出 (With Replacement)：標本を抽出した後、その標本を再び母集団に戻して次の抽出が行われます。
非復元抽出 (Without Replacement)：一度抽出された標本は母集団に戻されないため、次の抽出には影響しません。また、同じ要素が複数回選ばれることはありません。

主な2つの種類の層化抽出法

①：等配分法
各層から同じ数の観測値を抽出する方法。この方法は、全ての層が同じ重要度を持つと考える場合に使用されます。

②：比例配分法
各層の観測値の数が、母集団内での層の相対的なサイズに比例するように抽出される方法。これは、大きな層からは多くの観測値を抽出し、小さな層からは少ない観測値を抽出することになります。

等配分法がよい場合

層のサイズが異なる場合には、推定量の分散が単純無作為抽出の標本平均の分散よりも大きくなることがあります。また、各層の母平均の推定に精度を保ちたい場合に有効となります。

層内分散がほぼ同じ: 各層の内部分散がほぼ同じで、層のサイズが異なる場合、等配分法が効果的です。この場合、各層から同じ数のサンプルを抽出することで、全体的な推定精度が高まります。

特定の層の詳細な分析: 特定の層に焦点を当てた分析を行いたい場合、等配分法を使用して各層から同じサイズのサンプルを取得すると、その層の分析が容易になります。？

比例配分法がよい場合

層間分散が大きい: 各層の内部分散が大きく異なる場合、比例配分法が効果的です。この場合、大きな分散を持つ層からより多くのサンプルを抽出することで、推定の効率が向上します。

全体的な母平均の推定: 全体的な母集団の平均や合計の推定に関心がある場合、各層のサイズに比例してサンプリングすると、全体の推定がより正確になる可能性があります。

層のサイズが大きく異なる: 各層のサイズが非常に異なる場合、比例配分法は各層から適切な数のサンプルを取得し、全体的な推定のバイアスを減らすのに役立ちます。

最終的に、どちらの方法が適切かは、調査の目的、データの特性、資源の制約などを慎重に検討する必要があります。

有限母集団の非復元抽出

非復元抽出と有限修正は、実際のデータ収集や分析の際に考慮すべき重要な要因であり、母集団の大きさと標本の大きさがどのように相互作用するかを理解するのに役立ちます。特に、母集団が限られている場合や非復元抽出が行われる場合には、これらの概念が統計的推測の精度に重要な影響を及ぼす可能性があります。

Python3ではじめるシステムトレード：非復元単純無作為抽出参照

母集団の特性

各層$l$の母平均を$\mu_l$、母分散を$\sigma_l^2$とします。

母集団の全体の母平均を$\mu$、母分散を$\sigma^2$とします。

$$\mu=\sum_{l=1}^L\frac{N_l}{N}\mu_l ----------(1)$$

層内分散の和$\sigma_{l,w}^2$と各層の母分散$\sigma_l^2$との関係は

$$\sigma_{l,w}^2=\sum_{l=1}^L \frac{N_l}{N}\sigma_l^2 --------- (2)$$

ここで$N_l$は各層の大きさです。一方で、層間分散は
$$\sigma_{b}^2=\sum_{l=1}^L \frac{N_l}{N}(\mu_l-\mu)^2 --------- (3)$$

母分散、層内分散の和、層間分散の関係は
$$\sigma^2=\sigma_{l.w}^2+\sigma_{b}^2 ------------ (4)$$
となります。

有限母集団からの非復元抽出と標本平均の分散

大きさ$N$の母集団の単純無作為抽出で大きさ$n$の非復元抽出の標本平均の分散は

$$Var(\bar{Y})=\frac{N−n}{N−1}\frac{\sigma^2}{n}$$

となります。$\frac{N−n}{N−1}$を有限修正といいます。

標本分散(等配分法)

等配分法（Equal Allocation）における層化抽出は、各層から同じサイズの標本を抽出する方法です。この方法の特徴として、層間で分散が大きく異なる場合や層の大きさが大きく異なる場合に、単純無作為抽出に対して効率が悪くなることがあるという問題があります。

以下、証明を行います。

大きさ$N$の母集団をL層に分割し、各層から等しく$n_l=n/L$個ずつ抽出するとします。

等配分法における層内の標本平均の分散は、層が異なる場合の観測は独立とすると、以下のようになります。

Var(\bar{Y}_{l,w})=\sum_{l=1}^L\left(\frac{1}{L}\right)^2\frac{N_l-n/L}{N_l-1}\frac{\sigma_l^2}{n/L}

一方、単純無作為抽出の場合の標本平均の分散は

$$Var(\bar{Y})=\frac{N−n}{N−1}\frac{\sigma^2}{n}$$

となります。

ここで、等配分法と単純無作為抽出の標本分散の差を考えると、

$$dV=\sum_{l=1}^L \left(\frac{1}{L}\right)^2\frac{N_l-n/L}{N_l-1}\frac{\sigma_l^2}{n/L}-\frac{N−n}{N−1}\frac{\sigma^2}{n} -------(5)$$

となります。有限修正を無視すると

$$dV=\sum_{l=1}^L\left(\frac{1}{L}\right)^2\frac{\sigma_l^2}{n/L}-\frac{\sigma^2}{n}$$

となり

$$dV=\sum_{l=1}^L\frac{1}{L}\frac{\sigma_l^2}{n}-\sum_{l=1}^L \frac{N_l}{N}\frac{\sigma_l^2}{n}-\frac{\sigma_b^2}{n}$$
層の大きさが同じであれば
$$dV=\sum_{l=1}^L\frac{1}{L}\frac{\sigma_l^2}{n}-\sum_{l=1}^L \frac{1}{L}\frac{\sigma_l^2}{n}-\frac{\sigma_b^2}{n}$$
$$=-\frac{\sigma_b^2}{n}$$

よって、層の大きさが同じとき等配分法と単純無作為抽出との間で標本分散の差は主に層間分散の大きさに依存します。しかし、層の大きさが同じでなければ、等配分法による標本分散は単純無作為抽出によるものを上回る可能性があります。

各層の平均値に違いがある、層間分散が大きくなる場合、等配分法の方が非効率的である可能性もあります。

標本平均の精度(等配分法)

等配分法において各層の母平均の推定においても精度を保てることを証明します。ただし、「各層の母平均の推定においても精度を保てる」とはどのような意味でしょうか？

各層の母平均の推定においても精度を保てることとは、各層からの標本が全体の特性を適切に反映し、各層の特性を正確に推定できる状況を指します。層の分割において母集団内の異質性がうまくとらえられているとします。その際には

等配分法における各層内の標本平均の分散は以下のようになります。

$$Var(\bar{Y}_{l})=\frac{N_l-N_l/Nn}{N_l-1}\frac{\sigma_l^2}{N_l/Nn}=\frac{N_l-n/L}{N_l-1}\frac{\sigma_l^2}{n/L}------- (6)$$

この分散を最小にするためには$N_l/N=1/L$が必要です。これは母集団が特性値によると$L$個に分割することが適切であり、各層の大きさも等しいということをいっています。また、$n/L=N_l$とすればこの分散はゼロになりますがそれは$n=N$である必要性を示しています。つまり母集団の大きさと標本の大きさが等しいときに限り、標本平均と母平均は一致して、標本分散はゼロになります。各層の母分散は各層のデータ点の偏差の平方和を評価しています。したがってこの値が大きければ、標本平均の分散も必然的に大きくなります。しかし、有限母集団の非復元抽出の場合は有限修正された標本平均の分散については慎重な理解が必要です。まず、各層の母分散は期待値です。この期待値は有限母集団のものです。有限修正された式は$n/L=1$のとき有限修正項は１になり、

$$Var(\bar{Y}_{l})=\sigma_l^2------- (7)$$

になります。そして、標本の大きさを増やしていくと$Var(\bar{Y}_{l})$は小さくなります。これは標本の大きさが大きくなるに従い標本平均は各層の母平均に近づいていくので小さくなります。そして、$n=N$となった時に標本平均と母平均が一致して、求める分散はゼロになります。

したがって、標本平均の精度に与える要因は各層の大きさと$n$であり、各層の大きさは近ければ精度が上がり、また$n$が大きくても精度が上がり、各層の母分散が大きければ、差も大きくなる可能性があります。

標本分散(比例配分法)

比例配分法における層内の標本平均の分散は、層が異なる場合の観測は独立とすると、以下のようになります。

Var(\bar{Y}_{l,w})=\sum_{l=1}^L\frac{N_l}{N}\frac{N_l-nN_l/N}{N_l-1}\frac{\sigma_l^2}{nN_l/N}

一方、単純無作為抽出の場合の標本平均の分散は

$$Var(\bar{Y})=\frac{N−n}{N−1}\frac{\sigma^2}{n}$$

となります。

ここで、比例配分法と単純無作為抽出の標本分散の差を考えると、

$$dV=\sum_{l=1}^L \frac{N_l}{N}\frac{N_l-nN_l/N}{N_l-1}\frac{\sigma_l^2}{nN_l/N}-\frac{N−n}{N−1}\frac{\sigma^2}{n} -------(8)$$

有限修正を無視すると

$$dV=\sum_{l=1}^L\left(\frac{N_l}{N}\right)^2\frac{\sigma_l^2}{nN_l/N}-\frac{\sigma^2}{n}$$

となり

$$dV=\sum_{l=1}^L\left(\frac{N_l}{N}\right)\frac{\sigma_l^2}{n}-\sum_{l=1}^L \frac{N_l}{N}\frac{\sigma_l^2}{n}-\frac{\sigma_b^2}{n}$$
$$=-\frac{\sigma_b^2}{n}$$

よって、結果は配分方法が適切であれば、等配分法と同じになります。大きさ$n$の標本の比例配分法による標本分散と単純無作為抽出による標本分散の関係は層内・層間分散の大きさに左右され、それは層の大きさが異なれば異なるほど、また、各層の分散の大きさが異なれば異なるほど大きくなり、そのような場合には比例配分法の方が母分散を小さく見積もり、効率が悪くなります。

標本平均の精度(比例配分法)

比例配分法において各層の母平均の推定においても精度を保てることを証明します。

比例配分法における各層内の標本平均の分散は以下のようになります。

$$Var(\bar{Y}_{l})=\frac{N_l-nN_l/N}{N_l-1}\frac{\sigma_l^2}{nN_l/N}------- (9)$$

この分散の値が小さければ精度が高いといえます。$N_l=nN_l/N$のときゼロになります。これは母集団と標本が一致していることを表しています。その際には各層の標本平均は母平均と一致しているために標本平均の分散はぜるになるのです。$nN_l/N-1$のとき、各層の標本平均の分散は母平均と同じになります。標本の大きさ$n$が増えるにしたがい、標準誤差は小さくなります。

import numpy as np

def equalStr(N_l,sigma_l_sq,mu_l,n):
    L = len(N_l)
    N = sum(N_l) # 母集団の大きさ
    mu = sum(N_l * mu_l) / N # 母集団の全体の母平均
    print('層の数',L,'母集団の大きさ',N,'標本の大きさ',n)
    print('各層の大きさ',N_l)
    print('各層の分散',sigma_l_sq)
    print("各層の母平均",mu_l,"\n母平均",mu)
    
    # 等配分法でのサンプルサイズ
    n_l_eq = n / L

    # 層内分散の和と層間分散の計算
    sigma_lw_sq = sum(N_l / N * sigma_l_sq) # 式(2)
    sigma_b_sq = sum(N_l / N * (mu_l - mu)**2) # 式(3)
    # 母分散
    sigma_sq = sigma_lw_sq + sigma_b_sq # 式(4)

    print("母分散",sigma_sq,"層内母分散",sigma_lw_sq,"層間母分散",sigma_b_sq)

    # 等配分法での標本平均の分散
    Var_Y_lwe = sum((1/L) * (N_i - n/L) / (N_i - 1) * sigma / n \
               for N_i, sigma in zip(N_l, sigma_l_sq))

    # 比例配分法での標本平均の分散
    Var_Y_lwr = sum((N_i/N) * (N_i - N_i/N*n) / (N_i - 1) * sigma/n  \
               for N_i, sigma in zip(N_l, sigma_l_sq))
    # 単純無作為抽出法での標本平均の分散
    Var_Y_srs = (N - n) / (N - 1) * sigma_sq / n

    # 分散の差
    dVe = Var_Y_lwe - Var_Y_srs
    dVr = Var_Y_lwr - Var_Y_srs
    dVl=np.array([(N_i - n/L) / (N_i - 1) * sigma / n-sigma/n for N_i, sigma\
                 in zip(N_l, sigma_l_sq)])
    print("")
    print("等配分法での標本平均の分散:", Var_Y_lwe)
    print("比例配分法での標本平均の分散:", Var_Y_lwr)
    print("単純無作為抽出法での標本平均の分散:", Var_Y_srs)
    print('標本平均の分散について')
    print("等配分法　：分散の差:", dVe,"層間分散",Var_Y_lwe-Var_Y_srs)
    print("比例配分法：分散の差:", dVr,"層間分散",Var_Y_lwe-Var_Y_srs)
    print("等配分法が有利な場合:", dVe > 0)
    print("比例配分法が有利な場合:", dVr > 0)
    print('等配分法での各層の標本平均の精度')
    print('dV_l',dVl)

# 各層のサイズと分散
N_l = np.array([100, 100, 100, 100]) # 各層の大きさ
sigma_l_sq = np.array([1, 1, 1, 1]) # 各層の分散
n = 1 # 全体のサンプルサイズ
# 各層の母平均（ここではランダムに設定）
mu_l = np.array([0, 0, 0, 0])#np.random.normal(0, 1, L)

equalStr(N_l,sigma_l_sq,mu_l,n)

# 各層のサイズと分散
N_l = np.array([100, 100, 100, 100]) # 各層の大きさ
sigma_l_sq = np.array([1, 10, 100, 1000]) # 各層の分散
n = 399# 全体のサンプルサイズ
# 各層の母平均（ここではランダムに設定）
mu_l = np.array([1, 2, 3, 0])#np.random.normal(0, 1, L)

equalStr(N_l,sigma_l_sq,mu_l,n)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

pd.set_option('display.float_format', lambda x: '%.2f' % x)
# 各層のサイズ、平均、標準偏差を設定
N = [200, 100, 100]
means = [0, 0, 0]
std_devs = [1, 1, 1]
total_N = sum(N)
# サンプルサイズを設定
n = 150

# 各層からデータを生成
strata = [np.random.normal(loc=means[i], scale=std_devs[i], size=N[i]) for i \
in range(len(N))]

# 等配分法
equal_allocation = [np.random.choice(stratum, n//len(N), replace=False) for stratum\
 in strata]
equal_allocation_mean = np.mean([np.mean(stratum) for stratum in equal_allocation])

# 比例配分法
proportional_allocation/
 = [np.random.choice(stratum, int(n * Ni/total_N), replace=False)\
 for stratum, Ni in zip(strata, N)]
proportional_allocation_mean = np.mean([np.mean(stratum) * Ni/total_N for stratum, Ni\
 in zip(proportional_allocation, N)])

params=pd.DataFrame([N,means,std_devs],index=["N","means","std"])
print(params)

print("等配分法の平均:", equal_allocation_mean)
print("比例配分法の平均:", proportional_allocation_mean)
# 等配分法の有限修正
finite_correction_equal = sum([(Ni - n/len(N)) / (Ni - 1)\
    * np.var(stratum) / (n/len(N)) for stratum, Ni \
              in zip(equal_allocation, N)]) if n < total_N else 0

# 比例配分法の有限修正
finite_correction_proportional = sum([(Ni - n*Ni/total_N) / (Ni - 1)\
                 * np.var(stratum) / (n*Ni/total_N) for stratum, Ni\
                       in zip(proportional_allocation, N)]) if n < total_N else 0

print("等配分法の有限修正:", finite_correction_equal)
print("比例配分法の有限修正:", finite_correction_proportional)
plt.hist(equal_allocation)
plt.show()