gepuro.net
gepulog

データ分析エンジニアによる備忘録的ブログ

Rstanでブログが検索結果に出ない話

pic1

で表示されないけど、

pic2

でなら表示される。

Rstan単独で検索してブログが表示される人もいるので、検索エンジンに嫌われてる可能性が高くて涙目です。

ちなみに、「stan ワイブル」で検索すれば、上位にヒットしました。

Rstanで重回帰分析

ワイブルばかりだったので、メジャーな重回帰分析をしてみた。Rstanにハマりつつある。

linear_regression.stan

data {
  int<lower=1> J; // number of data
  int<lower=1> K; // number of covariate
  real y[J];
  matrix[J,K] x;
}
parameters {
  real<lower=0> sd0;
  vector[K] beta;
  real beta0;
}
model {
  for(j in 1:J){
    increment_log_prob(normal_log(y[j], x[j] * beta + beta0, sd0));
  }
}

linear_regression.R

library(rstan)
n <- 100
x1 <- rnorm(n, mean=1, sd=0.1)
x2 <- rnorm(n, mean=1, sd=0.1)
beta1 <- 1
beta2 <- 2
beta0 <- 3
y <- beta1 * x1 + beta2 * x2 + beta0
data.stan <- list(J = n,
                  K = 2,
                  y = y,
                  x = matrix(c(x1, x2), ncol=2))

fit <- stan(file = 'linear_regression.stan', data = data.stan, 
            iter = 1000, chains = 4)

推定結果

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

           mean se_mean    sd    2.5%     25%     50%     75%   97.5% n_eff  Rhat
sd0        0.00    0.00   0.0    0.00    0.00    0.00    0.00    0.00     2  5.75
beta[1]    1.00    0.00   0.0    1.00    1.00    1.00    1.00    1.00     2  6.20
beta[2]    2.00    0.00   0.0    2.00    2.00    2.00    2.00    2.00     2  4.98
beta0      3.00    0.00   0.0    3.00    3.00    3.00    3.00    3.00     2  6.53
lp__    2244.37   95.91 136.4 2073.04 2104.57 2223.56 2344.35 2440.01     2 17.83

Samples were drawn using NUTS(diag_e) at Wed Dec 10 22:45:14 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).

Rstanでワイブル回帰の推定(説明変数複数バージョン)

matrixやvectorの使い方を覚えた。

weibull_regression.stan

data {
  int<lower=1> J; // number of data
  int<lower=1> K; // number of covariate
  int status[J];
  real<lower=0> tt[J];
  matrix[J,K] x;
}
parameters {
  real<lower=1, upper=3> m;
  vector[K] beta;
  real beta0;
}
transformed parameters {
  real<lower=0> eta[J];
  for (j in 1:J)
    eta[j] <- exp(x[j]*beta + beta0);
}
model {
  for(j in 1:J){
    if(status[j] == 1){
      increment_log_prob(weibull_log(tt[j], m, eta[j]));
    }
    if(status[j] == 0){
      increment_log_prob(weibull_ccdf_log(tt[j], m, eta[j]));
    }
  }
}

weibull_regression.R

set.seed(0)
n <- 1000 # 台数
m <- 2 # 形状パラメータ
beta1 <- -1 # 係数
beta2 <- -2 # 係数
beta0 <- 10 # 切片
x1 <- rnorm(n)
x2 <- rnorm(n)
data.df <- data.frame(tt = rweibull(n, shape = m, scale=exp(beta1 * x1 + beta2 * x2 + beta0)), 
                      x1 = x1,
                      x2 = x2,
                      status=rep(1,n))
data.df[data.df$tt > 500,]$status <- 0
data.df[data.df$tt > 500,]$tt <- 500
1-sum(data.df$status)/nrow(data.df) # 打切り率 0.921

data.stan <- list(J = nrow(data.df),
                  K = 2,
                  tt = data.df$tt,
                  x = matrix(c(data.df$x1, data.df$x2), ncol=2),
                  status = data.df$status)

fit <- stan(file = 'weibull_regression.stan', data = data.stan, 
            iter = 1000, chains = 4)

結果

    mean    se_mean sd  2.50%   25% 50% 75% 97.50%  n_eff   Rhat
m   1.92    0.01    0.18    1.58    1.79    1.92    2.04    2.26    315 1.01
beta[1] -0.92   0   0.09    -1.12   -0.98   -0.92   -0.86   -0.76   343 1.01
beta[2] -1.83   0.01    0.16    -2.17   -1.93   -1.82   -1.73   -1.55   342 1
beta0   9.69    0.02    0.32    9.12    9.45    9.67    9.89    10.4    306 1
eta[1]  8827.87 144.01  2637.51 5109.09 6978.65 8373.4  10139.67    15307.5 335 1.01
eta[2]  747.41  3.98    84.09   607.53  688.78  738.27  796.95  940.81  447 1.01
eta[3]  6494.17 95.31   1764.01 3938.58 5243.86 6200.7  7374.96 10705.78    343 1
eta[4]  70601.7 1978.38 35808.67    27994.1 46483.16    62399.01    84052.98    164023.29   328 1.01
eta[5]  188588.63   6537.63 114920.76   64135.67    113845.27   159425.4    226058.4    499006.32   309 1.01
・・・

Rstanでワイブル回帰の推定

weibull_regression.stan

data {
  int<lower=1> J; // number of data
  int status[J];
  real<lower=0> tt[J];
  real x[J];
}
parameters {
  real<lower=1, upper=3> m;
  real beta;
  real beta0;
}
transformed parameters {
  real<lower=0> eta[J];
  for (j in 1:J)
    eta[j] <- exp(beta * x[j] + beta0);
}
model {
  for(j in 1:J){
    if(status[j] == 1){
      increment_log_prob(weibull_log(tt[j], m, eta[j]));
    }
    if(status[j] == 0){
      increment_log_prob(weibull_ccdf_log(tt[j], m, eta[j]));
    }
  }
}

mを1から3ぐらいにしておかないと、推定値が遥か彼方に飛ばされることがある。ワイブル分布は厄介やね。

weibull_regression.R

library(rstan)

set.seed(0)
n <- 1000 # 台数
m <- 2 # 形状パラメータ
beta <- -1 # 係数
beta0 <- 10 # 切片
x <- rnorm(n)
data.df <- data.frame(tt = rweibull(n, shape = m, scale=exp(beta * x + beta0)), 
                      x = x, 
                      status=rep(1,n))
data.df[data.df$tt > 3000,]$status <- 0
data.df[data.df$tt > 3000,]$tt <- 3000
1-sum(data.df$status)/nrow(data.df) # 打切り率 0.929

data.stan <- list(J = nrow(data.df),
                  tt = data.df$tt,
                  x = x,
                  status = data.df$status)

fit <- stan(file = 'weibull_regression.stan', data = data.stan, 
            iter = 1000, chains = 4)

結果

    mean    se_mean sd  2.50%   25% 50% 75% 97.50%  n_eff   Rhat
m   1.94    0.01    0.22    1.53    1.8 1.93    2.08    2.39    311 1.01
beta    -1.03   0.01    0.13    -1.3    -1.11   -1.02   -0.94   -0.8    303 1.01
beta0   10.13   0.02    0.27    9.64    9.94    10.1    10.3    10.7    292 1.01
eta[1]  6901.33 47.51   887.55  5532.96 6292.8  6791.62 7381.7  8862.99 349 1.01
eta[2]  36854.06    763.24  12739.71    19890.77    28340.35    33948.03    42868.47    67365.57    279 1.01
eta[3]  6438.14 41.13   784.76  5214.99 5898.46 6343.24 6861.3  8225.53 364 1.01
eta[4]  6833.67 46.56   872.2   5486.28 6236.14 6725.15 7303.96 8769.99 351 1.01
eta[5]  16780.77    231.8   3926.93 11085.77    14076.74    16048.53    18780.13    26145.32    287 1.01
・・・・

Rstanでワイブル分布のパラメータ推定(打切りのあるデータで)

weibull.stan

data {
  int<lower=1> J; // number of data
  int status[J];
  real<lower=0> tt[J];
}
parameters {
  real<lower=1> m; 
  real<lower=1> eta;
}
model {
  for(j in 1:J){
    if(status[j] == 1){
      increment_log_prob(weibull_log(tt[j], m, eta));
    }
    if(status[j] == 0){
      increment_log_prob(weibull_ccdf_log(tt[j], m, eta));
    }
  }
}

weibull.R

サンプルデータの生成

> set.seed(0)
> data.df <- data.frame(tt = rweibull(1000, shape = 2, scale = 100),
+                       status = rep(1,1000))
> data.df[data.df$tt > 30,]$status <- 0
> data.df[data.df$tt > 30,]$tt <- 30
> 1-sum(data.df$status)/nrow(data.df) # 打切り率
[1] 0.904

Rstan用データと推定

data.stan <- list(J = nrow(data.df),
                  tt = data.df$tt,
                  status = data.df$status)

fit <- stan(file = 'weibull.stan', data = data.stan, 
            iter = 1000, chains = 4)

推定結果

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

        mean se_mean    sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
m       1.96    0.01  0.19    1.61    1.82    1.94    2.09    2.37   349 1.01
eta    99.49    0.70 13.03   77.71   89.63   98.45  108.25  128.27   351 1.01
lp__ -622.54    0.04  0.96 -625.07 -622.94 -622.28 -621.84 -621.53   566 1.00

Samples were drawn using NUTS(diag_e) at Thu Dec 04 03:30:48 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).

Rstanでワイブル分布のパラメータ推定

Rstanでワイブル分布のパラメータ推定をしてみる。打切りを考慮するにはどうしたら良いのだろうか? Stan Modeling Language User's Guide and Reference Manual, v2.5.0の83ページ辺りを読めば出来そう。

weibull.stan

data {
  int<lower=0> J; // number of data 
  real<lower=0> tt[J]; // time
}
parameters {
  real<lower=0> m; 
  real<lower=0> eta;
}
model {
  for(j in 1:J)
    tt[j] ~ weibull(m, eta);
}

weibull.R

library(rstan)
data.ls <- list(J = 1000,
                tt = rweibull(1000, shape = 2, scale = 100))

fit <- stan(file = 'weibull.stan', data = data.ls, 
            iter = 1000, chains = 4)

結果

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

         mean se_mean   sd     2.5%      25%      50%      75%    97.5% n_eff Rhat
m        2.02    0.00 0.05     1.93     1.99     2.02     2.05     2.12   946    1
eta    101.42    0.05 1.67    98.32   100.24   101.39   102.60   104.82  1007    1
lp__ -5197.22    0.04 0.97 -5199.81 -5197.58 -5196.91 -5196.55 -5196.28   735    1

Samples were drawn using NUTS(diag_e) at Wed Dec 03 23:23:05 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).