rstanで自分で定義した確率分布からサンプリングする：Johnson's SU 分布

R Rpackage Rを使いこなす rstan 統計

本当はこの通りにしたかったが、自作関数のサンプリングが遅すぎたので先に正規分布normalからのサンプリングがvectorかひとつひとつかtargetかで変わるのかを検証していた。
結論から言うと組み込みのvector型サンプリングは速いが、自作関数はひとつひとつサンプリングすることになるので、圧倒的遅さ。
mikuhatsune.hatenadiary.com

これをやっていたが、実際のデータはやはり単純にガンマ分布ではあてはまりが悪かった。
mikuhatsune.hatenadiary.com
というわけで原点から最頻値が遠くて、かつ、左右非対称な分布をゴリ押ししようと思っていたら、Gumbel 分布（gumbel）でできるようだが、もっと歪ませたいと思っていたら、Johnson's SU 分布とうものがあるらしい。これは、適当な変換を挟んで正規分布に従うようにして、確率密度関数は
\frac{\delta}{\sqrt{2 \pi \lambda^2}}\frac{1}{\sqrt{1+(\frac{x-\xi}{\lambda})^2}}e^{-\frac{1}{2}*(\gamma+\delta\sinh^{-1} (\frac{x-\xi}{\lambda}))^2}
となる。これをrstanで定義すると
delta/(2*pi()*lambda^2)^0.5 * 1/(1+((x-xi)/lambda)^2)^0.5 * exp(-0.5*(gamma+delta*asinh((x-xi)/lambda))^2)
となる。

Johnson 関数群自体はSuppDistsパッケージで使える。適当に正の範囲（といっても定義域は負もある）で最頻値が50程度で左右非対称な分布を作る。パラメータを与えると、sJohnson関数で平均や分散など得られる。

sJohnson(parms)

$title
[1] "Johnson Distribution"

$gamma
[1] -5.5

$delta
[1] 2

$xi
[1] 7.5

$lambda
[1] 5

$type
SU 
 3 

$Mean
[1] 51.63246

$Median
[1] 46.44676

$Mode
[1] 37.23034

$Variance
[1] 561.298

$SD
[1] 23.69173

$ThirdCentralMoment
[1] -23119.93

$FourthCentralMoment
[1] 2784825

$PearsonsSkewness...mean.minus.mode.div.SD
[1] 0.6078967

$Skewness...sqrtB1
[1] -1.738587

$Kurtosis...B2.minus.3
[1] 5.83916


推定した結果はそこそこよい。しかし時間がクッソかかった。
f:id:MikuHatsune:20200314214536p:plain

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

library(SuppDists)
parms <- list(gamma=-5.5, delta=2, xi=7.5, lambda=5, type="SU")
n <- 3000
hoge <- rJohnson(n, parms)

code <- "
  functions{ //手作り正規分布を定義
    real johnsonSU_lpdf(real x, real gamma, real delta, real xi, real lambda){
      real tmp;
      tmp = delta/(2*pi()*lambda^2)^0.5 * 1/(1+((x-xi)/lambda)^2)^0.5 * exp(-0.5*(gamma+delta*asinh((x-xi)/lambda))^2);
      return log(tmp);
    }
  }
  data {
    int<lower=0> N;
    real<lower=0> y[N];
  }
  parameters {
  real gamma;
  real<lower=0> delta;
  real xi;
  real<lower=0> lambda;
  }
  model {
    for(i in 1:N){
    y[i] ~ johnsonSU(gamma, delta, xi, lambda);
    }
  }
"

m0 <- stan_model(model_code=code)
standata <- list(N=length(hoge), y=hoge)
fit <- sampling(m0, standata, iter=1000, warmup=500, chain=4)
ex <- extract(fit, pars=head(fit@model_pars, -1))
pa <- c(lapply(ex, median), type="SU")


x <- seq(0, 200, length=300)
d0 <- dJohnson(x, parms)
d1 <- dJohnson(x, pa)
hist(hoge, nclass=100, freq=F, main="Johnson SU", col="yellow")
lines(x, d0, col=2, lwd=4)
lines(x, d1, col=4, lwd=4)