サンプル分位数の漸近正規性
サンプル分位数の漸近正規性で示した通り次が成り立ちます.
$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}
\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}_{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}_{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)$です.
よく一致していることがわかります.