ほくそ笑む

R言語と統計解析について

【速報】神ハドリーの名著 "Advanced R" の翻訳が出ます

ビッグニュースです。

統計言語 R の業界において神とあがめられているハドリー・ウィッカム
彼は ggplot2 や dplyr などの数々の超便利パッケージを作成したことで知られていますが、その R に対する深い知識をもとに書かれた "Advanced R" という名著を昨年秋に出版しました。

Advanced R (Chapman & Hall/CRC The R Series)

Advanced R (Chapman & Hall/CRC The R Series)

本書は、これまで複雑でわかりにくいとされてきた R 言語を、エッセンスにしぼり、体系的にここまで綺麗にまとめることができたのかと、思わずうなるほどの出来で、アメリAmazon で星 4.8 という超高評価をつけられている名著です。

この本の日本語訳が、ついに出版されます!

R言語徹底解説

R言語徹底解説

訳者はこれまで数々の R 関連書を翻訳してきた石田基広先生と、R 界隈で知らない人はいないお三方。
翻訳の品質については折り紙付きです。
値段はややお高いですが、はっきりいって、この本一冊あればもう他の R 解説本はいらないんじゃないかと思えるほどの内容なので、むしろお得です。

このニュースに、R 界隈の著名人たちが一斉に反応しています。

これを読まずして R を語ることは許されません。
今後 10 年、この本が R ユーザーのバイブルとなることはほぼ間違いないでしょう。

というわけで、みなさん、今すぐ予約しましょう! ぼーっとしてたらすぐに売り切れてしまいますよ!

R言語徹底解説

R言語徹底解説

可視化で理解する「負の二項分布」

みどりぼんでカウントデータの過分散対策のために使われると書かれている負の二項分布ですが、Wikipediaの説明を読んでもよく分かりません。

そこでおススメなのが、このスライドです。

ようするに、負の二項分布は、\lambda がガンマ分布に従うようなポアソン分布だと思えばだいたい OK みたいです。

今日はこれを可視化してみます。

負の二項分布(Negative Binomial Distribution)

負の二項分布はパラメータを 2つ持ちます。成功回数を表す r と成功確率を表す p です。

統計言語 R には負の二項分布に従う乱数を生成する関数 rnbinom() があり、これらのパラメータはそれぞれ引数 size と prob に対応しています。

したがって、r=4, p=0.3 の負の二項分布は次のようにして描画できます。

negative_binom_values <- rnbinom(10000, size = 4, prob = 0.3)
barplot(table(negative_binom_values)[1:21])

これに対応するのは、\lambda がガンマ分布  {\rm Gamma}(r, p/(1-p)) に従うポアソン分布です。
このガンマ分布から \lambda を生成し、各 \lambda に従うポアソン分布から乱数を生成という手順でやってみましょう。

gamma_values <- rgamma(10000, 4, 0.3/(1-0.3))
poisson_values <- sapply(gamma_values, function(lambda) rpois(1, lambda))
barplot(table(poisson_values)[1:21], main="λがガンマ分布に従うポアソン分布")

見事に同じような分布の形状になりました。

少しわかりにくいので、この様子をアニメーションにしてみましょう。*1

最下段がガンマ分布であり、このガンマ分布からサンプリングされた値を赤い線で示しています。
この値をパラメータ  \lambda として持つポアソン分布が中段のグラフです。
このポアソン分布からサンプリングされた値を赤い線で示し、この値からなるヒストグラムを最上段に描いています。
毎回異なるポアソン分布から値がサンプリングされており、それらの値が負の二項分布を形成していく様子がわかります。

混合ポアソン分布

同じ状況を別の角度から見てみましょう。
ガンマ分布からサンプリングされた値をパラメータ \lambda として持つポアソン分布を複数作成します。

library(ggplot2)

df <- Reduce(rbind, Map(function(lambda) data.frame(lambda, x=0:20, y=dpois(0:20, lambda)), gamma_values[1:16]))
ggplot(df, aes(x=x, y=y, fill=lambda)) + geom_bar(stat="identity") + facet_wrap(~ lambda, nrow=4)

これらの分布を合成したものが負の二項分布です。

ggplot(df, aes(x=x, y=y, fill=lambda)) + geom_bar(stat="identity") + theme(legend.position = "none")

上図は 16個のポアソン分布を合成したものですが、1000個ほど合成すればほとんど負の二項分布と一致することがわかります。

library(gridExtra)

df2 <- Reduce(rbind, Map(function(lambda) data.frame(lambda, x=0:20, y=dpois(0:20, lambda)), gamma_values[1:1000]))
g1 <- ggplot(df2, aes(x=x, y=y, fill=lambda)) + geom_bar(stat="identity") + theme(legend.position = "none") + xlab("混合ポアソン分布(n=1000)") + ylab("")
g2 <- qplot(x=0:20, y=dnbinom(0:20, 4, 0.3)) + geom_bar(stat="identity") + xlab("負の二項分布") + ylab("")
grid.arrange(g1, g2, ncol=1)

この様子も無駄にアニメーションにしてみましょう。


まとめ

みどりぼんでは、個体差を持つカウントデータに対処するために GLMM(一般化線形混合モデル)を導入しています。

この本では、GLMM の例として、R の glmmML パッケージを使って「パラメータが正規分布に従う二項分布」の推定を行っています。
これに対して、負の二項分布は「パラメータがガンマ分布に従うポアソン分布」だと考えれば同じ枠組みでとらえることができると思います。

R で負の二項分布を用いた GLMM は、MASS パッケージの glm.nb() 関数によってできます。

負の二項分布がパラメータによってどのような形状をとるかの感覚を知りたい場合は、次のサイトが便利です。

みなさまの理解の一助になれば幸いです。

追記

ホクソエムの切り込み隊長 teramonagi さんが glm.nb() の使い方について詳しく書いてます。

*1:アニメーションの作成コード:https://gist.github.com/hoxo-m/3060e2fe269ff1fdf4a3

Stan でパラメータに大小関係の制約をつける

EC サイトなどでは、1 ページの読み込み時間の平均値が短いユーザほどコンバージョン率が高くなるという話があります。
この原因については、サクサク読み込みが終わる速いサイトの方がユーザは気前よく購入する、というわけではなく、単に速い回線を持つほどお金持ちなので購入しやすいという身も蓋もない話もあります。
今回はこの原因についてではなく、制約条件が付いた時のパラメータ推定の話をします。
とりあえず擬似データを作成しましょう。

true_cvr <- c(0.209, 0.126, 0.096, 0.093, 0.086, 0.077, 0.067, 0.057)

load_time <- c("0-1", "1-3", "3-7", "7-13", "13-21", "21-35", "35-60", "60+")
session <- c(1000, 6000, 4000, 1500, 700, 500, 200, 150)
set.seed(71)
cv <- unlist(Map(function(n, p) rbinom(1, n, p), session, true_cvr))

data.frame(load_time, cv, session)
  load_time  cv session
1       0-1 208    1000
2       1-3 769    6000
3       3-7 366    4000
4      7-13 142    1500
5     13-21  54     700
6     21-35  35     500
7     35-60  17     200
8       60+  11     150

ロード時間ごとに真のコンバージョン率(true_cvr)を設定して、コンバージョン数(cv)を二項分布からの乱数(rbinom)として生成しています。

このときのコンバージョン率を Stan で推定してみましょう。

// load_time.stan
data {
  int<lower=0> N;
  int<lower=0> cv[N];
  int<lower=0> session[N];
}
parameters {
  real<lower=0, upper=1> cvr[N];
}
model {
  cv ~ binomial(session, cvr);
}
library(rstan)
model <- stan_model("load_time.stan")
standata <- list(N = length(cv), cv = cv, session = session)
result <- sampling(model, data = standata, chains = 3)

この結果を見やすく可視化してみます。

library(tidyr)
df <- rstan::extract(result, pars="cvr")$cvr %>%
  data.frame %>%
  setNames(load_time) %>%
  gather(load_time, cvr)

cvr_hat <- colMeans(rstan::extract(result)$cvr)
point_df <- data.frame(load_time, true_cvr, cvr_hat) %>%
  gather(load_time, cvr) %>%
  setNames(c("load_time", "cvr_type", "cvr"))

library(ggplot2)
ggplot(df, aes(x = load_time, y = cvr)) + 
  geom_violin() + ylim(0, NA) +
  geom_point(data = point_df, aes(x=load_time, y=cvr, color=cvr_type))

パラメータの推定値(cvr_hat)と真値(true_cvr)には開きがあります。

制約条件の追加

ここからが本題です。
今、「読み込み時間が短いユーザほどコンバージョン率が高くなる」という条件が分かっています。
この条件を考慮してパラメータを推定することはできないでしょうか?
Stan では、パラメータを ordered 型で定義することによって、パラメータに順序関係の制約を加えることができます。
このような制約条件を付けた場合の Stan コードは次のようになります。

// load_time_ordered.stan
data {
  int<lower=0> N;
  int<lower=0> cv[N];
  int<lower=0> session[N];
}
parameters {
  ordered[N] cvr_rev;
}
transformed parameters {
  real<lower=0, upper=1> cvr[N];
  for(i in 1:N) {
    cvr[i] <- inv_logit(cvr_rev[N - i + 1]);
  }
}
model {
  cv ~ binomial(session, cvr);
}
library(rstan)
model <- stan_model("load_time_orderd.stan")
standata <- list(N = length(cv), cv = cv, session = session)
result <- sampling(model, data = standata, chains = 3)

非常にきれいにパラメータを推定することができました。

注意点

まず、parameters ブロックにおいて、ordered 型のパラメータを作成しています。

  ordered[N] cvr_rev;

この際、ordered 型は、「小さい順」という制約であり、求めたいパラメータは大きい順なので、これを逆順にする必要があります。
この逆順にするという操作を行っているのが transformed parameters ブロックです。
ここで、新しいパラメータ cvr に cvr_rev を逆からひとつひとつ代入する、つまり、

  real<lower=0, upper=1> cvr[N];
  for(i in 1:N) {
    cvr[i] <- cvr_rev[N - i + 1];
  }

のようにすれば良いように思いますが、それではうまくいきません。
ordered 型には lower や upper のような最小値・最大値の制約がつけられないため、−∞〜+∞までの値を取ります。
しかし、cvr には 0〜1 という制約が付いているで、cvr_rev を適当にサンプリングすると、ほとんどがはみ出してしまいます。
このような場合に便利なのが inv_logit() です。
これは、logit 関数の逆関数(logistic 関数!)で、−∞〜+∞ の値を 0〜1 にマッピングしてくれます。
したがって、

  real<lower=0, upper=1> cvr[N];
  for(i in 1:N) {
    cvr[i] <- inv_logit(cvr_rev[N - i + 1]);
  }

と書けば、cvr_rev がどのような値を取ろうが、cvr は 0〜1 の値となり、cvr の制約条件に違反しなくなります。

まとめ

Stan でパラメータに大小関係の制約を付ける方法を紹介しました。
このように、パラメータ間の制約を柔軟に加えることができるのも、ベイズ統計モデリング(Stan)の魅力かなと思います。
Enjoy!

協調フィルタリングについて発表しました

某所で協調フィルタリングについて発表してきました。資料を公開します。
機械学習とかの知識を全く持たないエンジニアさん向けの導入資料です。


協調フィルタリングは、レコメンド(推薦)に使われるアルゴリズムの一つです。

レコメンドは適用できる場面が多く、効果も分かりやすいので、実務では広く使われています。

レコメンドシステムの概要については、神嶌先生の論文が非常に分かりやすく、参考にさせていただきました。

また、協調フィルタリングについては、下記 2 つの記事を参考にさせていただきました。

みなさまのご理解の一助となれば幸いです。

Shiny アプリを収集するサイトを作った

RStudio の作成した Shiny パッケージは、R 言語で Web アプリを簡単に作成するためのフレームワークを提供します。

また、RStudio は shinyapps.io というホスティングサービスを提供しており、Shiny で作成した Web アプリを簡単にアップロードして動かすことができます。

現在、shinyapps.io には様々な Shiny アプリがアップロードされています。
例えば、次のアプリは、入力中の文章の次の単語を予測するアプリです。

テキストエリアに "to be or not to be that is the" まで入れてみると、"question" と出てきます。

また、id:ksmzn さんが作った、次のアプリも Shiny 製です。

いろいろな確率分布の形状を、パラメータを変更しながら確認できるため、非常に便利です。


このように、shinyapps.io には、日々様々なアプリが作成されてはアップロードされています。
しかし、shinyapps.io のサイトでは、自分のアプリを管理することはできても、他の人がどんなアプリを作ったのかを知ることはできません。
せっかく作られたアプリケーションも、誰にも知られずにウェブの藻屑となって消えてしまっている状況です。
また、気に入ったアプリがあっても、あとで探すときに記憶があいまいで苦労するということもあります。

  • 他の人がどんなアプリを作っているか知ることはできないでしょうか?
  • shinyapps.io に新たに追加されたアプリを知る方法はないでしょうか?
  • 気に入った Shiny アプリを記録してまとめておくことはできないでしょうか?


というわけで、作りました。

shinyapps.io に登録されたアプリを収集するサイトです。

原理的には、Google 検索にひっかかった Shiny アプリを片っ端から登録してます。


トップページはアプリが見つかった順番で並んでおり、新しいものほど上に表示されます。
画像をクリックすると、アプリの詳細ページに行き、そこのリンクからアプリにアクセスすることができます。

また、右上のログインボタンから Google アカウントでログインすると、アプリをお気に入りに入れることができるようになります。
お気に入りに入れるには、詳細ページの星マークを押します。

お気に入りに入れたアプリはマイページ(ログインしたら右上にボタンが出ます)で一覧表示できます。


これにより、

  • shinyapps.io に新しく登録されたアプリを知ることができる
  • 作者名をクリックすれば、その人が作ったアプリの一覧が見れる
  • 気に入ったアプリを登録して自分だけの一覧を作れる

ということができるようになりました。
やったね!


また、このサイトに登録された新規 Shiny アプリの情報を、Twitter で流しています。

ぜひフォローして最新情報を追ってください!


というわけで、Shiny アプリを収集するサイト、Shinyapps Gallery を作りました。ぜひ使ってみてください。

それではよき Shiny ライフを!

データ分析のプロを目指すエンジニア必読の書

福島真太朗『データ分析プロセス』を読みました。

データ分析プロセス (シリーズ Useful R 2)

データ分析プロセス (シリーズ Useful R 2)

「データ分析」とひとことで言っても、色々な人たちが色々な考え方で「データ分析」をやっていると思います。

その一大勢力として挙げられるのが「Excelで集計だけしてる人たち」です。これは特にマーケター出身の人が多いのではないでしょうか。*1

一方、最近のデータ分析界隈に増加していると思われるのが「機械学習ツールにデータを入れるだけ」の人たちです。
特にエンジニア出身の方が、上司に「データ分析が流行ってるみたいだから今日からデータ分析やって」と言われて泣きながらやっているケースが多いのではないかと。

そういう人たちは、機械学習についてある程度理解しており、機械学習ツールにデータを入れて結果を出すことはできるのだけれど、分類精度が悪かったり、うまくレコメンド結果が出ないなどで悩む人も多いかと思います。

そしてふと、「そもそもデータ分析ってこれでいいの?」ということに突き当たります。
「そもそもデータってそのまま機械学習に入れていいの?」
「精度が悪いときにはどうすればいいの?」
「そもそもどうやって精度を測ればいいの?」
「みんなどうやってるの?」

つまり、「データ分析のプロセスってどういうものが一般的なの?」ということが彼らの知りたいことです。

これに対する答えが、福島真太朗『データ分析プロセス』に載っています。

データ分析プロセス (シリーズ Useful R 2)

データ分析プロセス (シリーズ Useful R 2)

この本は、データ分析の実務家である著者が、データの集計・前処理から、予測モデルの構築、データマイニングまで、データ分析のプロセスを一通り行うための知識をまとめたものです。

例えば、「4.1 予測モデルの構築」では、

  1. 予測問題の設定
  2. 特徴量の構築
  3. ハイパーパラメータの最適化
  4. 予測モデル構築

のように、予測モデル構築までのフローを提示し、各プロセスで何をやるのかについてコード付きで解説してくれます。


特に私が本書を読んで感嘆したのは、データハンドリングと前処理だけで、本書の半分を費やしているということです。
機械学習ツールにデータを入れる際、そのままのデータを入れてもうまく精度は出ません。
「Garbage in, garbage out(がらくたを入れてもがらくたしか出てこない)」という言葉があるように*2、どんなに高性能な機械学習アルゴリズムを使っても、データがごみならごみのような結果しか出すことはできません。
そこで重要なのが、欠損値の処理、外れ値の検出、変数の離散化などの前処理です。

  • 欠損のある観測値を単純に削除していませんか?
  • 外れ値を考えずにデータをそのまま入力していませんか?
  • 連続変数をそのまま使っていませんか?

本書では、これらのやり方について、様々な手法の紹介をしています。
著者は恐ろしいほどの量の文献を調査し、伝統的な手法から現代的な手法まで、バランスよく厳選された手法の紹介を行っています。


また、本書の後半では、caret を使った予測モデルの構築について書かれています。
caret は、R 言語から使える機械学習フレームワークで、クロスバリデーションを用いたパラメータチューニングや属性選択(feature selection)など、モデル構築を行うための便利な機能をたくさん持っています(並列計算もできます)。
まだ caret を使ったことがないという方は、かなり損してると思います。


一方、この本では機械学習手法そのものについては、あまり説明していません。
取り扱っている機械学習手法は、ランダムフォレストとサポートベクタマシンくらいです。
これは、この本が主眼を置いているのは、

  • 欠損処理、外れ値処理、離散化などの前処理を行い、機械学習に入力するデータを作ること
  • パラメータチューニングや属性選択手法を行い、精度の良いモデルを作ること

であり、これらは「どんな機械学習手法を使うとしても、やり方は同じ」だからと考えられます。
つまり、この本で学べる知識は、今は存在していない新しい機械学習アルゴリズムが現れたとしても、全く同じように適用可能です。


というわけで、この本をおススメしたいのは「コードが読めるエンジニア」であり、「欠損や外れ値の処理を行わず、そのまま機械学習に入れちゃってる人」です。
そういう人たちが「実務で通用する」データ分析を学ぶためには、必読の書と言えるでしょう。
また、実務で機械学習を使ったことがない分析者志望の人にとっても、実際の現場におけるデータ分析の流れを、具体例を通して知ることのできる貴重な本だと思います。

『データ分析プロセス』、ぜひ手に取ってみてください。

関連

本書の各章でどのような手法が紹介されているかは、下記の書評が参考になります。

*1:そういう人たちにおススメなのはこちらです:http://d.hatena.ne.jp/hoxo_m/20141010/p1

*2:Wikipediaにも載ってます: https://en.wikipedia.org/wiki/Garbage_in,_garbage_out

Spark Meetup 2015 で SparkR について発表しました #sparkjp

「初めてのSpark」刊行記念 Spark Meetup 2015で発表してきました。
タイトルは「はじめての SparkR」です。

SparkR は R 言語から Spark を使うためのパッケージで、Spark 1.4.0 から正式にサポートされました。
「RStudio から使える」「dplyr ライクなデータ操作」など、R ユーザーには嬉しい設計になっていますが、機能としてまだ不十分なところが多いです。

発表では、最新の Spark 1.4.1 の SparkR を元にお話しましたが、Meetup 当日に Spark 1.5.0 がリリースされたため、一瞬にして内容が古くなるというミラクルが起きました。

また、発表中に「Parquet とか Avro とかの読み方がわからない」と言ったら、親切にも shiumachi さんが教えて下さいました。

本 Meetup 全体の内容については、下記まとめレポートが詳しいです。

Spark やってみたいという方は、『初めてのSpark』買いましょう!

初めてのSpark

初めてのSpark