一日一膳(当社比)

RとJavaと時々数学

【R】listにlistを要素として追加しようとして1時間溶かした

R 小ネタ

Rでlistにlistを追加しようとしただけなのに1時間溶かしてしまったので、始末書がてら記事化します。

1. 行いたかったこと

つぎのようなlist l1, l2があるとします。


l1 <- list(
    list(x=1, y=2),
    list(x=3, y=4)
)

l2 <- list(x=5, y=6)

これらを結合してつぎのようなlistが作りたい。


list(
    list(x=1, y=2),
    list(x=3, y=4),
    list(x=5, y=6)
)

listの追加ならappendでうまくいくはず。

append(l1, l2)

が、ダメっ...！ l2は展開されてしまうという結果に。


list(
    list(x=1, y=2),
    list(x=3, y=4),
    x=5,
    y=6
)

2. 解決法

a. nativeで書く方法

list1に要素を仮置きで追加したlistを作成し、そのあとl2へ参照させれば欲しいlistが得られます。


add1 <- function(l, new) {
    tmp <- append(l, 0)
    tmp[[length(l) + 1]] <- new
    return(tmp)
}

add1(l1, l2)
b. rlang::execを使う

rlangパッケージのexec関数を使うと、 listを複数の引数に展開して関数に適用させることができます。


require(rlang)
add2 <- function(l, new) {
    return(rlang::exec(list, !!!l1, l2))
}

add2(l1, l2)

なんか自分のブログのアクセス数が木曜日に増えるっぽい

 

自分のブログのアクセス数を見ていて、木曜日に若干アクセスが伸びている気がした。

そのことについての雑感を書きます。

 

木曜日にアクセス数が増えてるっぽい

まず昨年の11/19日から今年1/17までの60日間のアクセス数をグラフにしてみた。

f:id:kimigayoseishou:20200119232359p:plain

アクセス数推移

木曜日を青色でハイライトしてみると、木曜日のアクセス数が多い気(?)がしてきますね。そこで平均アクセス数を木曜日 / 木曜日以外で集計してみると、

曜日 平均アクセス数
木曜日以外 4.9
木曜日 8.6
全体 5.5

と、木曜日が全体と比較して5割強程度アクセスが多いようでした。

状態空間モデルで曜日別アクセス数を推定

統計モデルを使って曜日ごとにアクセス数に違いがあるか検証してみます。

{\sf access}_tを t日目のアクセス数とします。

曜日ごとのアクセス数の違い s_tと、定常部分 w_tを導入して、

つぎのような状態空間モデルを用いてモデル化します。

 

w_t \sim {\sf Normal}(\mu, \sigma_w ^ 2), \quad ({\rm i.i.d.})
s_t \sim {\sf Normal} \bigl(-\sum_{j = 1}^6 s_{t-j}, \sigma_s^2\bigr),
{\sf access}_t \sim {\sf PoissonLog}(\lambda_t),\quad \lambda_t = w_t + s_t,

 

ただし、 {\sf PoissonLog}は対数Poisson分布でつぎの確率質量分布をもつ離散確率変数

  {\sf PoissonLog}(x | \lambda ) = \frac{1}{x!} \exp(\lambda x - \exp(\lambda)), \quad (x = 0, 1, 2, \ldots)

です。

このモデルで \exp(s_t)を推定し曜日ごとについてバイオリンplotにしてみた図が以下になります(横線は四分位点)。木曜日(day_id = 4)で中央値が約1.6と、木曜日のアクセス数が他の曜日より多いことがわかります。

 

f:id:kimigayoseishou:20200120000134p:plain

結局木曜日にアクセス数が多い理由は？

f:id:kimigayoseishou:20200120001957p:plain

アメリカからのアクセスが木曜日に増えてるっぽい

Google Analyticsで調べてみたら、どうやらアメリカからのアクセスが木曜日に増えてるみたいだ。しかもそのアクセスのほとんどがイリノイ州からによるものでした。

f:id:kimigayoseishou:20200120003933p:plain

USからのアクセスでは、そのほとんどがイリノイ州からのアクセス。

 

不思議に思って調べてみるとつぎのような記事を発見。

www.dj-mope.com

引用した記事によると、イリノイ州からアクセスしているユーザーはmicrosoftのbotによるアクセスなのだとのこと。ここから推測するに、自分の場合もこの方と同様に

botによるアクセスが発生しており、そのスケジューリング設定(?)の影響で木曜日にアクセスされているのではと思うなど。。。

 

以上、取り止めのない内容ですが今日はここで失礼いたします。

【R】WAICを用いた予測モデル評価を混合Poisson分布で試す

R Stan 統計学

こちらの本の6, 7章の汎化損失・WAICの内容をまとめてみました。

社会科学のための ベイズ統計モデリング (統計ライブラリー　)

社会科学のための ベイズ統計モデリング (統計ライブラリー　)

 

 理論と数値計算(R + RStan)がバランスよく書かれているのでおすすめです。

予測モデルの評価と汎化損失

真の分布q(X)が未知であり、真の分布からの独立なn個のサンプルの実現値 \underline{x} = (x_1, \ldots, x_n)が与えられた状況を考えます。この実現値を用いて、予測モデル p^\ast (X)を作成したとき予測パフォーマンスのよさを測る指標としてつぎの汎化損失があります。

汎化損失
G_n := {\mathbb E}_{q(X)}[ - \log p^\ast (X)] = - \int q(X)\log p^\ast (X) dX

真の分布が未知の場合、汎化損失の計算は直接計算が行えないことに注意が必要です。

以下、予測モデルとしてベイズモデルを使う状況を考えます。

  • モデル(尤度) : p(\underline{x} | \theta)
  • 事前分布 :  \phi(\underline{x})

このモデルの予測分布(predictive distribution)は、つぎの期待値で定義されます。

p^\ast (X) := {\mathbb E}_{p(\theta |\underline{x})}[p(X |\theta ] = \int p(X | \theta) p(\theta |\underline{x}) d\theta

特に、ベイズ予測モデルについての汎化損失はつぎの式で表せることになります。

G_n = - \int q(X) \log \Bigl( \int p(X | \theta) p(\theta |\underline{x}) d\theta \Bigr) dX

WAIC

ここではベイズ予測モデルの汎化損失について、真の分布が未知の場合にも計算が可能な推定量となるWAICを導入します。

WAIC
\underline{x} = (x_1, \ldots, x_n)を真の分布のサンプルの実現値とする。
経験損失(training loss)T_nと汎関数分散(general function variance)を以下で定義する
T_n := - \frac{1}{n} \sum_{i = 1}^n {\mathbb E}_{p(\theta | \underline{x})} [\log p(x_i|\theta)]
V_n = \sum_{i = 1}^n \{ {\mathbb E}_{p(\theta | \underline{x})} [ (\log p(x_i|\theta))^2 ] - {\mathbb E}_{p(\theta | \underline{x})} [ \log p(x_i|\theta) ]^2 \}.
このとき、WAIC(Watanabe-Akaike Information Criterion)は次で定義される
{\sf WAIC} := T_n + \frac{V_n}{n}

定義からわかるようにWAICはサンプルの実現値の出方によってばらつきがあります。しかし、つぎの結果によりサンプルの実現値に関する平均をとったときには、WAICは近似的に汎化損失の推定量を与えることが知られています。

q(\underline{X}) \sim \prod_i^n q(X_i)を真の分布のサンプルとする。このとき次が成り立つ。
{\mathbb E}_{q(\underline{X})}[{\sf G}_n] = {\mathbb E}_{q(\underline{X})}[ {\sf WAIC}_n ] + o(1/n)
(ただし、 o(-)はランダウの記号を表す。)

混合Poissonモデルで数値実験

以下混合Poisson分布モデルを題材に、WAICを用いた予測モデル評価をRを使って行います。

真の分布として、2つのPoisson分布 {\sf Poisson}(\lambda_i), i = 1, 2を比率 \mu, ( \mu \in (0, 1))で混合して得られる分布 q(X)を考えます。このとき確率質量関数は

q(X) = \mu {\sf Poisson}(X | \lambda_1) + (1 - \mu) {\sf Poisson}(X | \lambda_2)

となります。以下 \mu = 0.5, \lambda_1=2, \lambda_2=8としてシミュレーションを行います。

まず、真の分布から適当な数のサンプルを作成します。


# 真の分布のパラメーター
target.param <- list(
    mu = 0.5,
    lambda = c(2, 8)
)

# 真の分布からsampler
MixPois.Sampler <- function(param) {
    return(function(n) {
        choices <- rbinom(n, 1, param$mu)
        return(mapply(function(choice, x, y) {
            return(ifelse(choice == 1, x, y))},
            choices, rpois(n, param$lambda[1]), rpois(n, param$lambda[2]))
        )
    })
}

# サンプリング
N <- 200L
target.Sampler <- MixPois.Sampler(target.param)
target.sample <- target.Sampler(N)
hist(target.sample)

f:id:kimigayoseishou:20200101122444j:plain

真の分布からのサンプル(N = 200)

このサンプルに対して、混合Poisson分布を用いたモデリングを行います。

N = 200をサンプル数, \underline{x} = (x_1, \ldots, x_N)サンプルとして

K > 1個のクラスターをもつ混合Poisson分布の尤度は

p(\underline{x} | \underline{\lambda}, \underline{\mu}) = \prod_{n =1}^N \sum_{k = 1}^K \mu_k {\sf Poisson}(x_n | \lambda_k)

で定義されます。ただし、\underline{\lambda}, \underline{\mu}は以下のようなハイパーパラメーターとします。

  • \lambda_k \sim {\sf Gamma}(a, b),\quad a, b \in (0, \infty) , k = 1, \ldots, K
  • \underline{\mu} \sim {\sf Dirichlet}(\underline{r}), \quad \underline{r}=(r_1, \ldots, r_K) \in {\mathbb R}_{+}^K

以下、この予測モデルについてクラスター数 Kを変化させたときの予測精度の変化をWAICで評価してみます。

WAICを求める上で事後分布 p(\underline{\lambda}, \underline{\mu} |\underline{x})が必要ですが、この分布は解析的に解くことが難しいためRStanを使ったモンテカルロ法で処理します。

Kが2以上の場合

 以下stanコードになります。


// Reference
// http://rpubs.com/kaz_yos/fmm2
// https://mc-stan.org/docs/2_21/stan-users-guide/soft-k-means.html
data {
  int<lower=1> K;//number of Cluster 
  int<lower=1> N;//data size
  int<lower=0> x[N];//sample 
  real<lower=0> a;
  real<lower=0> b;
  vector[K] r;//Dirichlet parameter
}

parameters {
  positive_ordered[K] lambda;
  simplex[K] mu;
}

transformed parameters {
  real lp[N, K];
  for(n in 1:N) {
    for(k in 1:K) {
      lp[n, k] = log(mu[k]) + poisson_lpmf(x[n]|lambda[k]); 
    }
  }
}

model {
  //prior 
  mu ~ dirichlet(r);
  
  for(k in 1:K) {
    lambda[k] ~ gamma(a, b);
  }
  
  //model
  for(n in 1:N) {
    target+= log_sum_exp(lp[n]);
  }
}

generated quantities {
  real lp_sample[N, K];
  real log_lik[N];
  for(n in 1:N) {
    for(k in 1:K) {
      lp_sample[n, k] = log(mu[k]) + poisson_lpmf(x[n]|lambda[k]); 
    }
    log_lik[n] = log_sum_exp(lp_sample[n]);
  }
}

 

このコードをRからキックして、事後分布からの乱数が得られます。

以下は K = 2の場合のRコードになります。


.model = rstan::stan_model(file = './mixture_of_poisson.stan')
model.data <- list(
    K = 2L,# クラスター数
    N = N,# サンプルサイズ
    x = target.sample,# サンプル
    a = 1, b = 1,
    r = c(1, 1) # K-vector 
)

# 各chain 試行回数2000 & warmup = 1000で4回サンプリング
iter <- 4000L
.fit = rstan::sampling(.model, data = model.data, iter=iter, chains=4)

# Rhat < 1.1に収まっているか収束確認
summary <- rstan::summary(.fit)$summary
summary[summary[,'Rhat'] >= 1.1,] %>% View

# 収束を目視で確認
.pars <- c(str_glue('lambda[{c(1:model.data$K)}]'),
           str_glue('mu[{c(1:model.data$K)}]'))
rstan::stan_hist(.fit, pars = .pars)

 

WAICは発生させたサンプルからつぎのように計算できます。


waic_mcmc <- function(log_lik) {
    train_loss <- mean( - log(colMeans(exp(log_lik))))
    gen_fun_var.mean <- mean(apply(log_lik, MARGIN = 2, var))
    waic <- train_loss + gen_fun_var.mean
    return(waic)
}
model.sample <- rstan::extract(.fit)
model.waic_mcmc <- waic_mcmc(model.sample$log_lik)

 

K=1の場合

この場合はモデルは単にPoissonモデルになりますので、stanコードは以下のようになります。 サンプリング & WAICのシミュレーションは K > 1の場合と同様なので省略します。


// Reference
data {
  int<lower=1> N;//data size
  int<lower=0> x[N];//sample 
  real<lower=0> a;
  real<lower=0> b;
}

parameters {
  real<lower=0> lambda;
}

model {
  //prior 
  lambda ~ gamma(a, b);
  
  //model
  target+= poisson_lpmf(x | lambda);
}

generated quantities {
  real log_lik[N];
  for(n in 1:N) {
    log_lik[n] = poisson_lpmf(x[n]|lambda);
  }
}

 

予測分布のシミュレーション

\{\theta_s \}_{s = 1, ... S}を事後分布 p(\theta | \underline{x})からのサンプルとしたとき、モンテカルロ法により予測分布はつぎのように近似的に求められます

p(X | \underline{x}) = \int p(X | \theta ) p(\theta | D) d\theta \approx \frac{1}{S} \sum_{s = 1}^S p(X | \theta_s).

今回の数値実験ではR側で事後分布の計算を行うことにします。


model.pred_pmf <- function(x, sample) {
    return(x %>% sapply(function(.x) {
        if (is.null(model.sample$mu)) {
            # K = 1の場合
            model.sample$mu <- rep(1, length(x))
        }
        pred_sample <- model.sample$mu %*% t(dpois(.x, model.sample$lambda))
        return(mean(pred_sample))
    }))
}
# 事後分布p(X | \underline{x})をX = 0, ..., 20についてシミュレーション
pred_pmf <- model.pred_pmf(c(0:20L), model.sample)

K = 1~4の各モデルについて予測分布をシミュレーションした結果がつぎの図になります。図を見ると、 K=2, 3, 4のモデルでは真の分布のサンプルに見られた2つのピークをうまく表現できていのに対し、 K=1のモデルではピークをうまく表現できていないことがわかります。

f:id:kimigayoseishou:20200102201445p:plain

各モデルの予測分布と真の分布からのサンプルの度数の比較

 

WAICのシミュレーション結果

K=1, 2, 3, 4でWAICのシミュレーションした結果が以下の表になります。

K WAIC
1 2.98
2 2.60
3 2.60
4 2.60

 

この表から以下のような観察が可能です。

  • K=1のときのモデル(Poisson分布)よりも、(真の分布のクラスター数である) K=2のモデルの方がWAICが小さく、予測のパフォーマンスが高い。
  • K= 2, 3, 4のモデル間ではWAICの変化が見られず、クラスター数を増やすことによる予測パフォーマンスの改善は(K < 5 の範囲では)見られない。

Macでrstanの実行環境を用意する

Stan R

 

パソコン買い換えた時にrstanの実行環境を用意した時のメモです。

想定するOS/Rのバージョン

  • OS　Catalina 10.15.1
  • R 3.6.1 

手順1.C++コンパイラの導入

Stanの実行に必要なのでc++のコンパイラをPCにいれてない人はこの手順1を行ってください。(他にXcodeというMac向け開発環境からc++を入れる方法もあるようなのですが、自分は失敗しました。)

CRANのサイトからclang-7.0.0.pkgをダウンロードしてインストールしてください。

cran.r-project.org

以上で手順1は完了です。

手順2.rstanのダウンロード

RStan公式のgitリポジトリに従ってrstanのインストールを行います。

github.com

古いバージョンのrstanを消します。RStudioのコンソールから次のコマンドを実行しましょう。

remove.packages("rstan")
if (file.exists(".RData")) file.remove(".RData")

 

最後にRstanのパッケージをダウンロードします。

install.packages("rstan", repos = "https://cloud.r-project.org/", dependencies = TRUE)

以上で、rstan環境が用意完了です。試しにpoisson分布のサンプリングでテストしてみます。


library("rstan")
Poisson.model <- "
data {
    int N;
    int<lower=0> x[N];
}
parameters {
    real<lower=0> lambda;
}
model {
    for (n in 1:N) {
        x[n] ~ poisson(lambda);
    }
}    
"
N <- 100
lambda <- 4
x <- rpois(N, lambda)
.data <- list(N = N, x = x)
.fit <- stan(
    model_code = Poisson.model ,
    data=.data,
    iter=2000,
    warmup=200
)
stan_hist(.fit)

f:id:kimigayoseishou:20191116184452j:plain

サンプリング結果


サンプリングできていれば、成功です。もしrstan実行時に次のようなエラーメッセージが吐かれている場合、c++導入が正しくできてない可能性があるので手順1の方法を試してみると良いかもしれません。


    compileCode(f, code, language = language, verbose = verbose) でエラー: Compilation ERROR, 
    function(s)/method(s) not created! /bin/sh: /usr/local/clang7/bin/clang++: No such file or directory

「蜜蜂と遠雷」というキーワードと一緒にツイートされる作曲家をRで集計

R 小ネタ

「蜜蜂と遠雷」に関するネタ分析を行ってみた。

「蜜蜂と遠雷」について

「蜜蜂と遠雷」は恩田陸による長編小説です. ピアノコンクールの出場者を題材とした作品で、 10/4から映画が公開されています. 大まかな内容はというと、ピアノコンクールの出場者を題材(Wikipedia)に、 クラシック音楽(特にピアノ曲)が登場する作品となっているようです。

そのため、「蜜蜂と遠雷」というキーワードを含むツイートの中には、クラシック音楽の作曲家の名前が含まれると推測できます。 今回の記事では、この推測について簡単に検証してみることにします。

分析結果

「蜜蜂と遠雷」を含む直近1万ツイートを取得し(Twitter APIのRラッパー{rtweet}を利用)、 そのツイート本文中に含まれる固有名詞を抽出してみたところ出現回数Top20は以下のようになった.

f:id:kimigayoseishou:20191013191622p:plain
「蜜蜂と遠雷」を含むツイート(直近1万ツイート)」内に含まれる固有名詞トップ20

キーワードを見てみると、次のようなことがわかる.

  • トップは「https」で、Webページリンクが頻繁にツイートされていると思われる.
  • 映画版の出演者の苗字である「松岡」(松岡茉優)と「松坂」(松坂桃李)がそれぞれ2番目と5番目にランク入り.
  • 原作者名に関連した「恩田」と「陸」がそれぞれ3番目と4番目にランク入り.

また、作曲者の名前としてはロシアの作曲家「プロコフィエフ」が16番目にランクインしていた。

その他「蜜蜂と遠雷」というキーワードと一緒にツイートされる作曲家

トップ20には入らなかった作曲家も追加で調べてみた. 作曲家名一覧のデータが欲しかったが、あいにく手持ちの良いデータがなかったため Wikipedia)から拝借した.

抽出結果は次のようになった。 抽出方法があまり良くないため解釈が困難なワードが散見されるが、最も出現回数の多い「プロコフィエフ」の他に、 「ドビュッシー」や「ショパン」などのピアノ曲で有名な作曲家の名前が一緒にツイートされているようだ.

f:id:kimigayoseishou:20191013191839p:plain
「蜜蜂と遠雷」というキーワードと一緒にツイートされた作曲家

Rコード

分析に使ったRコードを貼っておきます. もっと手際の良い方法などありましたら、教えてくださると管理人は泣いて喜びます.


#'R+Twitter
#'rtweet R documentation
#'https://cran.r-project.org/web/packages/rtweet/rtweet.pdf
library(rtweet)
library(tidyverse)
library(RMeCab)
library(rvest)

#tweet取得(日本標準時の列created_at.JST(1時間刻み)を追加している)
data <- search_tweets("蜜蜂と遠雷", n = 10000, include_rts = FALSE) %>%
    dplyr::mutate(
        created_at.JST = as.POSIXct(stringr::str_sub(created_at, 1, -1), 
                                 format="%Y-%m-%d %H", tz = "Japan"))

#ファイル保存
data$text %>% write(file = '蜜蜂と遠雷_tweets_text.txt')

#「蜜蜂と遠雷」を含んだツイートに現れる固有名詞
freq <- data.frame(as.matrix(RMeCabFreq('蜜蜂と遠雷_tweets_text.txt')))
top_words <- freq %>% 
    dplyr::filter(Info1 == '名詞', Info2 == '固有名詞') %>%
    dplyr::arrange(desc(Freq))

bar_plot.Top20 <- top_words[1:20,] %>% 
    ggplot(mapping = aes(x = reorder(Term, as.numeric(Freq)), y = as.numeric(Freq))) + 
    geom_bar(stat = "identity") + 
    xlab('キーワード') + ylab('出現回数') + 
    ggtitle('「蜜蜂と遠雷」を含む直近1万ツイート内に含まれる\n固有名詞(上位20)') + 
    theme_bw(base_family = 'HiraKakuProN-W3') + 
    coord_flip()
print(bar_plot.Top20)

#作曲家一覧をWikipediaのWebページから取得
#https://ja.wikipedia.org/wiki/クラシック音楽の作曲家一覧_(五十音順)
#{ミドルネーム}・{ラストネーム}が抽出できる
composers <- read_html('https://ja.wikipedia.org/wiki/クラシック音楽の作曲家一覧_(五十音順)') %>%
    html_nodes(css = '#mw-content-text > div > ul > li > a') %>%
    html_attr('title')

#作曲家のラストネームを取得
composers_lastName <- composers %>%
    sapply(function(name) unlist(str_split(name, pattern = "・"))[2]) %>%
    as.character()

#「蜜蜂と遠雷」を含んだツイートに現れる作曲家
top_words_composers <- top_words %>% 
    dplyr::filter(Term %in% composers_lastName)

bar_plot.TopComposers <-　top_words_composers %>% 
    ggplot(mapping = aes(x = reorder(Term, as.numeric(Freq)), y = as.numeric(Freq))) + 
    geom_bar(stat = "identity") + 
    xlab('キーワード') + ylab('出現回数') + 
    ggtitle('「蜜蜂と遠雷」を含む直近1万ツイート内に含まれる作曲家') + 
    theme_bw(base_family = 'HiraKakuProN-W3') + 
    coord_flip()
print(bar_plot.TopComposers)

隠れマルコフ連鎖モデルの変分推論(Rで実装)

R 統計学

今読んでいるベイズ推論のテキストベイズ推論による機械学習入門 の内容より、 隠れマルコフモデルの変分推論についてRで実装したので記事化してみます。

1.隠れマルコフモデルにおける変分推論

モデルの確認

K個の状態 {1, 2, ..., K}の間を移りあうマルコフ連鎖を考えます。 この状態が潜在変数で表されているモデルを隠れマルコフモデルといいます。

いくつか記号を用意します。 状態を表す潜在変数 s_{n, k}を時点 nにおける状態がkの場合1, それ以外の場合は0であると定めます。また、 \underline{s_n} = (s _ {n, 1}, ..., s _ {n, K})^ {\sf T}, \quad {\sf S} = (\underline{s_1},..., \underline{s_N}) とまとめて書きます。

2つのハイパーパラメーター

\pi = (\pi_1, ..., \pi_K) , \quad (\sum _ {i = 1} ^ K \pi_k = 1) : \underline{s_1}の分布

{\sf A} = (A _ {i, j}) _ {1 \leq I, j \leq K}, \quad (\sum _ {l = 1} ^ K A _ {l, k}) = 1 : 遷移行列

を与えると、

( \underline{s_n}) _ {1 \leq n \leq N}はマルコフ連鎖に従うので、同時確率 p({\sf S} | \pi, {\sf A})が次の式で与えられます。

p({\sf S} | \pi, {\sf A}) = p( \underline{s _ 1} | \pi ) \prod _ {n = 2} ^ N p( \underline{s _ n} | \underline{s _ {n - 1}}, {\sf A})

また、 \underline{s_n}は直接観測されない代わりに以下の式で与えられる確率変数 x_nが観測データとして与えられるとします。

p(x_n | \underline{s_n}, \underline{\theta}) = \prod _ {k = 1} ^ K {\rm Poisson}(x_n | \theta_k) ^ {s _ {n, k}}

ただし、 \underline{\theta} = (\theta _ 1, ..., \theta _K)はハイパーパラメーターとします。

以上の設定のもと、観測データ {\sf X} = (x_1, ..., x_N)に対する同時分布が次の式で与えられることになります。

p({\sf X}, {\sf S}, \underline{\theta} , \pi, {\sf A})

\quad = p({\sf X} | {\sf S}, \underline{\theta} ) p({\sf S} | \pi, {\sf A}) p(\underline{\theta})p(\pi)p({\sf A})

\quad = p(\underline{\theta})p(\pi)p({\sf A}) p(x_1 | \underline{s _ 1}, \theta _ 1 ) p( \underline{s _ 1} | \pi _ 1) \prod _ {n = 2} ^ N p(x_n |\underline{s _ n} , \underline{\theta}) p( \underline{s _ n} | \underline{s _ {n - 1}}, {\sf A})

完全分解変分推論

隠れマルコフモデルの事後分布 p({\sf S}, \underline{\theta} , \pi, {\sf A} | {\sf X}) を近似的に求める方法として、変分推論を適用することを考えます。

完全分解型変分推論( {\sf completely \ factorized \ variational \ inference})では事後分布を次のような分布の積によって近似を行います。

p({\sf S}, \underline{\theta} , \pi, {\sf A} | {\sf X}) \approx q( \underline{\theta}, {\sf A}, \pi ) \prod _ {n = 1} ^ N q( \underline{s _ n})

テキストでは次の事前分布に対する仮定のもと、完全分解型変分推論に対する更新アルゴリズムの導出が載っていますので気になった方はご参照ください。

\pi \sim {\rm Dirichlet} (\underline{\alpha}), \quad \underline{\alpha} = (\alpha _ 1, ..., \alpha _ K)

{\sf A} \sim \prod _ {k = 1} {\rm Dirichlet} (\underline{ \beta _ k }), \quad \underline{ \beta _ k } = (\beta _ {1, k}, ..., \beta _ {K, k})

\theta _ k \sim {\rm Gamma} (a, b)

2.Rによる実装


library(tidyverse)
library(glue)
#'@param N
#'@param K 
#'@param X array([N])
#'@param prior list(alpha[K], beta[K, K], theta list(a[K], b[K]))
#'@param init list(alpha[K], beta[K, K], theta list(a[K], b[K]))
#'@param iterMax
CFVI <- function(N, K, X, prior, init, iterMax) {

    softmax <- function(s) exp(s) / sum(exp(s))
    entropy.Gamma <- function(a, b) - digamma(a) + log(b)
    entropy.Dirichlet <- function(alpha) - digamma(alpha) + digamma(sum(alpha))
    
    hyp.params = list(
        S = array(rep(0, K * N), dim = c(K, N)),
        pi = init$alpha,
        beta = init$beta,
        theta = list(a = rep(init$theta$a, K), b = rep(init$theta$b, K)))
    
    for (iter in 1:iterMax) {
        #update of hyperparams of S
        hyp.params$S <- c(1:N) %>% 
            map_dfc(function(n) {
                entropy.beta <- hyp.params$beta %>% 
                    as.data.frame %>%
                    map_dfc(~ {entropy.Dirichlet(.)}) %>% 
                    as.matrix
                entropy.theta <- entropy.Gamma(hyp.params$theta$a, hyp.params$theta$b)
                exp.theta <- hyp.params$theta$a / hyp.params$theta$b
                
                if (n == 1) {
                    return(c(- X[1] * entropy.theta - exp.theta - 
                               entropy.Dirichlet(hyp.params$pi) -
                               hyp.params$S[,2] %*% entropy.beta))
                } else if (n == N) {
                    return(c(- X[N] * entropy.theta - exp.theta - 
                                hyp.params$S[,N - 1] %*% entropy.beta))
                } else {
                    return(c(- X[n] * entropy.theta - exp.theta -
                               hyp.params$S[,n - 1] %*% t(entropy.beta) - 
                               hyp.params$S[,n + 1] %*% entropy.beta))
                }
            } %>% softmax) %>% as.matrix
        
        #update of hyperparams of alpha
        hyp.params$pi <- hyp.params$S[,1] + prior$alpha
        
        #update of hyperparams of beta
        hyp.params$beta <- hyp.params$S[,2:N] %*% t(hyp.params$S[,1:(N-1)]) + prior$beta
        
        #update of hyperparams of theta
        hyp.params$theta <- list(
            a = c(hyp.params$S %*% X) + prior$theta$a,
            b = rowSums(hyp.params$S) + prior$theta$b)
    }
    return(hyp.params)
}

#test###########################################################################
test <- function(){
    #'generate sample from hidden markov model with 2 - states = {1, 2}
    test.N <- 100L
    test.K <- 2L
    #'transition matrix
    test.A <- array(c(0.2, 0.8, 0.4, 0.6), dim = c(2, 2))
    
    nextState <- function(current, A) {
        sample(c(1, 2), size = 1, replace = TRUE, prob = A[,current])
    }
    test.state <- purrr::reduce(1:test.N, ~ {append(.x, nextState(tail(.x, 1), test.A))})
    
    #test data
    test.theta <- c(2, 10)
    test.X <- purrr::map_dbl(test.state, ~ {rpois(1, test.theta[.x])})
    
    #目視チェック
    plot <- ggplot(data.frame(x = c(1:test.N), val = test.X), aes(x = x, y = val)) + 
        geom_point() + 
        labs(title='plot of test data')
    print(plot)
    
    #事前分布のパラメーターを与える
    test.prior <- list(alpha = c(1, 1), 
                       beta = array(rep(0.5, 4), dim = c(2, 2)), 
                       theta = list(a = 1, b = 1))
    test.init <- list(alpha = c(0.1, 1), 
                      beta = array(rep(0.5, 4), dim = c(2, 2)), 
                      theta = list(a = 1, b = 1))
    
    estim <- CFVI(test.N, test.K, test.X, test.prior, test.init, 20L)
    
    #'状態を推定し可視化
    estim.state <- ifelse(estim$S[1,] < 0.5, 1, 2)
    estim.plot <- ggplot(data.frame(x = c(1:test.N), val = test.X, state = as.factor(estim.state)), 
                         aes(x = x, y = val, colour = state)) + 
        geom_point() + 
        labs(title='estimation of states')
    print(estim.plot)
    
    #'正解率([推定された状態 = 真の状態]であるデータの個数 / 全データ数)
    print(glue('正解率 : {sum(ifelse(estim.state == test.state, 1, 0)) / test.N}'))
}

test()

推定結果による判別の正解率({\sf accuracy \ rate})を

{\sf accuracy \ rate} = \frac{{\sf 真の状態 = 推定状態となったデータの個数}}{{\sf テストデータの個数}}

で定めて評価を行ったところ、反復回数10程度で約0.90ほどに達するようです。

f:id:kimigayoseishou:20190923220620p:plain
推定結果による2-状態判別

3.参考

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

  • 作者:須山 敦志
  • 出版社/メーカー: 講談社
  • 発売日: 2017/10/21
  • メディア: 単行本（ソフトカバー）

Rで複数のcsvファイルからデータをまとめて読み込む

小ネタ R

このところ暑い日が続いてますが，本気で20度くらい気温下がってほしいですよね．

普段の業務で表形式のデータの集計なんかしている方でしたら， csvを複数一気に読み込まないといけないなんていう状況があると思います．

今日はそんな場面で使えるRコードをご紹介します．

#'@param filePaths 読み込むCSVファイルのパスのリスト
#'@param reader csvファイルからデータを読む関数
mergeCSV <- function(filePaths, reader) {
    csv_list <- lapply(filePaths, reader)
    return(Reduce(function(x, y) rbind(x, y), csv_list))
} 

#csvフォルダ下のfrutes_1.csv, frutes_2.csv, frutes_3.csvを読み込む
filePaths <- paste(rep('./csv/fruits_', 3), c(1:3), rep('.csv', 3), sep='')

data <- mergeCSV(filePaths,read.csv)

View(data)

一応コメントしておきますと，csvから読み込んだ複数のデータをReduce関数を利用して 1つにまとめればよいという感じですね．他のReduceの使用例は例えば次のようなものです．

#Reduceの例, 1~100までの整数の和
list <- 1:100
#Reduce(function(current, result), list)で，listの先頭から2つずつfunctionにより評価を行う．
list.sum <- Reduce(function(current, result) current + result, list)
#5050