ベイズマックスへようこそ。
統計に関するあれやこれやを自由気ままに投稿します。
全記事一覧はこちらから。

2020年10月02日

一般化ガンマ分布をStanに実装する

 最近は生存時間解析 (survival analysis; 継続時間解析とも) を使った研究に勤しんでいるのですが,パラメトリックモデルについて勉強していたら,ガンマ分布・ワイブル分布・指数分布を包含する一般化ガンマ分布 (generalized gamma distribution) なる面白い分布が登場したので,これをStanに実装してみました。一般化ガンマ分布は非負の値を返す確率分布で,確率密度関数は以下の式で表されます。

ggamma_formula.png

この式ではrateパラメータをλ (逆数の1/λがscaleパラメータ),2つの形状パラメータをg,wと置いています。gはガンマ分布の形状パラメータ,wはワイブル分布の形状パラメータに対応しています。g=1のときワイブル分布,w=1のときガンマ分布,g=w=1のとき指数分布と一致することが特徴です。他のパラメータ化の方法もありますが,他の分布との関係が分かりやすいのでこのパラメータ化を採用しています。分布間の関係は次のように図示できます。

ggamma_family.png

一般化ガンマ分布を使えば,データがガンマ分布・ワイブル分布・指数分布のどれに従うのかが分からない場合でも分布形に当たりを付けることができたりします。もちろん,情報量規準などを用いて4つの分布のどれが適切かを判断することもできます。便利ですね。

 一般化ガンマ分布を使ったモデリングができるようにStanコードを書いてみました。

ggamma.stan:一般化ガンマ分布モデル
functions{
  //for文で1つずつ計算するときはこっち
  real ggamma_lpdf(real x, real rate, real g, real w){
    real lpdf;
    lpdf = log(w) + g*w*log(rate) + (g*w-1)*log(x) - (x*rate)^w - lgamma(g);
    return(lpdf);
  }
  
  //データをベクトルで渡すときにはこっち
  real sum_ggamma_lpdf(vector x, real rate, real g, real w){
    real sum_lpdf;
    sum_lpdf = sum(log(w) + g*w*log(rate) + (g*w-1)*log(x) - exp(w*log(x*rate)) - lgamma(g));
    return(sum_lpdf);
  }
}

data{
  int N;
  vector[N] Y;
}

parameters{
  real rate;
  real g;
  real w;
}

model{
  //事前分布(シミュレーションしやすいように狭めにしています)
  rate ~ uniform(0,10);
  g ~ uniform(0,10);
  w ~ uniform(0,10);
  
  //モデル
  target += sum_ggamma_lpdf(Y | rate, g, w);
  
  ////以下のような書き方でも同じ
  //for(i in 1:N) Y[i] ~ ggamma(rate, g, w);
}

 上記コードではggamma_lpdf()とsum_ggamma_lpdf()という2種類の関数を定義していますが,サンプルを要素ごとにreal型としてfor文で渡すときにはggamma_lpdf()を,サンプルをvectorとしてまとめて渡すときにはsum_ggamma_lpdf()を使います。後者では対数尤度をまとめて計算してその和を求めているという訳です。推定結果はどちらを使っても同じですが,データの数が多い時にはおそらくsum_ggamma_lpdf()のほうが計算が速いです(ちゃんと検証した訳ではありませんが)。では,このコードでうまくパラメータが推定できるか確認してみましょう。今回は,一般化ガンマ分布(λ = 0.5,g = 3,w = 5),ガンマ分布(λ = 0.5,g = 3),ワイブル分布(λ = 0.5,w = 5),指数分布(λ = 0.5)からの乱数を2,000個ずつ発生させて4つのデータを作り,それらに一般化ガンマ分布を当てはめてパラメータを推定してみます。Rコードと推定結果は以下の通りです。なお,一般化ガンマ分布をRで扱うにはggammaというパッケージが利用できます(gamlssというパッケージもありますが,パラメータ化の方法が異なります)。

ggamma.R:一般化ガンマ分布モデルの確認
#R 3.6.3を使用
library(tidyverse) #version 1.2.1
library(rstan) #version 2.19.2
library(ggamma) #version 1.0.1

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

## 架空のデータを作る ----

N <- 2000
rate <- 0.5
shape_g <- 3
shape_w <- 5

set.seed(8823)
Y_ggamma <- rggamma(N, a = 1/rate, b = shape_w, k = shape_g)
Y_gamma <- rgamma(N, rate = rate, shape = shape_g)
Y_weibull <- rweibull(N, scale = 1/rate, shape = shape_w)
Y_exp <- rexp(N, rate = rate)

## Stanに渡す ----
sm <- stan_model("stanmodel/ggamma.stan")

warmup <- 1000
iter <- 1000 + warmup
t0<-proc.time()
fit_ggamma <- sampling(sm,data=list(N=N,Y=Y_ggamma)
                ,seed=8823,chains=4,iter=iter,warmup=warmup
)

fit_gamma <- sampling(sm,data=list(N=N,Y=Y_gamma)
                       ,seed=8823,chains=4,iter=iter,warmup=warmup
)

fit_weibull <- sampling(sm,data=list(N=N,Y=Y_weibull)
                       ,seed=8823,chains=4,iter=iter,warmup=warmup
)

fit_exp <- sampling(sm,data=list(N=N,Y=Y_exp)
                       ,seed=8823,chains=4,iter=iter,warmup=warmup
)

## 一般化ガンマ分布 (rate = 0.5, g = 3, w = 5) からの乱数の結果 ----
s_ggamma <- summary(fit_ggamma)$summary
round(s_ggamma,3)

#          mean se_mean    sd     2.5%      25%      50%      75%    97.5%   n_eff  Rhat
# rate    0.511   0.002 0.037    0.458    0.484    0.504    0.532    0.598 488.966 1.008
# g       3.271   0.030 0.667    2.276    2.777    3.161    3.663    4.825 492.162 1.008
# w       4.992   0.023 0.511    3.982    4.628    4.996    5.364    5.974 493.785 1.009
# lp__ -332.568   0.040 1.214 -335.789 -333.102 -332.264 -331.696 -331.168 920.305 1.001

## ガンマ分布 (rate = 0.5, g = 3, w = 1) からの乱数の結果 ----
s_gamma <- summary(fit_gamma)$summary
round(s_gamma,3)

#           mean se_mean    sd      2.5%       25%       50%       75%     97.5%   n_eff  Rhat
# rate     0.873   0.031 0.595     0.372     0.569     0.731     0.979     2.186 370.082 1.017
# g        3.695   0.039 0.830     2.426     3.142     3.575     4.111     5.644 464.461 1.012
# w        0.878   0.004 0.099     0.686     0.813     0.877     0.941     1.081 496.439 1.010
# lp__ -5125.935   0.038 1.201 -5129.165 -5126.483 -5125.623 -5125.061 -5124.577 982.829 1.004

## ワイブル分布 (rate = 0.5, g = 1, w = 5) からの乱数の結果 ----
s_weibull <- summary(fit_weibull)$summary
round(s_weibull,3)

#         mean se_mean    sd      2.5%       25%       50%       75%     97.5%    n_eff  Rhat
#rate     0.502   0.001 0.017     0.475     0.491     0.501     0.512     0.540  515.125 1.007
#g        1.036   0.006 0.131     0.814     0.942     1.027     1.110     1.327  515.630 1.008
#w        4.999   0.017 0.388     4.234     4.738     4.987     5.254     5.772  529.559 1.007
#lp__ -1076.843   0.040 1.281 -1080.196 -1077.396 -1076.519 -1075.901 -1075.407 1020.983 1.003

## 指数分布 (rate = 0.5, g = 1, w = 1) からの乱数の結果 ----
s_exp <- summary(fit_exp)$summary
round(s_exp,3)

#           mean se_mean    sd      2.5%       25%       50%       75%     97.5%   n_eff  Rhat
# rate     0.574   0.004 0.103     0.425     0.506     0.560     0.624     0.814 553.282 1.005
# g        1.066   0.006 0.136     0.838     0.975     1.059     1.144     1.362 581.750 1.005
# w        0.982   0.003 0.076     0.837     0.932     0.980     1.030     1.132 598.027 1.004
# lp__ -3293.357   0.042 1.288 -3296.725 -3293.945 -3293.026 -3292.422 -3291.911 932.671 1.002

 サンプルサイズや乱数の関係でEAP (mean) がやや真値から離れているように見える箇所もありますが,いずれの場合でも95%確信区間に真値が含まれていることが確認できます。特に,ガンマ分布ではw≃1,ワイブル分布ではg≃1,指数分布ではg≃w≃1と正しく推定できていることが分かります。というわけで,ぜひみなさんも一般化ガンマ分布で遊んでみてください。※この記事で紹介したコードを研究等に使用される場合は自己責任でお願いします。

posted by mutopsy at 14:05 | 統計