RとStanを用いた指数-ガウス分布の回帰モデル
意味のある外れ値に対するモデリング【試論】


ExGaussPreface1.png ExGaaussRegPred.png ExGaussPreface2.png

はじめに

本稿では、反応時間データへのフィッティングなどに用いられる指数-ガウス分布(Ex-Gaussian distribution, Exponential Modified Gaussian distribution)と その回帰モデルについてご紹介します。

指数ガウス分布とは、指数分布と正規分布が重なったような連続型確率分布です。使い所は、”右に裾が長く”,”データの中に外れ値がふくまれている”と仮定できる場合でしょうか。

右に伸びる分布の裾を指数分布のパラメータで捉え、 メインの山の部分を正規分布で表現します。言い換えれば、「外れ値の部分」と「正規分布の部分」を分けて考えることができる分布です。

(外れ値の部分に、指数分布を仮定していますので、外れ値生成メカニズムの最頻値が0という仮定になることはお忘れなく。その辺の仮定が合うかどうかはデータ次第)

この性質をうまく使えば、新しい研究仮説や解析アイデアにつながるかもしれません。色々と、コードを書いて感触を掴みましょう。

指数-ガウス分布

exGaussImage.png
Figure 指数-ガウス分布の一例(mu =10, sigma = 2, rate=0.09(nu = 1/0.09))

指数-ガウス分布は、指数分布と正規分布の混合分布。精確には、指数分布と正規分布の積を畳み込み積分した確率分布です。

正規分布のパラメータである平均と分散、指数分布のrateパラメータ(stanの引数ではlambda)の3つのパラメータを持ちます

その性質として、指数分布と正規分布から生成された乱数値の和は、指数ガウス分布からの乱数と一致します。

少し、数値実験をしてみます。正規乱数と指数乱数の和と、指数-ガウス乱数からの乱数を比較してみます。

この時、いくつかの関数やパッケージで、指数分布のrateパラメータの逆数を引数にする場合などがありますのでご留意ください。

library(gamlss.dist) #指数-ガウス分布を扱うパッケージ
library(tidyverse) #データ成形、可視化関連パッケージ群
library(rstan)

#stanを並列化する
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

#正規分布と指数分布の乱数を生成
df_norm<-data.frame(norm = rnorm(2000, mean = 10, sd =1), exp = rexp(2000,rate = 0.5))

#指数-ガウス分布からの乱数を生成
df_exGauss<-data.frame(exGauss = rexGAUS(4000,mu=10,sigma=1,nu=1/0.5)) #rateと逆数を指定

#データフレームを結合して、正規分布と指数分布の和のデータを追加
df_bind<-cbind(df_norm,df_exGauss)%>%
    dplyr::mutate(norm_exp = norm + exp)%>%
    tidyr::gather()

#可視化
df_bind2<-transform(df_bind,key = factor(key,levels = c("norm","exp","norm_exp","exGauss")))
ggplot(df_bind2,aes(x=value,colour=key,fill=key))+#ggplot2で描画
    geom_histogram(alpha=0.5)+
    geom_vline(xintercept = 0,linetype=2)+
    theme_bw()+facet_wrap(~key)

よっこらせっ

exGaussRandom.png

上段が正規分布と指数分布。

その和の分布が左下。

おなじパラメータの値で、指数ガウス分布からの乱数を生成したのが、右下です。

下段の左右の分布がとても類似していることが確認されます。

指数-ガウス分布のパラメータ推定

Stanでは、exp_mod_normal()関数で指数-ガウス分布を扱うことができます。

#stan modelの定義
model_exGauss <- "
data{
  int N;
  real Y[N];
}
parameters{
  real mu;
  real<lower=0> lambda;
  real<lower=0> sigma;
}
model{
  for(i in 1 : N)
  Y[i] ~ exp_mod_normal(mu,sigma,lambda);
}

"

#data to be passed on stan
datastan = list(N = length(df_exGauss$exGauss),Y = df_exGauss$exGauss)

#stan fit
fit_exGauss<-stan(model_code = model_exGauss, data = datastan, seed = 123, iter = 2000, warmup = 1000)

結果の確認

比較的短時間で推定が完了したのではないかと思います。

早速結果を確認してみましょう。

新しいモデルを試すときは、乱数の初期値と推定値が一致しているか、

確かめておくと安心ですね。

#推定値
formattable((summary(fit_exGauss)$summary)[,c(1,4,8,9,10)]%>%data.frame())

おりゃっ

exGaussTable.png

95%確信区間(信用区間)内に、乱数生成に用いた真値が含まれていますね。よかったよかった。

<2017/10/17追記> あ、sigmaがリカバリできてないですね。これは要検討。改めて、後日記事修正します。

#事後分布
stan_hist(fit_exGauss)

てぇぇぇいぃ

exgaussFigure.png

各パラメータの事後分布です。

指数-ガウス分布の回帰モデル

指数-ガウス分布は、正規分布と指数分布の混合分布です。二つの要素が混合しているため、回帰には2つの方法があります。

まずは、正規分布側の平均を線形モデルで予測する方法。

もう一つは、指数分布側の期待値(1/lambda)に対して、対数リンクを介した非線形モデルを通して予測していく方法です。

今回は、後者のモデルを実装したいと思います。

正規分布のパラメータは、今回固定です。

乱数データの生成

  1. 説明変数Xは正規分布からの適当な乱数を出して標準化
  2. muとsigmaはそれぞれ固定
  3. 切片と傾きの値(真値)を決めて、対数リンクの逆リンク関数(指数変換)を介させる
  4. 指数分布のパラメータに変換後の線形モデルを入れ、ex-Gaussian乱数を生成
#乱数生成
N = 1000 #乱数の数(観測数)
X = c(scale(rnorm(N,15,1))) #説明変数の標準化した正規乱数
beta0 = 1 #切片
beta1 = 0.5 #傾き
sigma = 1 #誤差(標準偏差=分散=1)

#指数分布パラメータに線形モデルを
nu=exp(beta0+beta1*X)

#指数ガウス乱数の生成
Y=rexGAUS(N,mu=0.25,sigma=sigma,nu=nu)

#分布の可視化
df_exGauss_reg<-data.frame(Y=Y,X=X)
ggplot(df_exGauss_reg,aes(x=X,y=Y))+
    geom_point()+theme_bw()

とぅぅぅりゃ

exGaussScatter.png

これが今回のデータ。X軸が大きくなるほど、指数分布の期待値が上昇して幅が広がっていっているのがわかります。

stanモデルと推定の実行

Stanではexp_mod_normal()関数で指数ガウス分布を利用することができますが、指数分布のパラメータlambdaに対してそのまま非線形モデルを定義する方が計算が安定します(平均を通るように1/lambda = exp(beta0+beta1*X)としない)。

これをすると、回帰係数と切片パラメータの推定結果が正負逆転して得られますので、結果報告の際には十分にご留意ください。

結果の確認に使用するため、 generated quantities{}ブロックで、ex-gaussian全体の事後予測分布と、muとsigmaの正規分布の部分を乱数生成しておきます。

model_exGauss_reg<-"
data{
  int N;
  real Y[N];
  vector[N] X;
}
parameters{
  real beta0;
  real beta1;
  real<lower=0> sigma;
  real<lower=0> mu;
}
transformed parameters{
  vector<lower=0>[N] lambda;
  for (i in 1 : N)
    lambda[i] = exp(beta0+beta1*X[i]);
}
model{
  for (i in 1 : N)
    Y[i] ~ exp_mod_normal(mu,sigma,lambda[i]);
}
generated quantities{
  vector[N] pred_Y;
  vector[N] pred_lambda;
  real pred_norm;
  pred_norm = normal_rng(mu,sigma);
  for (i in 1:N){
    pred_Y[i] = exp_mod_normal_rng(mu,sigma,lambda[i]);
    pred_lambda[i] = exponential_rng(lambda[i]);
  }
}

"

datastan = list(N = N, Y=Y, X=c(scale(X)))
fit_exGauss_reg<-stan(model_code = model_exGauss_reg, data = datastan, seed = 123, iter = 2000, warmup = 1000)

結果の確認

これも、比較的短時間で推定が完了したのではないかと思います。

切片と傾きパラメータのパラメータリカバリは、うまくできているでしょうか。

確認してみましょう。

#推定値
formattable((summary(fit_exGauss_reg)$summary)[c("beta0","beta1","sigma","mu"),c(1,4,8,9,10)]%>%data.frame())

どっせい!

exGaussRegTable.png

95%区間内に真値が入ってますね。良い感じです。

#事後分布
stan_hist(fit_exGauss_reg,pars = c("beta0","beta1","sigma","mu"))

うりゃうりゃうりゃー

exGaussRegFigure.png

指数-ガウス分布を使った回帰のパラメータ事後分布です。

推定結果の可視化①

まずは、ex-gaussian全体の事後予測分布を可視化したいと思います。

平均のトレンドを緑の太線で表現し、帯を95%区間に設定しています。

#sampleの取り出しと整理
#事後予測分布の95%区間の計算
df_exGauss_res<-rstan::extract(fit_exGauss_reg)$pred_Y%>%
    data.frame()%>%
    gather()%>%
    dplyr::mutate(id = rep(c(1:N),each=4000))%>%
    group_by(id)%>%
    dplyr::summarize(pred_EAP = mean(value),
                     pred_lower = quantile(value,0.025),
                     pred_upper = quantile(value,0.975))%>%
    ungroup()%>%
    cbind(df_exGauss_reg)

#結果の可視化
ggplot(df_exGauss_res)+
    geom_point(aes(x=X,y=Y),alpha=0.95)+
    geom_ribbon(aes(x=X,ymin=pred_lower,ymax = pred_upper),alpha=0.1,fill="darkgreen")+
    geom_line(aes(x=X,y=pred_EAP),colour = "darkgreen",alpha=0.8,size=1.3)+
    theme_bw()

ぷろっとぉぉぉぉぉぉぉぉぉぉ

ExGaussPreface2.png

推定結果の可視化②

続いて、lambdaパラメータで生成した指数分布の乱数と、muとsigmaで生成した正規分布の乱数を取り出して、 その95%区間と共に可視化してみます。

df_exGauss_res2<-rstan::extract(fit_exGauss_reg)$pred_lambda%>%
    data.frame()%>%
    gather()%>%
    dplyr::mutate(id = rep(c(1:N),each=4000))%>%
    group_by(id)%>%
    dplyr::summarize(lambda_EAP = mean(value),
                     lambda_lower = quantile(value,0.025),
                     lambda_upper = quantile(value,0.975))%>%
    ungroup()%>%
    cbind(df_exGauss_res)

#muの取り出し
mu<-mean(rstan::extract(fit_exGauss_reg)$mu)
pred_norm_up<-quantile(rstan::extract(fit_exGauss_reg)$pred_norm,0.975)
pred_norm_lo<-quantile(rstan::extract(fit_exGauss_reg)$pred_norm,0.025)

ggplot(df_exGauss_res2)+
    geom_hline(yintercept = mu,colour = "darkred",size=1.3)+
    geom_ribbon(aes(x=X,ymin=pred_norm_lo,ymax = pred_norm_up),alpha=0.2,fill="darkred")+
    geom_point(aes(x=X,y=Y),alpha=0.95)+
    geom_ribbon(aes(x=X,ymin=lambda_lower,ymax = lambda_upper),alpha=0.1,fill="darkblue")+
    geom_line(aes(x=X,y=lambda_EAP),colour = "darkblue",alpha=0.8,size=1.3)+
    theme_bw()+
    scale_x_continuous(expand=c(0,0))+
    scale_y_continuous(expand=c(0,0))

ぴろろろろろ〜

ExGaaussRegPred.png

指数分布の部分と、正規分布の部分を別々に出力させました。これの和が、ex-gaussianというわけですね。

まとめ

指数-ガウス分布を使ったパラメータ推定と、「指数分布」の部分を対数リンクで回帰する方法をご紹介しました。

「どういうデータなのか」に依存する部分はもちろんありますが、外れ値(右の裾に効く部分)のトレンドを表現することができました。

以下は一例ですが、外れ値こそが重要な意味をもつ場面、あると思います。

たとえば、「電車の車掌さんが、信号を見てブレーキをかけるまでの時間」に外れ値が増えてくると、それって大変なことじゃないでしょうか。考えてみれば、ブレーキまでの「平均反応時間が遅くなる」ことよりも、「外れ値が増える」ことのほうが、危険です。

睡眠の質が低いほど、不注意から外れ値の成分が増えたりはしないでしょうか。このような時に、ex-gaussianの回帰モデルが使えるかもしれません。

外れ値の生成メカニズムに指数分布を仮定している点が、妥当かどうかがわかりませんゆえ、

適用場面は限局的かもしれません。

しかし、ex-gaussianを用いた回帰モデル、いくつか使い所がありそうです。

そんな試論でした。

みなさんも、データ分析に役立ててみてください。

Enjoy!!

Written on October 17, 2017