[統計] 信頼区間とは何か

統計 Python

はじめに

与えられたデータの性質を知るために、回帰分析を行って係数（パラメータ）を推定するということが行われる。

このとき、得られた推定結果をp値や信頼区間、AICといったさまざまな指標で評価する。

どの評価指標を使うべきかについては、これまでに数多くの議論がされてきた。ここでその詳細は追わないが、その議論の中で、p値よりも解釈の分かりやすい信頼区間を使うべきだとの指摘もあった。

では信頼区間とは何だったか。どうやって求めるのか。これが今回のテーマである。

点推定と区間推定

母集団のパラメータをどう推定するか。まず思いつくのは、最も可能性の高そうなある1つの値を示す方法である。これは点推定と呼ばれる。

点推定 point estimation：母集団のパラメータをある1つの値で指定する方法。

しかし、推定には常に誤差が伴うものである。点推定の結果が正しいとは限らないので、推定誤差の大きさも評価したい。

そこで、推定誤差の大きさも同時に評価できるような、区間推定という方法が考えられる。

区間推定 interval estimation：真のパラメータの値が入る確率がある値 1-\alpha 以上と保証される区間 [L, U] を示す方法。ここで 1-\alpha は信頼係数と呼ばれる。

区間推定に関する上記の説明は、直観的には理解しやすい。しかし、立ち止まってよくよく考えてみると、分からないところがある。

区間推定では「真のパラメータがある区間に入る確率」なるものを考えている。真のパラメータの値は、我々が調査する前からただ1つに決まっているはずなのに、である。

例えば、目の前に密閉されて中の見えない箱が1つあったとしよう。ある調査の結果、「その箱の中に猫が入っている確率は95%以上であると保証されます」と報告された。これは一体どう理解するべきだろうか。箱の中には最初から猫がいるかいないか決まっているはずであり、箱を開ける瞬間に猫がいるかどうか確率的に決まるとは考えにくい（少なくとも我々が通常認識している時空間的スケールでは）。では、同じ箱が100個あったら、95個には猫が入っているということだろうか？しかし、目の前には箱がたった1つあるだけである。

区間推定の定義

区間推定の正確な定義を見よう。『統計学入門』のp.225から引用する。

　同一の母集団から抽出した標本でも、標本ごとに信頼区間の推定値は変化する。\theta は未知ではあるが決まった定数である。したがって、一つの標本から信頼区間を具体的な数値として推定してやれば、これは信頼区間に含まれるか含まれないかのいずれかしかない。すなわち、具体的に数値として計算した現実の信頼区間に対して、”1-\alpha の確率で \theta を含む”ということはない。信頼区間の意味は、繰り返し多くの異なった標本について信頼区間をここで述べた方法によって何回も計算した場合、\theta を区間内に含むものの割合が 1-\alpha となるというころである。

母集団から取り出された n 個の標本に基づいて、母集団のパラメータ \theta を区間推定するということを、m 回繰り返すという状況を考えてみよう。

1回目の区間推定の結果を [L_1, U_1]、i 回目の区間推定の結果を [L_i,U_i] と書くことにしよう。

すると m 回の区間推定のうち、区間内に真のパラメータ \theta が含まれる回数（\theta \in [L_i, U_i] である回数）は、m が十分大きいとき、おおよそ (1 - \alpha)m 回になるということである。

区間推定の妥当性

区間推定の正確な定義を参照した上で、区間推定の妥当性について考えてみる。

ここでは信頼係数 1-\alpha を0.95としておこう。

区間推定では「95%の確率で真のパラメータはだいたいこの区間にあるよ」ということが示される。

では「95%の確率で」というのがどういう意味かというと、区間推定の定義に立ち戻れば、「同じ母集団から今回と同じ数だけデータを取ってきて同じように区間推定した場合、95%の区間推定は真の値を含んでいるよ」ということである。

したがって、今回の区間推定が的中しているのか失敗しているのかは分からないけれども、区間推定が失敗する確率（5%）が無視できるとすれば、真のパラメータは今回の区間推定で得られた区間に入っているはずである。

以上が定義に基づく区間推定の解釈だが、正直言ってなかなかややこしい。

以下では具体的に区間推定を行う方法を説明する。

正規分布に従う母集団の平均を区間推定する（t 統計量を用いる方法）

正規分布 N(\mu, \sigma^2) に従う母集団の平均を、母集団から取り出した n 個の標本 (X_1, X_2, \dots, X_n) から推定する場合を考える。

標本平均 \bar{X} = \frac{1}{n} \sum_{i=1}^{n} X_i は t 分布に従うことが知られている。この結果を応用すると、母平均 \mu の信頼係数 1-\alpha の信頼区間は

\displaystyle [\bar{X} - t_{\alpha/2} (n - 1) s / \sqrt{n}, \bar{X} + t_{\alpha/2} (n - 1) s / \sqrt{n}]

となる。ここで s は標本分散の平方根であり、

\displaystyle s^2 = \frac{1}{n-1} \sum_{i=1}^{n}(X_i-\bar{X})^2

から計算される。また、t_{\alpha} (n - 1) は自由度 n-1 の t 分布の上側確率 100\alpha %のパーセント点である。

では実際に計算してみよう。以下のPythonコードでは、平均 \mu=2.456 かつ標準偏差 \sigma=0.327 の正規分布に従う集団から、N=40 個の標本を取り出し、信頼係数 1-\alpha=0.95 で母平均の区間推定を行うという操作を M=100 回繰り返している。

import pandas as pd
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt


# Parameters of a population
mu = 2.456
sigma = 0.327

# The number of samples drawn from the population at once
N = 40

# The number of resample and estimation
M = 100

# The confidence coefficient is (1 - alpha)
alpha = 0.05


def make_data(n):
    return pd.DataFrame({'x': np.random.normal(loc=mu, scale=sigma, size=n)})


def calc_confint(df):
    n = len(df)
    x_bar = df['x'].mean()
    t_alpha_half = stats.t.ppf(q=1 - alpha / 2, df=n - 1)
    s2 = np.sum((df['x'].values - x_bar * np.ones(n)) ** 2) / (n - 1)
    s = np.sqrt(s2)
    lower = x_bar - t_alpha_half * s / np.sqrt(n)
    upper = x_bar + t_alpha_half * s / np.sqrt(n)
    return (lower, upper)


def estimation():
    df = []
    for i in range(M):
        samples = make_data(n=N)
        confint = calc_confint(samples)
        df.append(confint)
    df = pd.DataFrame(df, columns=['lower', 'upper'])
    df.to_csv('estimation.csv', sep=',')
    return


def n_hit():
    df = pd.read_csv('estimation.csv')
    n_hit = 0
    n_miss = 0
    for lower, upper in zip(df['lower'], df['upper']):
        if (lower < mu) and (mu < upper):
            n_hit += 1
        else:
            n_miss += 1
    return n_hit


def plot():
    df = pd.read_csv('estimation.csv')
    fig, ax = plt.subplots()
    for i, (lower, upper) in enumerate(zip(df['lower'], df['upper'])):
        ax.vlines(i + 1, ymin=lower, ymax=upper, color='C0')
    ax.hlines(mu, xmin=0, xmax=M + 1, color='red')
    ax.text(x=0.95, y=0.95, s='{:d}/{:d}'.format(n_hit(), len(df)), transform=ax.transAxes, ha='right')
    fig.savefig('estimation.jpg', dpi=300)
    return


if __name__ == '__main__':
    estimation()
    plot()

結果の一例を図にして示す。

f:id:cyanatlas:20201014000003j:plain

垂直な青線で各回の区間推定の結果を表示している。水平な赤線は真の母平均を示す。青線が赤線と交点を持つ場合、その回の区間推定は真のパラメータを含んでいた（的中した）ということになる。逆に、青線が赤線と交点を持たない場合は、その回の区間推定は真のパラメータを含んでいなかった（失敗した）ということになる。

この例では、100回の区間推定のうち、95回では区間推定で得られた区間に真の値が含まれていた。この結果は、信頼係数 1 - \alphaと対応している（乱数によってばらつきがあるので、必ず100回中95回というわけではない）。

ではもっと回数を増やしてみよう。M=1000 回の例を示す。

f:id:cyanatlas:20201014002510j:plain

この場合は、1000回中967回で区間推定が的中していた。

もっと増やしてみよう。最後は10000回。

f:id:cyanatlas:20201014005904j:plain

10000回中9537回で的中していた。おおよそ95%である。

正規分布に従う母集団の平均を区間推定する（ブートストラップ法）

母集団が正規分布に従う場合、標本平均は t 分布に従うことがあらかじめ分かっていた。上の例では、それを利用して、比較的簡単な計算で区間推定を行うことができた。

しかし、標本から得られた統計量の従う分布が分かっていない場合も多い。また、標本数が多い場合は漸近的に正規分布に従うが、標本数が少ないときは正規分布に近似すると精度が悪いという場合もある。

そこで、統計量の従う分布が数学的に分かっていなくても、コンピュータで乱数を生成することで注目する統計量の分布を生成してしまおう、というのがブートストラップ法である。

ここでは、先程と同じように、正規分布 N(\mu, \sigma^2) に従う母集団の平均を、母集団から取り出した n個の標本 (X_1, X_2, \dots, X_n) から推定する事例を、ブートストラップ法を使って考えたい（もちろん、本来であれば、この事例には上で説明したような t 分布を使う方法の方が良い）。

以下のPythonコードでは、平均 \mu=2.456 かつ標準偏差 \sigma=0.327 の正規分布に従う集団から、N=40 個の標本を取り出し、信頼係数 1-\alpha=0.95 で母平均の区間推定を行うという操作を M=100 回繰り返している。

ただし、区間推定の方法は先程とは異なる。得られた N=40 個の標本の平均と（標本）分散を持つ正規分布から N=40 個の乱数を生成してその平均を計算する、という操作を N_{bootstrap} = 2000 回繰り返すことによって、N=40 個の標本の平均が従う確率分布を計算し（この分布は t 分布と同じになる、はず）、この分布に基づいて区間推定を行っている。

実行には時間がかかるので注意。

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


# Parameters of a population
mu = 2.456
sigma = 0.327

# The number of samples drawn from the population at once
N = 40

# The number of resample and estimation
M = 100

# The confidence coefficient is (1 - alpha)
alpha = 0.05


def make_data(n):
    return pd.DataFrame({'x': np.random.normal(loc=mu, scale=sigma, size=n)})


N_bootstrap = 2000


def calc_confint(df):
    samples = [np.mean(np.random.normal(loc=df['x'].mean(), scale=df['x'].std(), size=len(df))) for _ in range(N_bootstrap)]
    lower = np.quantile(samples, q=alpha / 2)
    upper = np.quantile(samples, q=1 - alpha / 2)
    return (lower, upper)


def estimation():
    df = []
    for i in tqdm(list(range(M))):
        samples = make_data(n=N)
        confint = calc_confint(samples)
        df.append(confint)
    df = pd.DataFrame(df, columns=['lower', 'upper'])
    df.to_csv('estimation.csv', sep=',')
    return


def n_hit():
    df = pd.read_csv('estimation.csv')
    n_hit = 0
    n_miss = 0
    for lower, upper in zip(df['lower'], df['upper']):
        if (lower < mu) and (mu < upper):
            n_hit += 1
        else:
            n_miss += 1
    return n_hit


def plot():
    df = pd.read_csv('estimation.csv')
    fig, ax = plt.subplots()
    for i, (lower, upper) in enumerate(zip(df['lower'], df['upper'])):
        ax.vlines(i + 1, ymin=lower, ymax=upper, color='C0')
    ax.hlines(mu, xmin=0, xmax=M + 1, color='red')
    ax.text(x=0.95, y=0.95, s='{:d}/{:d}'.format(n_hit(), len(df)), transform=ax.transAxes, ha='right')
    fig.savefig('estimation.jpg', dpi=300)
    return


if __name__ == '__main__':
    estimation()
    plot()

結果の一例を示す。おおよそ95%で的中している。

f:id:cyanatlas:20201014005200j:plain

M=1000 回の場合。こちらもおおよそ95%で的中している。

f:id:cyanatlas:20201014013759j:plain

※ M=10000 回で計算してみたところ、93.99%で的中していた。上の二つの例も踏まえると、95%よりも少ない方にバイアスがかかっているかもしれない。原因が分かったら追記したい。

まとめ

  • 区間推定で得られた区間の厳密な解釈には注意が必要（95%の”確率”で区間は真のパラメータを含んでいる）
  • 区間推定には、あらかじめ統計量の従う確率分布が分かっている場合（t 分布など）には、その分布を利用して区間推定が行われる。
  • 確率分布が分かっていない場合でも、ブートストラップ法によって区間推定を行うことができる。ただし、計算時間は長くなる。

参考資料・参考文献

関連記事

cyanatlas.hatenablog.com

cyanatlas.hatenablog.com