stanにBurr Type Xii Distributionを実装してみる

タイトル通り、stanにBurr Type Xii Distributionを実装してみる。
(以下簡略化のためBurr分布と略する。また話の都合の上、{actuar}パッケージの表記に従うものとする)
Burr分布[0,+∞)にて連続かつ単調減少の裾の重い分布。
三つのパラメータ$\alpha$,$\gamma$,$\theta$からなっており、その確率密度関数$f(x;\alpha,\gamma,\theta)$は

$$
f(x;\alpha,\gamma,\theta)=
\frac{\alpha\gamma(\frac{x}{\theta})^{\gamma}}
{x[1+(\frac{x}{\theta})^{\gamma}]^{\alpha+1}}
$$

となっている。
n個の観測値がこのBurr分布より得られたものとすると、その対数尤度を$l$とすれば、
$$
l=n(log(\alpha)+log(\gamma)-log(\theta))+(\gamma-1)\sum_{i=1}^n log(\frac{x_i}{\theta})-(\alpha+1)\sum_{i=1}^n log(1+(\frac{x_i}{\theta})^\gamma)
$$
となる。
このBurr分布をstanに実装してみる。

burr_code <- "
functions{
 real Burr_type12_lpdf(real y, real alpha, real gamma, real theta){
  return( log(alpha)+log(gamma)+gamma*log(y/theta)-log(y)-(alpha+1)*log(1+(y/theta)^gamma) );
 }
}

data { 
  int<lower=0> N; 
  real y[N]; 
} 
parameters { 
  real<lower=0> alpha;
  real<lower=0> gamma;
  real<lower=0> theta;
} 
model { 
  for(t in 1:N){
   y[t] ~ Burr_type12(alpha,gamma,theta); 
  }
} 
"

このコードを使って実際に{actuar}パッケージrburrから得られたデータを基に各パラメータを推定してみよう。

install.packages("actuar")
set.seed(3141592)
y<-rburr(500,2,3,1)
data<-list(y=y,N=length(y))

burr<-stan(model_code=burr_code,data=data,seed=3141592,chains=4,warmup=1000,iter=3500)

print(burr)

まずまず良い結果が得られたのではないだろうか。

参照

・wikipedia Burr Distribution
・Package 'actuar' - CRAN.R-project.org