サンプル分位数の漸近正規性を確かめる

$\newcommand{\lnl}{\\[8pt]}$ $\newcommand{\Lnl}{\\[18pt]}$ $\newcommand{\delt}{\mathrm{d}}$ $\newcommand{\comb}{\mathrm{C}}$ $\DeclareMathOperator*{\ssum}{\Sigma}$ $\DeclareMathOperator*{\sprod}{\Pi}$

サンプル分位数の漸近正規性

サンプル分位数の漸近正規性で示した通り次が成り立ちます.

$0 < p < 1$とし, 確率分布関数$F$が密度関数$f$を$Q_p$の近傍で持ち, $f$が$Q_p$で正であり連続ならば,
\begin{align}
\hat{Q}_{pn} \xrightarrow{d} \mathrm{N}\left(Q_p,\frac{p(1-p)}{f(Q_p)^2n}\right)
\end{align}

となる.

今回は, この定理が本当に成り立っているのかRを使ってシミュレーションしていきましょう.

シミュレーション

$F$が正規分布の場合

$F$を$\mathrm{N}(50,10^2)$に従う場合でシミュレーションしてみます.
$n= 100 , p = 0.3$とします.
この場合, $Q_p \fallingdotseq 44.756 , f(Q_p) \fallingdotseq 0.03477$なので,

\begin{align}
\hat{Q}_{0.3} \xrightarrow{d} \mathrm{N}\left(44.756,1.737\right)
\end{align}

となるはずです.シミュレーションにより$\hat{Q}_p$を独立に$10000$個生成し上記の正規分布に従うか確認してみます.

Rのコードは次のようになります.

loop <- 1:10000 # Qpを計算する回数
n <- 100 #各回で利用するサンプル数
p <- 0.3 #分位点の位置

mean <- 50 # Fの従う正規分布の平均
sd <- 10 # Fの従う正規分布の標準偏差

qp <- qnorm(mean = mean,sd = sd,p = p) # 真のQpの値
f_qp <- dnorm(mean=mean,sd=sd,x = qp) #真のf(Qp)の値

qps = c() # 各回のQpを格納
for(i in loop){
  # サンプルをn個作成
  x <- rnorm(n,mean=mean,sd=sd)
  z = sort(x)[ceiling(n*p)] #標本分位点
  qps <- rbind(qps,c(z))
}
plot(density(qps)) #得られたQpの分布から密度関数を推定しプロット
# 収束先の分布をプロット
curve(dnorm(x,mean=qp,sd=sqrt(p*(1-p)/n)/f_qp),add = T,col="red")

実行結果:

黒線がシミュレーションで生成した$\hat{Q}_p$で赤線が収束先の分布, つまり$\mathrm{N}\left(44.756,1.737\right)$です.
よく一致していることがわかります.

$F$がカイ二乗分布の場合

$F$を$\chi^2(8)$に従う場合でシミュレーションしてみます.
$n= 200 , p = 0.7$とします.
この場合, $Q_p \fallingdotseq 9.524 , f(Q_p) \fallingdotseq 0.0769$なので,

\begin{align}
\hat{Q}_{0.7} \xrightarrow{d} \mathrm{N}\left(9.524,0.1775\right)
\end{align}

となるはずです.シミュレーションにより$\hat{Q}_p$を独立に$10000$個生成し上記の正規分布に従うか確認してみます.

Rのコードは次のようになります.

loop <- 1:10000 # Qpを計算する回数
n <- 200 #各回で利用するサンプル数
p <- 0.7 #分位点の位置

df <- 8 # Fの従うカイ二乗分布の自由度

qp <- qchisq(df = df,p = p) # 真のQpの値
f_qp <- dchisq(df = df,x = qp) #真のf(Qp)の値


qps = c() # 各回のQpを格納
for(i in loop){
  # サンプルをn個作成
  x <- rchisq(n,df = df)
  
  z = sort(x)[ceiling(n*p)] #標本分位点
  qps <- rbind(qps,c(z))
}
plot(density(qps)) #得られたQpの分布から密度関数を推定しプロット
# 収束先の分布をプロット
curve(dnorm(x,mean=qp,sd=sqrt(p*(1-p)/n)/f_qp),add = T,col="red")

実行結果:

黒線がシミュレーションで生成した$\hat{Q}_p$で赤線が収束先の分布, つまり$\mathrm{N}\left(9.524,0.1775\right)$です.
よく一致していることがわかります.