Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

スコア検定

概要

スコア関数

分布f(x∣θ)f(x|\theta)からランダムに得られたnn個のサンプルをX=(X1,⋯ ,Xn)X = (X_1, \cdots, X_n)とする。

XXの同時確率関数（尤度）を∏i=1nf(xi∣θ)\prod_{i=1}^n f(x_i|\theta)をfn(x∣θ)f_n(x|\theta)と表す。

対数尤度l(θ,X)=log⁡fn(X∣θ)l(\theta, X) = \log f_n(X|\theta)の導関数

ln′(θ,X)=ddθlog⁡fn(X∣θ)l'_n(\theta, X) = \frac{d}{d \theta} \log f_n(X|\theta)
(1)

を スコア関数 （score function）という。

フィッシャー情報量

スコア関数の2乗の期待値

In(θ)=E[{ln′(θ,X)}2]=E[{ddθlog⁡fn(X∣θ)}2]I_n(\theta)=E\left[\left\{l'_n(\theta, \boldsymbol{X})\right\}^2\right]=E\left[\left\{\frac{d}{d \theta} \log f_n(\boldsymbol{X} \mid \theta)\right\}^2\right]
(2)

を フィッシャー情報量 という。

スコア関数の性質

E[ln′(θ,X)]=0Var[ln′(θ,X)]=−E[ln′′(θ,X)]=In(θ)=nI1(θ)E[l'_n(\theta, X)] = 0\\ Var[l'_n(\theta, X)] = -E[l''_n(\theta,X)]= I_n(\theta) = n I_1(\theta)
(3)

（I1(θ)I_1(\theta)はサンプル1つのフィッシャー情報量）

となるので、中心極限定理により、θ\thetaが真値のとき

l′(θ,X)/nI1(θ)→dN(0,1)l'(\theta, \boldsymbol{X}) / \sqrt{n I_1(\theta)} \rightarrow_d \mathcal{N}(0,1)
(4)

という漸近正規性があるので、これを使って検定ができる

TODO: 導出を追う

l1′(θ,Xi)={(d/dθ)f(x∣θ)}/f(x∣θ)l'_1(\theta, X_i) = \{(d/d\theta) f(x|\theta)\} / f(x|\theta)と書けるらしいので（？）

E[l1′(θ,Xi)]=∫{(d/dθ)f(x∣θ)f(x∣θ)}f(x∣θ)dx=∫ddθf(x∣θ)dx=ddθ∫f(x∣θ)dx=ddθ1=0\begin{aligned} E\left[l'_1\left(\theta, X_i\right)\right] & =\int\left\{\frac{(d / d \theta) f(x \mid \theta)}{f(x \mid \theta)}\right\} f(x \mid \theta) d x \\ & =\int \frac{d}{d \theta} f(x \mid \theta) d x=\frac{d}{d \theta} \int f(x \mid \theta) d x=\frac{d}{d \theta} 1=0 \end{aligned}
(5)

スコア検定

H0:θ=θ0v.s.H1:θ≠θ0H_0: \theta = \theta_0 \quad \text{v.s.} \quad H_1: \theta \neq \theta_0
(6)

に対して、

R={x∈X∣∣ln′(θ0,x)∣/nI1(θ0)≥zα/2}R=\left\{\boldsymbol{x} \in \mathcal{X} \mid | l'_n\left(\theta_0, \boldsymbol{x}\right) \mid / \sqrt{n I_1\left(\theta_0\right)} \geq z_{\alpha / 2}\right\}
(7)

を棄却域とする検定を スコア検定 （score test） という。

例

参考：Julia and real estate statistics

ガンマ分布っぽいデータがあったとする。

Source
import numpy as np

np.random.seed(0)
n = 1000
alpha = 1.1  # 知りたいパラメータ
scale = 1    # 今回は既知とする
x = np.random.gamma(shape=alpha, scale=scale, size=n)

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=[4,2])
ax.hist(x, bins=30)
ax.set(xlabel="Value", ylabel="Frequency")
fig.show()
<Figure size 400x200 with 1 Axes>

ガンマ分布のshapeパラメータα\alphaは、α=1\alpha=1のとき指数分布と一致する。

データが指数分布（α=1\alpha=1）に従っているのか、かガンマ分布（α≠1\alpha \neq 1）なのか検定する（ガンマ分布のscaleパラメータは既知とする）

最尤推定量αML\alpha_{ML}について以下の仮説に対する検定を行う

H0:αML=1v.s.H1:αML≠1H_0: \alpha_{ML} = 1 \quad \text{v.s.} \quad H_1: \alpha_{ML} \neq 1
(8)

最尤推定量αML\alpha_{ML}を求める

from scipy.stats import gamma

def log_likelihood(alpha):
    return np.sum(gamma.logpdf(x, a=alpha, scale=scale))

def neg_log_likelihood(alpha, scale, data):
    # scipyはminimize_scalar しかないので、対数尤度ではなく負の対数尤度にする
    return -np.sum(gamma.logpdf(data, a=alpha, scale=scale))

from scipy.optimize import minimize_scalar
result = minimize_scalar(neg_log_likelihood, args=(scale, x), bounds=(0.01, 10), method='bounded')
alpha_hat = result.x

# 最尤推定量
print(f"{alpha_hat=:.3f}")
alpha_hat=1.075

E[−ln′′(θ,X)]=In(θ)E[-l''_n(\theta,X)]= I_n(\theta)という関係があるので、対数尤度の2階微分でフィッシャー情報量を求める

import numdifftools as nd
grad = nd.Derivative(log_likelihood, n=1)
hess = nd.Derivative(log_likelihood, n=2)
print(f"{grad(alpha_hat)=}, {hess(alpha_hat)=}")
grad(alpha_hat)=array(0.0002217), hess(alpha_hat)=array(-1480.84211289)

なぜフィッシャー情報量で割るのか → 対数尤度のピークが鋭ければl′′(θ)l''(\theta)は大きな値をとるためらしい

Source
fig, axes = plt.subplots(figsize=[6, 6], nrows=4, sharex=True)

alphas = np.linspace(0.01, 3, 100)
axes[0].plot(alphas, [log_likelihood(a) for a in alphas])
axes[0].axvline(alpha_hat, color="darkorange")
axes[0].text(alpha_hat, log_likelihood(alpha_hat)*2, r" $\hat{\alpha}$", color="darkorange")
axes[0].set(ylabel=r"log likelihood $\ell(\alpha)$",
           title=fr"log likelihood and derivatives (true $\alpha$ = {alpha})")

scores = [grad(a).tolist() for a in alphas]
axes[1].plot(alphas, scores)
axes[1].set(ylabel=r"$\ell'(\alpha)$")

infos = [hess(a).tolist() for a in alphas]
axes[2].plot(alphas, infos)
axes[2].set(ylabel=r"$\ell''(\alpha)$")

test_stats = [(abs(score)/np.sqrt(-info)) for score, info in zip(scores, infos)]
axes[3].plot(alphas, test_stats)
axes[3].set(xlabel=r"$\alpha$", ylabel=r"$\frac{|\ell'(\alpha)|}{ \sqrt{-\ell''(\alpha)}}$")

fig.show()
/usr/local/lib/python3.10/site-packages/numpy/lib/nanfunctions.py:1563: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a,
/usr/local/lib/python3.10/site-packages/numdifftools/limits.py:150: UserWarning: All-NaN slice encountered
  warnings.warn(str(msg))
<Figure size 600x600 with 4 Axes>

∣ℓ′(α)∣−ℓ′′(α)\frac{|\ell'(\alpha)|}{ \sqrt{-\ell''(\alpha)}}の値が大きいので標準正規分布に従っているのか不安だが、概ねそれっぽい値になってる

α=1\alpha = 1とする（α0=1\alpha_0=1とする）という帰無仮説のもとでの検定統計量を計算

# 検定統計量
alpha_0 = 1

I = -hess(alpha_0)  # I = E[-l''(θ)]
score = abs(grad(alpha_0)) / np.sqrt(I)
print(f"|l'(α0)| / √I(α0) = {score:.3f}")
|l'(α0)| / √I(α0) = 2.894

これは標準正規分布上では次のような値になる

Source
from scipy.stats import norm

fig, ax = plt.subplots(figsize=[4, 2])

# 検定統計量が従う標準正規分布
x_for_plot = np.linspace(-3, 3, 100)
z = norm.pdf(x_for_plot, loc=0, scale=1)
ax.plot(x_for_plot, z)
clevel = 0.05 / 2
for a in [clevel, (1 - clevel)]:
    ppf = norm.ppf(a)
    ax.axvline(x=ppf, color="steelblue", linestyle="--")
    if ppf < 0:
        ax.text(ppf - 0.1, norm.pdf(ppf) + 0.01, "$Z_{a/2}$", color="steelblue", horizontalalignment="right")
    else:
        ax.text(ppf + 0.1, norm.pdf(ppf) + 0.01, "$Z_{1-a/2}$", color="steelblue", horizontalalignment="left")
ax.set(xlabel="x", ylabel="density", title=f"standard normal distribution $N(\mu=0, \sigma=1)$")

# 検定統計量
ax.axvline(score, color="darkorange")
ax.text(score, 0.3, f"{score=:.3f}", color="darkorange", horizontalalignment="right")

fig.show()
<Figure size 400x200 with 1 Axes>
from scipy.stats import norm

pvalue = 2 * (1 - norm.cdf(score))
pvalue
0.0038054754122531786
統計学
Wald検定、Score検定、尤度比検定の関係
統計学
Bootstrap