R: Skew normal distribution

Skew normal distributionをあつかうテスト。

Rではsnパッケージであつかえる。Skew normal distributionの乱数を発生させて、RStanでパラメーターを推定してみる。

library(sn)
library(rstan)

set.seed(123)
x <- rsn(100, xi = 1, omega = 1, alpha = 3)
hist(x, breaks = seq(0, 3.5, 0.5))

model.text <- "
data {
  int<lower=0>  N;
  real          x[N];
}
parameters {
  real          mu;
  real<lower=0> sigma;
  real          alpha;
}
model {
  x ~ skew_normal(mu, sigma, alpha);
  mu ~ normal(0, 100);
  sigma ~ uniform(0, 100);
  alpha ~ normal(0, 100);
}
"

fit <- stan(model_code = model.text,
            data = list(N = length(x), x = x))
print(fit)

ヒストグラム

結果

Inference for Stan model: model.text.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

      mean se_mean   sd  2.5%   25%  50%  75% 97.5% n_eff Rhat
mu    0.96    0.01 0.12  0.75  0.88 0.95 1.02  1.21   452 1.00
sigma 0.96    0.00 0.11  0.75  0.89 0.96 1.04  1.19   590 1.00
alpha 3.25    0.05 1.32  1.17  2.34 3.04 3.92  6.50   714 1.00
lp__  0.29    0.04 1.24 -2.81 -0.32 0.61 1.21  1.77   823 1.01

Samples were drawn using NUTS(diag_e) at Mon Jul 14 22:31:16 2014.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

タグ：STAn RStan R