パーセンタイル(四分位点など)から確率分布のパラメータを推定する方法の学習記録

【学習動機】

文献値のパーセンタイル(四分位点や中央値(メディアン)など)の情報から確率分布のパラメータを推定する必要があった. 使えるパーセンタイルの数とパラメータの数が一致していて連立方程式が解けるならば, その解を推定値とすればよいが, 使えるパーセンタイルの数がパラメータの数より多い場合, どのように推定するか悩んでいた. 色々とググっていたら, {rriskDistributions}パッケージ(https://cran.r-project.org/web/packages/rriskDistributions/rriskDistributions.pdf) を見つけた. 分位点情報からパラメータを推定できるようだ. このパッケージの関数の中身を見て, やっていることは理解したので, その推定の仕方を真似しながら, 推定の精度がどのようなものになるのかをシミュレーションで実験してみることにする. {rriskDistributions}の使い方は簡単なので, ここでは書かない.

【学習記録】

宇宙からパッケージ読み込み.

library(tidyverse)  # パッケージ

パラメータ\(\lambda=10\)の指数分布からサンプルサイズ\(N=50\)の乱数を生成.

N <- 50
lambda <- 10

# 指数分布
x <- rexp(n=N, rate=lambda)
summary(x)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.003677 0.033519 0.060534 0.088263 0.114248 0.343835

分位点情報からパラメータを推定するために損失関数を設定し, これを最小化してパラメータを推定する.

\((p_1, q_1), (p_2, q_2), ...\)のように\(p_i * 100\)%点が\(q_i\)であるという情報があり,
\(P_{exp}(q; \lambda)\)が\(q\)のパラメータ\(\lambda\)の指数分布における確率点を表すものとして,
損失関数は以下のように定義し,

これを最小化することでパラメータを推定する. (過剰識別のケースでGMMって流れに似てる?)

# 損失関数の設定
lossfunc <- function(lambda) {
  s <- stats::pexp(q = q, rate=lambda) - p
  sum(s^2)
}

以下で5パターンの推定を試す.

① 25%点, 50%点, 75%点の3点のパーセンタイル情報を使った推定

# 25%点, 50%点, 75%点
p <- c(0.25, 0.5 ,0.75)
q <- quantile(x, probs = p)

fit <- stats::optim(par = 1, lossfunc, 
                    method = "L-BFGS-B", lower = c(0.001, 0.001), 
                    upper = c(10000,10000))
fit$par

## [1] 11.07593

② メディアンのみの情報を使った推定

# 50%点
p <- c(0.5)
q <- quantile(x, probs = p)

fit <- stats::optim(par = 1, lossfunc, 
                    method = "L-BFGS-B", lower = c(0.001, 0.001), 
                    upper = c(10000,10000))
fit$par

## [1] 11.45052

③ 2.5%点, 50%点, 97.5%点の3点のパーセンタイル情報を使った推定

# 2.5%点, 50%点, 97.5%点
p <- c(0.025, 0.5 ,0.975)
q <- quantile(x, probs = p)

fit <- stats::optim(par = 1, lossfunc, 
                    method = "L-BFGS-B", lower = c(0.001, 0.001), 
                    upper = c(10000,10000))
fit$par

## [1] 11.05005

④ 2.5%点, 25%点, 50%点, 75%点, 97.5%点の5点のパーセンタイル情報を使った推定

# 2.5%点, 25%点, 50%点, 75%点, 97.5%点
p <- c(0.025, 0.25, 0.5, 0.75 ,0.975)
q <- quantile(x, probs = p)

fit <- stats::optim(par = 1, lossfunc, 
                    method = "L-BFGS-B", lower = c(0.001, 0.001), 
                    upper = c(10000,10000))
fit$par

## [1] 10.94071

⑤ 平均の情報のみを使ってモーメント法で推定(これはパーセンタイルの情報を使わないが比較のため)

# 平均の逆数(モーメント法)
1/mean(x)

## [1] 11.32973

さて, 上でやったようなことを10000回のシミュレーションで繰り返して, 推定の精度がどのようになっているのか確認する.

set.seed(20220309)
nSim <- 10000

est_quant = numeric(nSim)  # 四分位点
est_med = numeric(nSim)  # 中央値
est_95 = numeric(nSim)  # 95%区間
est_5p = numeric(nSim)  # 5点
est_mean = numeric(nSim) # 平均


for(i in 1:nSim){
  x <- rexp(n=N, rate=lambda)

  # (25%点, 50%点, 75%点)
  p <- c(0.25, 0.5 ,0.75)
  q <- quantile(x, probs = p)
  
  fit <- stats::optim(par = 1, lossfunc, 
                      method = "L-BFGS-B", lower = c(0.001, 0.001), 
                      upper = c(10000,10000))
  est_quant[i] <- fit$par
  
  # (50%点)
  p <- c(0.5)
  q <- quantile(x, probs = p)
  
  fit <- stats::optim(par = 1, lossfunc, 
                      method = "L-BFGS-B", lower = c(0.001, 0.001), 
                      upper = c(10000,10000))
  est_med[i] <- fit$par
  
  # (2.5%点, 50%点, 97.5%点)
  p <- c(0.025, 0.5 ,0.975)
  q <- quantile(x, probs = p)
  
  fit <- stats::optim(par = 1, lossfunc, 
                      method = "L-BFGS-B", lower = c(0.001, 0.001), 
                      upper = c(10000,10000))
  est_95[i] <- fit$par
  
  # (2.5%点, 25%点, 50%点, 75%点, 97.5%点)
  p <- c(0.025, 0.25, 0.5, 0.75 ,0.975)
  q <- quantile(x, probs = p)
  
  fit <- stats::optim(par = 1, lossfunc, 
                      method = "L-BFGS-B", lower = c(0.001, 0.001), 
                      upper = c(10000,10000))
  est_5p[i] <- fit$par
  
  # 平均の逆数(モーメント法)
  est_mean[i] <- 1/mean(x)
}

res_est <- data.frame(
  est_quant = est_quant,  # 四分位点
  est_med = est_med,  # 中央値
  est_95 = est_95,  # 95%区間
  est_5p = est_5p,  # 5点
  est_mean = est_mean # 平均
)

res_est_long <- res_est %>% 
  pivot_longer(cols = c(est_quant,est_med,est_95,est_5p,est_mean)) %>% 
  transform(name = factor(name, levels = c("est_quant","est_med", "est_95", "est_5p", "est_mean")))

上から順番に, 3つの四分位点, メディアンのみ, 95%区間の両端の2点, 3つの四分位点と95%区間の両端の5点, 平均を使った推定の結果の推定値のヒストグラムを示す.

ggplot(res_est_long) + geom_histogram(aes(x=value),bins = 50) +
  facet_wrap(facets = ~name, ncol = 1) +
  geom_vline(aes(xintercept=10), colour="red") +
  theme_bw()

上から順番に, 3つの四分位点, メディアンのみ, 95%区間の両端の2点, 3つの四分位点と95%区間の両端の5点, 平均を使った推定の結果の推定値の要約統計量を示す.

res_est_long %>% group_by(name) %>% 
  summarise(mean(value),
            var(value),
            min(value),
            max(value))

## # A tibble: 5 x 5
##   name      `mean(value)` `var(value)` `min(value)` `max(value)`
##   <fct>             <dbl>        <dbl>        <dbl>        <dbl>
## 1 est_quant          10.2         3.07         5.59         20.4
## 2 est_med            10.2         4.45         5.07         23.4
## 3 est_95             10.4         3.88         5.76         23.2
## 4 est_5p             10.2         2.92         5.89         20.4
## 5 est_mean           10.2         2.22         6.25         18.5

平均を使った推定が最も精度が良く, 5点のパーセンタイルを使ったものが次点で良かった. メディアンのみだと1点のみだから最も精度が悪かった.
2点のパーセンタイルを使っているものでは, 25%点と75%点を使った推定の方が, 95%区間の両端を使った推定より精度が良かった.
これは, 内側の順位統計量の方がパラメータの推定に使える情報を豊富に持っているということだろうか. このへんの理論や解釈はよく分かっていない. 詳しい人がいたら教えてほしい.

【学習予定】

今のところはとくになし. 分位点を使った推定にもう少し詳しくなったら何か書くかも.