RとStanでニューラルネットワークモデル

Mr.Unadon

2017/12/01


はじめに


Rで確率的プログラミング言語STANが利用できる{rstan}パッケージ。

私は普段、ベイズ統計モデリングを行うために使っています。

今回は、rstanを使って多層ニューラルネット(フルコネクションの中間層2層)に挑戦する機会があったので、

備忘録としてまとめておこうと思います。

kerasやmxnetでやるよりも圧倒的に時間がかかります。圧倒的です。

一方で、係数やニューロンの値を事後分布として求めることが出来ます。テストデータのラベルも事後予測分布で求めます。

モチベーションは、この先に考えている「ネットワークのベイズ統計モデリング」にあります。

認知モデリングと組み合わせれば、神経活動のベイズ統計モデリングなんかができちゃうかもしれませんね。

現実的な時間内で求まる程度の小規模データならば、ですが。

参照元は Stanでニューラルネットの実装に関する議論です。コードがそのまま公開されていたので、一部自分なりに使いやすく変更しています。

機械学習畑でも統計畑の方でも読めるように、言葉遣いをわざと平易なものにしているところがあります(統計よりです)。

よろしければ批判でもポエムでもアドバイスでもなんでも、教えていただければ幸いです。

Twitter(@MrUnadon)

その他のブログ記事


本稿の構成


  1. データの読込と成形
  2. Stanでのニューラルネットワークモデル
  3. 推定の実行と可視化
  4. 二重指数分布を係数の事前分布にした場合
  5. 結語

1. データの読込と成形


irisのSpecies(種類)であるsetosa(0)とversicolor(1)を被予測変数とした、分類判別を目的とします。

予測には花弁(Petal)とがく片(Sepal)の長さ(length)と幅(width)の4つの情報を使います。

library(rstan) #Rstan
library(formattable) #HTML Table by Kun Ren
library(tidyverse) #月がtidy
library(ggmcmc) #MCMCの事後処理
library(stringr) #文字列の処理

#stan parallel
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

#setseed
set.seed(1234)

#irisからsetosaとversicolorを取り出して01に変換し、先頭列に持ってくる。
#予測するための変数は標準化しておく
d <- iris%>%
    dplyr::filter(Species == "setosa" | Species == "versicolor")%>%
    dplyr::mutate(Petal.Length = c(scale(Petal.Length)),
                  Petal.Width = c(scale(Petal.Width)),
                  Sepal.Length = c(scale(Sepal.Length)),
                  Sepal.Width = c(scale(Sepal.Width)),
                  Species = as.integer(Species) -1
                  )%>%
    dplyr::select(Species,everything())

formattable(d[c(1:5,50:55),])
Species Sepal.Length Sepal.Width Petal.Length Petal.Width
1 0 -0.57815327 0.837617378 -1.0079000 -1.0368872
2 0 -0.88982620 -0.206793318 -1.0079000 -1.0368872
3 0 -1.20149912 0.210970961 -1.0768870 -1.0368872
4 0 -1.35733558 0.002088821 -0.9389130 -1.0368872
5 0 -0.73398974 1.046499517 -1.0079000 -1.0368872
50 0 -0.73398974 0.419853100 -1.0079000 -1.0368872
51 1 2.38273950 0.210970961 1.2686709 1.0864313
52 1 1.44772073 0.210970961 1.1306969 1.2633746
53 1 2.22690304 0.002088821 1.4066449 1.2633746
54 1 0.04519257 -1.668968292 0.7857619 0.9094881
55 1 1.60355719 -0.624557596 1.1996839 1.2633746

続いて、学習データとテストデータの分割を行います。

今回は不要ですが、データが変わったときのことも考えダウンサンプリング(0,1のバランスが1:1になるよう、片方のレコードを削る)を行っておきます。

#行をランダムに抜く
randomRows <- function(d, n) {d[sample(nrow(d),n),]}
#0と1のラベルの数が少ない方に、レコード数を合わせる
downSample <- function(df) {
    c1 <- d[d[,1] == 1,]
    c2 <- d[d[,1] == 0,]
    size <- min(nrow(c1), nrow(c2))
    rbind(randomRows(c1,size), randomRows(c2,size))
}
#関数の適用
balanced_d <- downSample(d)

#学習/テストデータの分割
splitdf <- function(df,rate) {
    index <- 1:nrow(df)
    trainindex <- sample(index, trunc(length(index)*rate))
    testset <- df[-trainindex, ]
    trainset <- df[trainindex, ]
    list(trainset=trainset,testset=testset)
}

#90%を学習に用いる(今回は90レコードになる)
splits <- splitdf(balanced_d, rate=0.9)
#学習・テストデータへの分割
train <- splits$trainset
test <- splits$testset

head(train)
##    Species Sepal.Length  Sepal.Width Petal.Length Petal.Width
## 98       1    1.1360478 -0.415675457    0.9927229   0.9094881
## 23       0   -1.3573356  1.046499517   -1.2838480  -1.0368872
## 89       1    0.2010290 -0.206793318    0.8547489   0.9094881
## 58       1   -0.8898262 -1.460086153    0.3028529   0.3786585
## 61       1   -0.7339897 -2.295614709    0.4408269   0.3786585
## 87       1    1.9152301  0.002088821    1.2686709   1.2633746
head(test)
##    Species Sepal.Length Sepal.Width Petal.Length Petal.Width
## 90       1   0.04519257  -1.2512040    0.7857619   0.9094881
## 79       1   0.82437488  -0.4156755    1.1306969   1.2633746
## 60       1  -0.42231681  -0.8334397    0.7167749   1.0864313
## 85       1  -0.11064389  -0.2067933    1.1306969   1.2633746
## 82       1   0.04519257  -1.4600862    0.5788009   0.3786585
## 95       1   0.20102904  -0.8334397    0.9237359   0.9094881

2. Stanでのニューラルネットワークモデル


活性化関数はシグモイド関数です。

重み付けの係数をbetaとし、切片をbiasをします。

下記stanコード内では、2つの関数を定義しています。

一つ目は各層の入力値を標準化する関数(scaled_inv_logit)で、

二つ目は、中間層のニューロンの値を計算し、出力一歩手前のニューロンの値をalphaで返すcalculate_alpha関数です。

Stanによる多層ニューラルネットワークモデル

//NeuralNetwork.stan
functions {
    matrix scaled_inv_logit(matrix X_train) {
        matrix[rows(X_train), cols(X_train)] res;
        for(i in 1:rows(X_train)){
            for(j in 1:cols(X_train)) {
            res[i,j] = inv_logit(X_train[i,j])*2-1;
            }
        }
        return res;
    }
    
    vector calculate_alpha(matrix X_train, vector bias, matrix beta_first, matrix[] beta_middle, vector beta_output) {
        int N_train;
        int num_nodes;
        int num_layers;
        matrix[rows(X_train),rows(beta_first)] layer_values[rows(bias)];
        vector[rows(X_train)] alpha;
        N_train = rows(X_train);
        num_nodes = rows(beta_first);
        num_layers = rows(bias);
        layer_values[1] = scaled_inv_logit(bias[1] + X_train * beta_first');
        
        for(i in 2:(num_layers-1)) 
            layer_values[i] = scaled_inv_logit(bias[i] + layer_values[i-1] * beta_middle[i-1]');
        alpha = bias[num_layers] + layer_values[num_layers-1] * beta_output;
        return alpha;
    }
}
data {
    int<lower=0> N_train;
    int<lower=0> V;
    int<lower=0> num_nodes;
    int<lower=1> num_middle_layers;
    matrix[N_train,V] X_train;
    int Y_train[N_train];
    int<lower=0> N_test;
    matrix[N_test,V] X_test;
}
transformed data {
    int num_layers;
    num_layers = num_middle_layers + 2;
}
parameters {
    vector[num_layers] bias;
    matrix[num_nodes,V] beta_first;
    matrix[num_nodes,num_nodes] beta_middle[num_middle_layers];
    vector[num_nodes] beta_output;
}

transformed parameters{
    vector[N_train] alpha;
        alpha = calculate_alpha(X_train, bias, beta_first, beta_middle, beta_output);
}

model{
    Y_train ~ bernoulli_logit(alpha);
    
    //priors
    bias ~ normal(0,1);
    to_vector(beta_output) ~ normal(0,1);
    to_vector(beta_first) ~ normal(0,1);
    for(i in 1:(num_middle_layers)) 
        to_vector(beta_middle[i]) ~ normal(0,1);
}
generated quantities{
    vector[N_test] predictions;
    {
        vector[N_test] pred_alpha;
        pred_alpha = calculate_alpha(X_test, bias, beta_first, beta_middle, beta_output);
        for(i in 1:N_test) 
            predictions[i] = inv_logit(pred_alpha[i]);
    }
}

【dataブロック】

  • N_train: トレーニングデータセットの行数(今回は90)

  • V: 予測変数の列数です(今回は4)

  • num_nodes: 中間層各層のニューロンの数です(今回は10)

  • num_middle_layers: 中間層の層の数です(今回は2)

  • X_train: 学習データの予測変数

  • Y_train: 学習データの被予測変数(ラベル)

  • N_test; テストデータの行数(今回は10)

  • X_test; テストデータの説明変数

【parameters ブロック】

  • bias 層の数だけ存在する切片

  • beta_first 入力層から第一中間層への重み付け係数

  • beta_middle 中間層から出力層までの各層への重み付け係数

  • beta_output 出力層の重み付け係数

【modelブロック】

  • 係数の事前分布に標準正規分布を仮定

【generated quantitiesブロック】

  • テストデータの非予測変数に対する事後予想分布を生成

3. 推定(学習)の実行と可視化


stanによる推定を行います。

推定するのは各係数の値。

テストデータの予測判別を事後予測分布で出力します。

中間層のニューロンの値は残しておりませんが、出力層のニューロンの値はalphaで取り出せます。

Y_test = test[,1] #テストデータの被予測変数
N_train = nrow(train) #トレーニングデータの行数
N_test = nrow(test) #テストデータの行数
V = ncol(train)-1 #予測変数の数

datastan=list(
    num_nodes=10, #ニューロン数は10に固定 
    num_middle_layers=2, #中間層は2層
    V=V,  #予測変数の数
    N_train=N_train,  #トレーニングデータのレコード数
    N_test=N_test,  #テストデータのレコード数
    X_train=train[,2:5],  #トレーニングデータの予測変数
    Y_train=train[,1],   #トレーニングデータの被予測変数(ラベル)
    X_test=test[,2:5] #テストデータの予測変数
    ) 

#stanモデルのc++コンパイルと保存
    #nnet_model <- stan_model("NeuralNetwork.stan")
    #saveRDS(nnet_model, "NeuralNet.rds")

#c++コンパイル済みモデルの読込
nnet_model<-readRDS("NeuralNet.rds")


#Stanの推定条件
iter = 2000 #サンプリング総数
warmup = 1000 #warmupによる初期サンプルの棄却数
chains = 4 #MCMC系列の数

#Stanによるパラメータ推定の実行
nn_fit <- sampling(nnet_model,
                   data = datastan,
                   iter = iter,
                   warmup = warmup,
                   chains = chains,
                   seed = 1234)

推定値の確認

推定値の要約量を確認します。

長いので、本稿では記事面に出力しません。

収束は事前に確認しております。

#結果の要約値取り出し
Summary<-summary(nn_fit)$summary[,c(1,3,4,8,9,10)]

#結果をデータフレームに
df<-data.frame(Summary)
#RhatのたまにあるNaNは数値認識されないのでゼロで置き換え
df$Rhat[is.nan(df$Rhat)]<-0
#小数点第3ケタまでにする
for(i in 1 : length(df)){df[,i]<-digits(df[,i],3)}
#lp___を取り除く
df<-df[-length(df$mean),]

#結果を表で描画
res_table<-formattable(df, list(
    mean = formatter("span",style = x ~ style(color=ifelse(x <0, "steelblue","tomato"))),
    sd = color_bar("orange"),
    X2.5. = formatter("span",style = x ~ ifelse(x <0, style(color = "steelblue", font.weight = "bold"),NA)),
    X97.5. = formatter("span",style = x ~ ifelse(x <0, style(color = "steelblue", font.weight = "bold"),NA)),
    n_eff = formatter("span",style = x ~ style(color = ifelse(rank(-x) <= c(length(df$n_eff)-1), "gray", "tomato")),
                      x ~ sprintf("%.2f (rank: %g)", x, rank(-x))),
    Rhat = formatter("span", x ~digits(x,2),
                     style = x ~ style(color = ifelse(x>=1.1, "tomato", "green")),
                     x ~ icontext(ifelse(x<1.1, "ok", "remove"), ifelse(x<1.1, "Yes", "No")))
))

#推定結果の出力
    #res_table

係数の可視化

全てのサンプルを取り出して全部可視化すると重たいので、

全係数の事後分布ついては、お手元でご確認いただければと思います。

ここでは出力層への係数だけを見ておきましょう。他の係数も、分布の位置と形状はほぼ同じような感じです。

beta_first<-stan_hist(nn_fit, pars = c("beta_first"))
beta_middle<-stan_hist(nn_fit, pars = c("beta_middle"))
beta_output<-stan_hist(nn_fit, pars = c("beta_output"))

#beta_first   #入力層から中間層第一層への係数
#beta_middle  #中間層から出力層まえまでの係数
beta_output  #出力層への係数

サンプルの取り出し、判別推定結果の可視化

判別結果を事後予測分布で見てみます。

samples<-ggs(nn_fit)%>%
    dplyr::filter(str_detect(Parameter, "predictions"))%>%
    dplyr::arrange(Parameter)%>%
    dplyr::mutate(Y_test = rep(Y_test,each = c((iter-warmup) * chains)))

g1<-samples%>%
    ggplot(aes(x = Parameter, y = value))+
    geom_violin(fill = "darkgreen", alpha = 0.5)+
    geom_point(aes(x = Parameter, y = Y_test))+
    theme_bw()
g1

なんか、それっぽい。


4. 二重指数分布を係数の事前分布とした場合

せっかくStanでやっているので、ベイジアンらしく事前分布を変えてみます。

この辺に事前知識を与えられるのは、なんだかワクワクします。

ここでは、適当に係数の事前分布に二重指数分布(Laplace分布)を置いてみます。

コストパラメータは適当に所与しています。理論とかありません。適当にやっても推定がうまくいくか試す目的です。

biasと係数の事前分布にLaplace分布を仮定したStanモデル

//NeuralNetwork_Laplace.stan
functions {
    matrix scaled_inv_logit(matrix X_train) {
        matrix[rows(X_train), cols(X_train)] res;
        for(i in 1:rows(X_train)){
            for(j in 1:cols(X_train)) {
            res[i,j] = inv_logit(X_train[i,j])*2-1;
            }
        }
        return res;
    }
    
    vector calculate_alpha(matrix X_train, vector bias, matrix beta_first, matrix[] beta_middle, vector beta_output) {
        int N_train;
        int num_nodes;
        int num_layers;
        matrix[rows(X_train),rows(beta_first)] layer_values[rows(bias)];
        vector[rows(X_train)] alpha;
        N_train = rows(X_train);
        num_nodes = rows(beta_first);
        num_layers = rows(bias);
        layer_values[1] = scaled_inv_logit(bias[1] + X_train * beta_first');
        
        for(i in 2:(num_layers-1)) 
            layer_values[i] = scaled_inv_logit(bias[i] + layer_values[i-1] * beta_middle[i-1]');
        alpha = bias[num_layers] + layer_values[num_layers-1] * beta_output;
        return alpha;
    }
}
data {
    int<lower=0> N_train;
    int<lower=0> V;
    int<lower=0> num_nodes;
    int<lower=1> num_middle_layers;
    matrix[N_train,V] X_train;
    int Y_train[N_train];
    int<lower=0> N_test;
    matrix[N_test,V] X_test;
}
transformed data {
    int num_layers;
    num_layers = num_middle_layers + 2;
}
parameters {
    vector[num_layers] bias;
    matrix[num_nodes,V] beta_first;
    matrix[num_nodes,num_nodes] beta_middle[num_middle_layers];
    vector[num_nodes] beta_output;
}

transformed parameters{
    vector[N_train] alpha;
        alpha = calculate_alpha(X_train, bias, beta_first, beta_middle, beta_output);
}

model{
    Y_train ~ bernoulli_logit(alpha);
    
    //priors
    bias ~ double_exponential(0,0.3);
    to_vector(beta_output)  ~ double_exponential(0,0.3);
    to_vector(beta_first)  ~ double_exponential(0,0.3);
    for(i in 1:(num_middle_layers)) 
        to_vector(beta_middle[i])  ~ double_exponential(0,0.3);
}
generated quantities{
    vector[N_test] predictions;
    {
        vector[N_test] pred_alpha;
        pred_alpha = calculate_alpha(X_test, bias, beta_first, beta_middle, beta_output);
        for(i in 1:N_test) 
            predictions[i] = inv_logit(pred_alpha[i]);
    }
}

推定の実行と可視化

上記と同様に推定を実行して推定値を確認します。

#モデルのコンパイルと保存
    #nnet_laplace_model <- stan_model("NeuralNetwork_Laplace.stan")
    #saveRDS(nnet_laplace_model, "NeuralNet_laplace.rds")

#コンパイル済みモデルの読込
nnet_laplace_model<-readRDS("NeuralNet_laplace.rds")

#stan fit
nn_laplace_fit <- sampling(nnet_laplace_model,
                   data = datastan,
                   iter = iter,
                   warmup = warmup,
                   chains = chains,
                   seed = 1234)



#結果要約値の取り出し
Summary_laplace<-summary(nn_laplace_fit)$summary[,c(1,3,4,8,9,10)]

#結果をデータフレームに
df_laplace<-data.frame(Summary_laplace)
#RhatのたまにあるNaNは数値認識されないのでゼロで置き換え
df_laplace$Rhat[is.nan(df_laplace$Rhat)]<-0
#小数点第3ケタまでにする
for(i in 1 : length(df_laplace)){df_laplace[,i]<-digits(df_laplace[,i],3)}
#lp___を取り除く
df_laplace<-df_laplace[-length(df_laplace$mean),]

#結果を表で描画
res_table_laplace<-formattable(df_laplace, list(
    mean = formatter("span",style = x ~ style(color=ifelse(x <0, "steelblue","tomato"))),
    sd = color_bar("orange"),
    X2.5. = formatter("span",style = x ~ ifelse(x <0, style(color = "steelblue", font.weight = "bold"),NA)),
    X97.5. = formatter("span",style = x ~ ifelse(x <0, style(color = "steelblue", font.weight = "bold"),NA)),
    n_eff = formatter("span",style = x ~ style(color = ifelse(rank(-x) <= c(length(df$n_eff)-1), "gray", "tomato")),
                      x ~ sprintf("%.2f (rank: %g)", x, rank(-x))),
    Rhat = formatter("span", x ~digits(x,2),
                     style = x ~ style(color = ifelse(x>=1.1, "tomato", "green")),
                     x ~ icontext(ifelse(x<1.1, "ok", "remove"), ifelse(x<1.1, "Yes", "No")))
))

#結果の確認
    #res_table_laplace

係数の事後分布を見ておきます。 重いのでoutputに関わる係数だけ確認します。

beta_first_laplace<-stan_hist(nn_laplace_fit, pars = c("beta_first"))
beta_middle_laplace<-stan_hist(nn_laplace_fit, pars = c("beta_middle"))
beta_output_laplace<-stan_hist(nn_laplace_fit, pars = c("beta_output"))


#beta_first_laplace
#beta_middle_laplace
beta_output_laplace

テストデータに対する予測、事後予測分布です。

samples_laplace<-ggs(nn_laplace_fit)%>%
    dplyr::filter(str_detect(Parameter, "predictions"))%>%
    dplyr::arrange(Parameter)%>%
    dplyr::mutate(Y_test = rep(Y_test,each = c((iter-warmup) * chains)))

g2<-samples_laplace%>%
    ggplot(aes(x = Parameter, y = value))+
    geom_violin(fill = "darkblue", alpha = 0.5)+
    geom_point(aes(x = Parameter, y = Y_test))+
    theme_bw()

g2

こっちもそれっぽい。


5. 結語


係数の値に「0より大きい」とか「ゼロより小さい」などの制約を置けば、

興奮性と抑制性の神経ネットワークのモデリングができるかもしれません(わかりません)。

しかし、irisを分けるのに5分とかは辛いです。

もっと大規模なデータで、

ニューロン数や層の数、事前分布の値などのチューニングとかを考え始めると、コンピュータお姉さんになりそうです。

どこかに、使い所があると良いのですが。

Enjoy!! Stan!!