ベイジアン線形回帰モデルの式変形とRでのギブスサンプリングの適用

今回は特に目新しい手法というわけでもなく、線形回帰モデルのギブスサンプリングについて忘備録として残しておきたいと思います。
ベイジアン線形回帰モデルはプログラミング言語で言う、Hello World!的なものなので、あえてブログで取り上げる必要があると考えていないのですが、導出をしては忘れの繰り返しが嫌なので自分のために残しておこうと考えました。加えて、Stanのありがたみを感じられ、Stanへのコミットメントが増すのではないかとも期待しています。

・モデル
・数式の展開
・Rのコードの紹介
・おわりに
・参考情報

モデル

東北大学の照井教授の『ベイズモデリングによるマーケティング分析』に載せられている表記に従い、以下のように記します。

説明変数の数がk個の正規線形モデル

を考える。その場合、尤度関数は

となる。

係数パラメータの事前分布や条件付きの誤差分散の事前分布は以下のように設定する。(βは正規分布に従い、σ2|βは逆ガンマ分布に従う。)

数式の展開

私が大学院生だった時に、数式の展開をどう進めるかを手っ取り早く知る方法としては、「ネットに上がっている海外の大学院の講義資料を漁る」という作戦を取っていました。こうすることで数学のセンスがそれほど高くなくても、理解し進めることができました。今回に関してもおそらく、わかりやすく解説している海外の研究者がいるはずだと思い漁ってみたところ、コロンビア大学の機械学習の講義資料を見つけることができました。

資料はこちらのPDF(Course Notes for Bayesian Models for Machine Learning)で126ページにもなっていますが、導出のステップなどが非常に丁寧に書かれています。

それでは、今回の講義ノートを参考にしながら、線形モデルにおいて、ギブスサンプリングを行うところまでの導出を行いたいと思います。

まず、同時事後分布を以下の左辺のように置き、ベイズの定理を用いて右辺のように表記する。

次に、条件付き確率の定義と先程の尤度関数から以下のようになる。

yが与えられたもとでのp(y)は一定のため、比例している分子だけを残すと以下のようになる。

同時事後分布に事前分布の関数を代入していくと、

となる。両辺について対数を取ると、

となる。ここでβやσ2についての事前分布の形状から、同時事後分布におけるβやσ2について整理するための目標となる形状を確かめる。
まず、βはp(β)の定義より、対数を取りβについて整理すると、

となる。つまり、1/B0や1/B0・β0に該当する表現を先程の対数を取った同時事後分布から得ることを目標とする。
他方、σ2についても同様に、p(β|σ2)の定義より対数を取りσ2について整理すると、

となる。つまり、ν0やδ0に該当する表現を、同じく対数を取った同時事後分布から得ることを目標とする。

以上のパラメータごとの目標とする形状になるように各々のパラメータについて、対数を取った同時事後分布を整理する。

まずはβについてまとめ、関係のない項をconst.にする。

先程もとめた目標の形状を当てはめると以下のようになる。

よって、βの事後分布は以下のようになる。

他方、σ2についても同様に、関係のない項をconst.にし、目標の形状にまとめると以下のようになる。

目標の形状と比較すると以下のようになる。

よって、σ2の事後分布は以下のようになる。

Rのコードの紹介

条件付き事後分布からβやσ2の従う分布の形状がわかったので、それらを使ってRでギブスサンプリングを行います。先日、たまたま見つけた線形回帰モデルのギブスサンプリングのRのソースコードを拝借しようと思います。

ギブスサンプリングでは、先程導出した条件付き分布からβ→σ2と交互にサンプリングしていきます。それを記述したRコードは以下の通りです。

先程導出したβの事後分布である正規分布からのサンプリングの後(15~17行目)、そのサンプリングしたβを用いて、同じく先程導出したσ2の事後分布である逆ガンマ分布からサンプリングし(19~21行目)、それを指定した回数だけ繰り返し、所定の数まではバーンインとして除外します。(25~27行目)こうして導出した数式と、ギブスサンプリングのコードを見比べると理解が捗ると思いました。

実際に、先程のGitHubのソースコードを回してみると、以下のようにギブスサンプリングのイタレーションのプロットや、パラメータの事後分布を確認できます。

おわりに

シンプルなモデルですらこれだけ導出に手間がかかるということからも、Stanなどの確率的プログラミング言語のありがたみは非常に大きいなと思いました。こうして残すことで今後忘れたとしてもすぐに思い出せる気がします。
しかしながら、Stanでは事前分布と尤度を指定してしまえば、事後分布を計算し、知りたいパラメータについて解いた条件付き分布からサンプリングしてくれるわけですから、研究者の寿命を伸ばしたと言っても過言ではないかもしれません。

参考情報

[1]John Paisley (2016), “Course Notes for Bayesian Models for Machine Learning”, Columbia University
[2]照井伸彦 (2008), 『ベイズモデリングによるマーケティング分析』, 東京電機大学出版局
[3]須山敦志 (2017), 『機械学習スタートアップシリーズ ベイズ推論による機械学習入門』, 講談社
[4]stablemarkets,BayesianTutorials/MultipleLinearReg/multiplelinearreg.r

[Stan]生存時間分析のコードと便利なデータセットについて

はじめに

仕事で生存時間分析を使うことは結構あるのですが、マーケティングの良いデータセットがない印象でブログにしにくいと感じていました。また、Stanでの生存時間分析の事例もあまり把握していません。そこで使えそうなデータセットやStanのコードを探して、そのデータに対して生存時間分析を適用してみたいと思います。

目次
・生存時間分析とは
・生存時間分析で使えるデータ
・生存時間分析をマーケティングで使う際の用途
・先行研究
・生存時間分析で使えるデータセット
・Stanでの実行例
・おわりに
・参考文献

生存時間分析とは

生存時間分析は、ある時点から興味のあるイベント(マーケティングだと解約など)が発生するまでの期間を分析対象としています。データを手に入れた時点で、すでに解約して、真の累積の契約期間が判明している人と、解約しておらず今後いつ解約するかわからない中での累積の契約期間が残されている人のようなデータを扱うことが多いです。ここでの後者をcensoring(打ち切り)されたデータと呼びます。

生存時間分析をマーケティングで使う際の用途

ブログなどを読み漁る限り、以下の用途で生存時間分析を活用できるようです。

  • 顧客のサービス離脱率の予測、離脱原因の特定
  • 顧客がマーケティングキャンペーンに反応するまでの期間の長さ
  • 故障率の予測、故障原因の特定

先行研究

Stanを用いた分析事例は、調べた限りですが以下のモデルがあるようです。

  • 指数分布のモデル
  • Weibull(ワイブル)分布による比例ハザードモデル
  • ハザードの対数値についてのランダムウォークモデル
  • 2階差分のマルコフ場モデル(生存時間の確率分布は正規分布)
  • 1階差分のランダムウォークモデル(生存時間の確率分布は正規分布)

生存時間分析で使えるデータセット

事例を調べる過程で見つけた、生存時間分析に適したデータセットは以下の通りです。

どうやら、マーケティング、HR、離婚、再犯と幅広いデータがオープンソースで手に入るようです。

Stanでの実行例

今回は、「Using Customer Behavior Data to Improve Customer Retention」のデータセットを用いて、先行研究のソースコードにより生存時間分析をしてみようと思います。データは電話会社の顧客の解約に関するもので、様々な顧客の履歴データなどが用意されています。
先行研究のソースコードはWeibull分布を想定した比例ハザードモデルです。今回は決済の電子化の有無と離脱の関係を確かめてみます。なお、今回の打ち切りデータは契約期間となります。

まずはStanのコードはこちらです。Xobs_bgに説明変数が来るようにデータを用意しておきます。

続いて、このStanコードをキックするためのRのソースコードです。 元のデータが7043件と多すぎるのでランダムサンプリングしています。サンプリング数を8000、チェイン数を4にして実行します。(なお、可視化のソースコードもあるので結構長くなっていますので。頑張ってスクロールしてください。)

Rhatは全て1.05以下になっています。

traceplotを見る限り、重なり合っているので問題なさそうです。

各パラメータごとの自己相関係数です。こちらも問題なさそうです。

推定したパラメータの分布です。

横軸は推定した継続期間です。決済の電子化をしていない消費者は、契約期間の短い際の確率密度が低い傾向があるようです。

どうやら離脱率に関して決済の電子化をしていない消費者は、そうでない消費者よりも低いようです。

こちらは95%で取りうる範囲をそれぞれプロットしたものです。

おわりに

Stanで生存時間分析を行うという事例はそんなに多くはないものの、業界の長たちが良いコードを作成してくれていました。また、面白そうなデータセットも見つけることができました。このようなデータがもっと広まっていけば、マーケティングにおける生存時間分析がより活発に行われるのかもしれません。

参考文献

[1] 豊田秀樹 (2017) 『実践 ベイズモデリング -解析技法と認知モデル-』 朝倉書店
[2]生存時間解析入門
[3]比例ハザードモデルはとってもtricky!
[4]Stanで生存時間解析(Weibull 回帰)
[5]生存時間分析をStanで実行してみた
[6]階層ベイズ生存解析を用いたwebサイトの訪問者分析に関するStanでの実装
[7]生存時間分析 – ハザード関数に時間相関の制約を入れる
[8]Bayesian Survival Analysis 1: Weibull Model with Stan
[9]Bayesian Inference With Stan ~062~
[10]生存時間解析について – 概要編
[11]Survival Analysis for Employee Attrition ※kaggleで提供されているHR系のデータをサバイバル分析に用いている。
[12]Survival Analysis with R※Random Forests Modelによる生存時間の推定がなされている。
[13]Survival Analysis with R and Aster ※服役後の犯罪に関する分析や、離婚の分析などをしている。
[14]Survival Analysis of Mobile Prepaid Customers Using the Weibull Distribution(ダウンロード注意)

[Stan]項目反応理論(IRT)の段階反応モデルでbaysemのアンケートデータの分析をしてみる

はじめに

stanのユーザーガイドを見ていて、項目反応理論(IRT)についての章があり気になりました。勉強会のLTなどで手法の名前をちらっと聞いたことはあったのですが、使い道について調べていませんでした。ビジネスにおける実活用もしやすそうだと思ったので、カジュアルに分析して備忘録として残したいと思います。

目次
・項目反応理論(Item Response Theory:IRT)とは
・ビジネスでの適用可能性について
・データ
・モデルの推定
・結果の解釈
・おわりに

項目反応理論(Item Response Theory:IRT)とは

関西学院大学の教授のブログによると、

項目反応理論とは、テストについての計量モデルで、問題に対する正解・不正解のデータから、問題の特性や、回答者の学力を推定するためのモデルです。

とあります。また、Wikipediaによると、TOEFLの問題の評価のために使われているそうです。

主に、バイナリーと順序変数のモデルがあるようで、以下の母数がモデルに想定されています。どちらもほぼ同じです。

回答が2値変数のモデル

  • 2母数のロジスティックモデル
    • 特性値(例えば、広告配信の満足度とか)
    • 識別度母数(項目特性曲線の傾き)
    • 困難度母数(項目特性曲線の切片)
    • 定数

回答が順序変数のモデル(まずい < まぁまぁ < おいしい)

  • 段階反応モデル
    • 特性値(例えば、広告配信の満足度とか)
    • 識別度母数(項目特性曲線の傾き)
    • 困難度母数(項目特性曲線の切片)

※項目特性曲線は横軸に特性値、縦軸に質問の正答率を取ったものです。

ビジネスでの適用可能性について

  • 顧客のアンケート結果の解釈
    • 異質な集団間の得点を比較可能
    • 異なる尺度間の得点を比較可能(昔のアンケートだと5段階、今のアンケートは7段階などの状況はビジネスデータでありうる。)
  • 人事評価のバイアスの統制
    • 採用面接時の個人特性の正当な評価
  • アンケート項目の項目削減によるアンケートコストの低減
    • 各アンケート項目が理解されたかどうかを分析し、一つ一つのアンケート項目の精度を高める

データ

今回扱うデータはbaysemパッケージに入っているデータセットです。Yellow Pagesの広告プロダクトにおける満足度サーベイの回答データで、全ての回答は1から10のスケールで点数が付けられています(1がPoorで10がExcellent)。質問数は10個で、回答数は1811件です。

各質問の内容(baysemパッケージのドキュメントに載っていました。)
q1:全体の満足度

価格について
q2:競争的な価格設定
q3:昨年と同じ広告の最小値に対しての価格の引き上げ
q4:消費者の数に対しての適切な価格設定

効果について
q5:広告の購入の潜在的な影響
q6:広告を通じて自身のビジネスへの集客ができたか
q7:多くの消費者にリーチしたかどうか
q8:年間を通じて消費者に対する長期での露出があったか
q9:多くの家計やビジネスを必要としている人に届いたかどうか
q10:ビジネスを必要としている地理上のエリアに届いたかどうか

今回のIRT適用における特性値は、「広告プロダクトに関する満足度の傾向」としてみたいと思います。

モデルの推定

今回は教科書にならって以下の段階反応モデルを用います。

ここでaは識別力(広告の満足度が高まりやすいかどうか)、bは境界パラメータ(回答カテゴリ間の境界値)、θは特性(回答者がどれだけ広告に満足しているか)を表しています。Dは定数項で、以下では1とおいています。cはアンケートの回答のカテゴリ番号です。今回の例では10段階の評価が入ることになります。最後に、uは反応を、jは質問の番号を表しています。

実践 ベイズモデリング -解析技法と認知モデル-

こちらの本のサポートサイトからダウンロードできるzipファイルにstanのコードやRコードがありますので、そちらを利用しています。

モデルですが、以下のような設定となっています。

こちらをキックするためのRコードです。

処理時間としては、2014年末モデルのMacbook Proのcorei5、メモリ8GBで数時間程度でした。(正確な時間はわかりませんが、寝て起きたら計算が終わっていました。)


どうやら収束してそうです。


Rhatも1.1未満におさまっています。

結果の解釈

まず、推定した特性値の値のユーザーごとの平均値を求めて、ヒストグラムを描いてみます。どうやら、上限周辺にやたらと高い評価をしてそうなユーザーがいるようです。

最後に、項目特性曲線を質問ごとに、そして回答ごとに描いてみようと思います。

質問1~10に関して、10段階の回答ごとの項目反応曲線を以下に描いています。上まで戻るのが面倒なので、質問内容を再掲します。

q1:全体の満足度
q2:競争的な価格設定
q3:昨年と同じ広告の最小値に対しての価格の引き上げ
q4:消費者の数に対しての適切な価格設定
q5:広告の購入の潜在的な影響
q6:広告を通じて自身のビジネスへの集客ができたか
q7:多くの消費者にリーチしたかどうか
q8:年間を通じて消費者に対する長期での露出があったか
q9:多くの家計やビジネスを必要としている人に届いたかどうか
q10:ビジネスを必要としている地理上のエリアに届いたかどうか

これらの傾向から、9〜10点を獲得するにはある程度は特性値が高まる必要がある質問としては、q1〜q6のように見えます。価格や購買など自身のビジネスに直結しそうな質問が多い印象です。逆にふわっとした質問であるq7~q10は特性値が低くても9〜10点を取れる可能性が高い傾向があります。

おわりに

Stanのユーザーガイドを読むことで、普段自分が業務で扱っているアプローチなどが如何に限定的であることが実感できました。今回はIRTのアンケートデータへの適用事例を知れ、そこから様々な文献や便利なコードに至ることができました。社内のアンケートデータへの適用は面白そうだと思いますので業務で使ってみようと思います。

参考情報

[1] 豊田秀樹 (2017) 『実践 ベイズモデリング -解析技法と認知モデル-』 朝倉書店
[2] Yoshitake Takebayashi (2015) 「項目反応理論による尺度運用」 SlideShare
[3] 持主弓子・今城志保 (2011) 「IRTの組織サーベイへの応用」
[4] 清水裕士 (2017) 「項目反応理論をStanで実行するときのあれこれ」 Sunny side up!
[5] 清水裕士 (2016) 「Stanで多次元項目反応理論」 Sunny side up!
[6] 小杉考司 (2013) 「項目反応理論について」
[7] Daniel C. Furr et al. (2016) “Two-Parameter Logistic Item Response Model”
[8] daiki hojo (2018) “Bayesian Sushistical Modeling” Tokyo.R#70
[9] abrahamcow (2017) 「[RStan]項目反応理論の応用でフリースタイルダンジョン登場ラッパーの強さをランキングしてみた」

Python/Rもくもく会をプライベートで開催するための参考図書・資料をまとめる

はじめに

社内で定時後に有志で勉強会というか、その場に集まってPythonやRをもくもくと勉強をするもくもく会を開きたいと考えています。目的としては分析スキルの向上や機械学習ができるようになりたいとかいう個々人の願いを叶えることです。
色々なスキルレベルのメンバーが参加することが予想されるので、皆を幸せにするためにもレベルに応じた良い教材が必要だと思いました。
ここでは、レベルに応じて適切な教材などを忘備録として残していきたいと思います。
(私自身、全てのレベルの対象者に適切な教材を網羅しているわけではないので、随時更新していこうと思います。)

受講対象について

受講対象(PythonやRをまともに触ったことがない人)は2軸で分けるとすると以下のようになると思います。

・プログラミング経験あり/経験なし
・数学の心得あり/心得なし

  • プログラミング経験なし&数学の心得あり(アルキメデス)
    理系出身の人がメインだと思います。学部・学科によっては全然扱わないですよね。数的な思考は得意だが、それを活かすスキルが不足しているような人でしょう。眼の前におかれた数学の問題を紙とペンで解くことはできるが、仕事で使えないという感じ。私も偉そうなことは言えないですが、コードが荒れがちなので周りに良い先生がいたほうが良いと思います。

  • プログラミング経験なし&数学の心得なし(葉っぱ隊)
    一番習得に時間がかかると思います。野球やったことないのに、野球選手になりたいという人に皆さんは違和感を感じるでしょう。イメージはそんな感じです。一番時間がかかるからこそ、挫折しないための教材選びが重要かもしれません。スキル的に全裸なので、葉っぱ隊と名付けましょう。

  • プログラミング経験あり&数学の心得あり(デーサイ候補)
    最も頼もしい存在です。教科書をお渡ししておけば勝手に成長すると思います。ある程度経験を積めば分析業務を任せても良いと思います。

  • プログラミング経験あり&数学の心得なし(進捗ありマン)
    各種手法の原理を知るまではそれなりに時間がかかると思いますが、手を動かして何ができるかをすぐに味わえるので、モチベーションを維持しながら学んでいきやすいと思います。コード自体は実行できるので進捗ありマンと名付けてみましょう。

この2軸でPythonとRに関する便利な資料を探したいと思います。
ただし、どの本に関してもどのレベルの人が買っても良いとは思います。ただ、数学の心得がない中で、テイラー展開とか平均値の定理とかラグランジュ未定乗数法などの表現を目にした際に、挫折してしまう可能性があるので、適した書籍から順次広げていくのが良いと思います。なお、今回はPCでもくもくと進めれそうな書籍を選んでいます。紙とペンで進める本も重要なのですが、そのようなかた向けの書籍は取り上げていません。

アルキメデス向けの教材

Python

R

  • みんなのR 第2版
    『Rによるデータサイエンス』と迷ったのですが、プログラムの実行結果がそのまま載っている印象だったので、こちらの本がプログラミング初心者には優しいと判断しました。ほとんど数式は出てこないのですが、一般化線形モデルや時系列解析などもカバーしてくれています。また、データの前処理に関する記述もこちらの本の方が手厚いです。

葉っぱ隊向けの教材

Python

  • Pythonスタートブック [増補改訂版]
    本当にプログラミングがはじめての人向けの本です。まずはプログラミング自体に慣れたほうが良いと思います。

  • プロゲートのPython入門講座
    妻におすすめされた講座です。無料枠でもある程度学びがあるようです。環境を構築しなくても良いという点が非常に葉っぱ隊に適しているとのことです。

R

  • Rによるやさしい統計学
    Rのインストールあるいは統計学の初歩のところから、応用まで幅広く説明している本です。数式はあまり出てきませんがコードが載っているので、手を動かすことができると思います。

読み物

デーサイ候補向けの教材

Python

R

進捗ありマン向けの教材

R

  • RStudioではじめるRプログラミング入門
    プログラミング経験のある進捗ありマンであれば、R言語の扱い方をまずは知りたいだろうと思います。関数の書き方やヘルプページの使い方、オブジェクトの説明、S3の話などが詳しく書かれています。

  • 新米探偵、データ分析に挑む
    R Studioのインストール方法なども載っているので、進捗ありマンなら最初から最後まで実践できると思います。数式もほとんど出てきません。色んな分析事例をRで取り組むことで分析業務のイメージも付いてくると思います。

  • RユーザのためのRStudio[実践]入門−tidyverseによるモダンな分析フローの世界−
    R言語について何となくつかめた進捗ありマンがモダンな記法であるtidyverseを効率よく学べる良い本です。データ整形・クロス集計・可視化がモダンな記法で書けるようになると結構楽しいと思います。

Python

今後について

そもそもPythonやRに触れたことがない人にとって、Tokyo.Rの初心者セッションは少し適していないのかなと思ったので、今回は取り上げていないですが、一通り使い方をわかってもらえたら初心者セッションの資料を使ったもくもく会も開きたいと思います。最終的にはKaggle部とかを作るとかになるのかもしれませんが、そこまで行けるか行けないか。

Uplift Modeling用のパッケージtools4upliftを使ってみた

はじめに

今回は、今後仕事で使いたいという思いもあり、RでUplift Modelingに関して便利なパッケージがないか探した結果、2019年に登場したばかりのtools4upliftの存在を知りました。アップリフトモデリングのモチベーションに関しても簡単に説明しながら、サンプルデータで実践してみようと思います。

・Uplift Modelingとはなにか
・Uplift Modelingの卑近な例え話
・Uplift Modelingのサンプルデータ
・tools4upliftについて
・tools4upliftでCriteoデータを試してみる
・『仕事ではじめる機械学習』の9章のコードをCriteoデータに試してみる
・おわりに
・参考文献

Uplift Modelingとはなにか

きちんとした説明は、あまりにも今更感があるので説明は端折りたいと思います。既出の文献がありますので、そちらを熟読ください。

Uplift Modelingの卑近な例え話

自分が吉野家のマーケティング担当だとしましょう。吉野家のアプリで割引クーポンを顧客にばらまくことができるとします。
マーケターとして重要なのは、割引クーポンを渡したことをきっかけとして吉野家に足を運び購入する顧客を増やせるかどうかになります。

マーケターの手元にあるのは、割引クーポンをばらまいた顧客とばらまかなかった顧客、そして吉野家で牛丼を食べたかどうかのデータです。
以前のマーケティング担当者がランダムにクーポンをばらまいていたことが重要なポイントです。

このデータから、顧客は以下の4分類に分かれます。

  • 無関心:割引クーポンをばらまこうが我関せず。そもそも吉野家に行く気はない。
  • 説得可能:普段、牛丼が安いすき屋にばかり行っているが、割高に感じている吉野家に負い目を感じている。割引クーポンで揺さぶられ来店する。
  • 天の邪鬼:吉野家コピペのように、割引クーポンを握りしめた家族連れに遭遇したくないので、割引クーポンをばらまかれたら来店しないような客。
  • 鉄板:毎日決まった時間に吉野家に行くことを心に決めている客。

マーケターは割引クーポンをばらまいた顧客と割引クーポンをばらまいていない顧客にデータを二分し、それぞれ機械学習のための訓練用データとテスト用データを用意します。

つまり、「割引クーポンをばらまいた顧客」の訓練用データとテスト用データと「割引クーポンをばらまいていない顧客」の訓練用データとテスト用データの計4つのデータセットを用意します。

まず、牛丼の購入の有無を教師とした訓練用データでロジスティック回帰モデルなどを推定します。
その結果、「割引クーポンをばらまいた顧客」から推定したモデルと、「割引クーポンをばらまいていない顧客」から推定したモデルが手元に残ります。

2つのテスト用データを1つにまとめて、先程推定したモデルを用いて、牛丼の購入確率を求めます。モデルは2つあるので、予測結果がテスト用データ1つに対して2つあることになります。

その予測結果の比(「割引クーポンをばらまいた顧客」モデルベースの予測値÷「割引クーポンをばらまいていない顧客」モデルベースの予測値)をアップリフトとみなします。

以下の図はこれまでの説明を図にしたものです。

アップリフトがどの程度の水準であれば、説得可能なユーザーが多いのかを探っていくことで、吉野家のアプリにおいて、どのユーザーに割引クーポンを発行するべきかがわかることになります。

Uplift Modelingのサンプルデータ

残念なことに吉野家のアプリのデータはありません。そこで今回は公開データを利用します。
以前より、The MineThatData E-Mail Analytics And Data Mining ChallengeのメールのデータがUplift Modelingで非常にしばしば取り上げられるデータでしたが、Twitterで他にデータないのかとぼやいたところ、2名の方にCriteo Uplift Prediction Datasetを紹介していただきました。

余談ですが、Criteo社と言えばディスプレイ広告のキング的な存在で、少し商品のリンクを踏んだだけであっという間に広告がレコメンドされますよね。自社で出稿用バナーを作っていましたが、CVRが高くなる良いクリエイティブを作ってきたのか、単にCriteo社のアルゴリズムが優秀なだけなのか非常に気になるところでしたね。

Criteo社が提供してくれている今回のデータは、2500万行に及ぶユーザーのデータで、プライバシー保護の観点から特徴量は復元できないような形式で提供されています。バイナリーのラベルとしては訪問やコンバージョンなどがあり、データ全体に占める処置群の割合は84.6%となっています。要は、吉野家で言う割引クーポンをばらまいた顧客が全体の84.6%に及ぶということです。

tools4upliftについて

2019年1月に公開されたRのUplift Modeling用のパッケージです。

  • 特徴量における連続値をカテゴリ変数にする際に、最適な階級値を求めてくれる関数
  • アップリフトモデリングの可視化する関数
  • アップリフトモデリングにおける特徴量選択ができる関数
  • アップリフトモデリングにおけるモデルのバリデーションを行う関数

などが提供されており、ちょいとRを触れるマーケターにとって、アップリフトモデリングにおける試行錯誤がかなりしやすくなる便利なパッケージだと思いました。
なお、このパッケージで扱っているモデルはロジスティック回帰になります。介入データをもとに推定したモデルの条件付き確率と非介入データをもとに推定したモデルの条件付き確率の差をアップリフトとして推定しています。

このパッケージの解説論文においては、アップリフトモデリングの評価指標としてQini曲線というものが提案されていました。Qini曲線はローレンツ曲線のようなもので、Qini曲線とランダムに割り当てた際のアップリフト量の差分の合計をQini係数と定義しています。

tools4upliftでCriteoデータを試してみる

こちらはアップリフト値の予測値の上位から右に並べた際のアップリフトの増大のグラフになります。20%あたりでピークになるようです。

こちらはアップリフト量の棒グラフです。20%の階級値を超えたらガクンと下がるのがわかります。

なお、Qini係数は0.03233551でした。

『仕事ではじめる機械学習』の9章のコードをCriteoデータに試してみる

tools4upliftの結果を鵜呑みにするのもあれなので、『仕事ではじめる機械学習』の9章のコードを使ってアップリフトモデリングを実践してみます。コードは丸パクリですが、謹んで掲載させていただきます。

こちらの図はアップリフト値の階級値ごとのCVRです。最上位のアップリフト値はCVRの差が大きいですが、上位40~50%程度のアップリフト値のときにCVRの差が最も大きいようです。

アップリフト値の順位とCVRの図です。順位が低くても処置群のほうがCVRがわずかに高いようです。

アップリフトのスコアとCVRの関係です。2未満であればCVRは処置群が上回っていますが、一様な傾向はなさそうです。

コンバージョンレートの差に対象群の人数を掛けることでliftを算出したものです。アップリフトスコアが1~2点であれば儲かるようです。

tools4upliftと出している指標が違うので比較ができないのが難点に思いました。tools4upliftはオートマチックな感じで便利なのですが、『仕事ではじめる機械学習』の9章を正義として進めたいので、どうにか揃えれるようにしていきたいと思います。

おわりに

tools4upliftというマーケターにとって銀の弾丸になりそうなパッケージの存在を知ることができ、実際に非常に便利そうな関数が用意されているのがわかりました。ただ、開発されたばかりのパッケージなのでそこまで結果を信じていません。『仕事ではじめる機械学習』本の結果と揃えたいなと思いました。その点がはっきりすれば業務で使ってみるのも良いですし、任意のマーケターに安心して共有できると思います。

参考文献

[1] 有賀康顕・中山心太・西林孝 (2018) 『仕事ではじめる機械学習』 オライリージャパン
[2] Mouloud Belbahri, Alejandro Murua, Olivier Gandouet, Vahid Partovi Nia (2019). “Uplift Regression: The R Package tools4uplift”, arXiv:1901.10867 [stat.AP]
[3] ohke (2019) 「Uplift modelingで施策が効く人を見極める」 け日記
[4] usaito (2018) 「Uplift Modelingで介入効果を最適化する」 Qiita

RのContextualパッケージをいじってみた際のメモ書き

はじめに

このブログの私の中での位置づけは、今後仕事で使いそうなものを調べて書き溜めるというところにあります。仕事で使っているものはブログに載せないというスタンスでもあるのですが、出来るだけ先回りしておきたいところです。今回は、昨年のJapan.RやTokyo.Rで紹介されていたcontextualパッケージを触ってみたというゆるふわな内容となっています。

目次

・バンディット問題とは
・マーケティング関連でバンディット問題が役に立つ場面
・バンディット問題で出てくる数学的な知識と方策
・Contextual Bandit問題とは
・Contextualパッケージでできること
・サンプル実行
・おわりに
・参考情報

バンディット問題とは

「選択肢の集合から1つの要素を選択して、その選択肢に対する報酬を得るものの、他の選択肢の報酬情報は得られないというプロセスを繰り返す設定において、報酬の合計値を最大化することを目指す逐次決定問題」とされています。バンディットは昔ながらのスロットマシンが客からお金をむしり取ること(盗賊)にちなんでいるそうです。胴元は盗賊ということなんでしょうか?

大学時代の知人は毎日パチンコ屋に行ってから講義に行っていましたが、出そうな台・出そうな店を転々としていましたが、あれはバンディット問題を彼なりに解いていたのでしょう。当時はサクラの台というのがあったらしく、3000円ほど投資すれば大当たりになるのだとか。そしてその大当たりに釣られて他の客が頑張るという意味で、サクラの台だそうです。

マーケティング関連でバンディット問題が役に立つ場面

私はマーケティング×データ分析を生業としているので、マーケティング方面にしか関心がないのですが、バンディット問題は役立つ可能性が十分にあるというか既に一部の企業ではバリューを出しています。

・インターネット広告配信:オレシカナイトでSpeeeの方がトンプソン抽出で精度を増していた。
・推薦システムにおけるコールドスタート問題:ネットフリックスが情報推薦の際にContextual Banditを適用

バンディット問題とは異なるものの、最適腕識別問題においては、クックパッドのクリエイティブ出し分けやGoogleのウェブテスト(旧Webサイトオプティマイザー)などで使われています。ちなみに、バンディット問題と最適腕識別問題は似て非なるものであるということを『バンディット問題の理論とアルゴリズム』で知りました。

また、マーケティングとは違いますが、株価のトレーディングの際にバンディットアルゴリズムを使っているという事例(Bandits and Stocks)が当然ながらあるようです。

バンディット問題で出てくる数学的な知識と方策

バンディット問題の書籍を読もうとすると、数理統計学の知識が必要です。

あるスロットを何回引くべきかという意思決定の際に、「神のみぞ知る真の報酬」と「あるスロットの報酬」がどれくらい外れているか、そしてそのハズレ具合は許容できるのかということが重要になります。
「神のみぞ知る真の報酬と、あるスロットの報酬がΔだけ外れている確率」の推論の精度に関心があるということです。

バンディット問題において、「その時のベストのスロットを引いた際のリターン」と「その時実際に選んだスロットのリターン」の差の期間合計値をリグレットとして、そのリグレットを小さくするようにスロットを選びます。
そのリグレットに対して理論的な下限を求める際に、数理統計学の知識が必要になります。

具体的には、ヘフディングの不等式、その前提となるマルコフの不等式やチェビシェフの不等式やチェルノフ限界、積率母関数やイェンセンの不等式などです。
それらを駆使しながら、様々な施策の中で、理論的な下限がより小さくなるようなものを探そうという流れのようです。

『バンディット問題の理論とアルゴリズム』を読む上で前提となっていそうな知識として、スタンフォード大学の講義資料(CS229 Supplemental Lecture notes Hoeffding’s inequality)を運良く見つけることが出来たので、これをもとに学ぶと理解が捗ると思います。

リグレットの下限を低めることを目指して、様々なアプローチが議論されます。

ε-貪欲法

  • 概要:スロットを回す回数のうち、一定割合(ε)をスロットの探索に当て、残りの期間を良いとされるスロットを回し続ける。
  • メリット:実装が容易でシステムに組み込み易い
  • デメリット:期待値が悪いスロットも良いスロットも同じ回数引いてしまうので性能が悪くなる。スロットの種類が多い際はより一層悪くなりやすい。

UCB(Upper Confidence Bound)方策

  • 概要:標本平均に補正項を足した、UCBスコアを各時点ごとに計算し、最もスコアが高いスロットを回す。なお、補正項は選択回数の少ないスロットに対して大きくなります。
  • メリット:ε-貪欲法と異なり、リグレットの上限がεなどの水準に左右されない。ハイパーパラメータが少ない。
  • デメリット:真の期待値についての信頼区間を求めることは本質的ではない。

KL-UCB

  • 概要:KLダイバージェンスを用いてUCBスコアを計算し、最もスコアが高いスロットを回す。
  • メリット:KLダイバージェンスを様々なモデルに応じて置き換えることができるなど、柔軟性がある。
  • デメリット:KLダイバージェンスの逆関数を計算する必要があり、毎回ニュートン法などを適用する必要がある。

MED(Minimum Empirical Divergence)方策

  • 概要:期待値最大である際の尤度が一定以上のスロットを回すという方策。
  • メリット:KLダイバージェンスの逆関数を計算する必要がない。
  • デメリット:KL-UCBよりも性能が悪い。IMEDという方策であればその弱点を克服している。

トンプソン抽出

  • 概要:期待値最大でないスロットの選択数の期待値を近似的に最小化するという取り組みを、ベイズ統計の枠組みで行ったもの。
  • メリット:経験的に高い性能となりやすい。
  • デメリット:?

Contextual Bandit問題とは

ある時点のあるスロットの報酬が、ユーザーの特徴量と誤差項により線形で表すことができるものを、線形バンディットと呼びます。
ユーザーの各行動の特徴量が時刻により異なる値を取ることを許すという設定を、文脈付きバンディット(Contextual Bandit)と呼びます。
つまり、Contextual Banditは時刻により異なるユーザーの特徴量が与えられたもとでの、利得の期待値の最大化問題となります。

具体的には、パチンコ店における期待値最大化の行動を考えるとすると、パチンコ台の大当たり確率は、午前か午後か、大当たりが既に他の台で出たか、その台がどれくらい回されているかなどの時間による文脈に左右されるという状況となります。

このContextual Banditにおいても、先程あげたようなリグレットを最小にするような様々な方策があります。LinUCB方策や、線形モデルのトンプソン抽出、ロジスティック回帰モデルのバンディットなどです。

Contextualパッケージでできること

こちらの資料にある通り、バンディットアルゴリズムのシミュレーションとオフライン評価が行えるパッケージです。
多様なバンディットアルゴリズムを試すことができます。
要となるデータですが、シミュレーションにより生成することもできれば、過去にランダムに出し分けたログなどのデータがあればそのデータをもとにアルゴリズムの検証をすることができます。

サンプル実行

さて、今回は完全に手抜きです。GitHubにあったサンプルコードを3つほど回すだけです。ただ、特徴量の突っ込み方などをサンプルコードから学べるので、ぜひ開発者のGitHubをご覧ください。

サンプル1:ABテストによる最適腕選択

パッケージのGitHub
にコードがありました。Bandit Algorithms for Website Optimizationという書籍に登場してきている例をRで実行できるサンプルです。
・ε-貪欲法を様々なεでシミュレーションして最適なスロットを見つける
・ソフトマックスによる方策に関しても様々なτに応じたシミュレーションをして最適なスロットを見つける
・UCB方策によりシミュレーションを行い、最適なスロットを見つける。ε-貪欲法やソフトマックスとの比較を行う
という実験ができます。シミュレーションの設定として、スロットごとの当たりの出る確率をベクトルで指定しています。

実行するのに10分くらいはかかるかもしれません。

ε-貪欲法

・最適なスロットを選んだ確率

・平均報酬額

 ・累積報酬額

ソフトマックスによる方策

・最適なスロットを選んだ確率

・平均報酬額

・累積報酬額

UCB方策

・最適なスロットを選んだ確率

・平均報酬額

・累積報酬額

サンプル2:文脈付きバンディット問題で映画のレーティングの最適化

同じGitHubにあるこちらのコードは、映画のデータセットに対して、文脈付きバンディット問題でオフラインテストをするためのコードです。映画のレーティングが4以上なら1そうでないなら0のデータを作り、特徴量として映画館で見たか家で見たか、一人で見たか家族と見たか、週末に見たかどうかなどの変数を7個ほど作成しています。方策としては、ランダムなもの、ε-貪欲法、トンプソン抽出、LinUCBをシミュレーションしています。

実行してから処理が止まるまで1時間程度はかかりましたが、LinUCBが累積の報酬が大きいようです。

サンプル3:文脈付きバンディット問題でMovieLensのTop50の作品における評価の最適化

こちらのコードは、MovieLensのデータセットにおいて、特徴量として過去にユーザーが評価した映画のカテゴリーの割合を19カテゴリ分用意して、ユーザーの見た映画の評価を最も高めるという、文脈付きバンディット問題です。こちらは実行して、30分程度で処理が終わりました。先程のサンプルと同じで、LinUCBが累積の報酬が大きいようです。

おわりに

2~3年前に、Tokyo Web Miningの懇親会でContextual Banditの論文いいぞとテラモナギさんが紹介していて、へー、そんなのあるんだと、「へー」の域を出なかったんですが、一歩前進した気がします。先人が切り開いた道を2~3年後に舗装されてから通るというのも遅いなと感じられるので、残業もっと減らして勉強時間増やしたいと思います。

参考情報

バンディット問題の理論とアルゴリズム (機械学習プロフェッショナルシリーズ)
Bandit Algorithms for Website Optimization: Developing, Deploying, and Debugging
Contextual package ~ Japan.R Shota Yasui
Package ‘contextual’
バンディットアルゴリズムの復習3:UCB(Upper Confidence Bound)

Rでオペレーションズ・リサーチ(OR)に関する情報をあさる/ コード付き

はじめに

私は基本的にデータ分析を生業としていますが、どうしても分析の案件が足りない時期は分析以外のものに手を染めることもあります。主に、RPAやクローリング、APIを用いたソーシャルリスニングなどです。今後も分析以外のことをやる時があるとしたら、レパートリーを増やしたいですよね。なので、ORについて調べてみることにしました。

ORとは

公益社団法人日本オペレーションズ・リサーチ学会による定義によると、

「現象を抽象化した数理モデルを構築し, モデル分析に基づいて種々の問題, とりわけ意思決定問題の解決を支援する方法論や技法の総称. 情報化社会の進展に伴って, 線形計画法に代表される最適化モデルや待ち行列理論に代表される確率的なモデル等, 多様なモデルに基づく分析が, 経営計画や生産・販売・財務等の企業意思決定や都市・公共システム等広く社会一般の問題解決に大きな役割を果たしている.」

とされています。うむ、お硬い感じの定義ですね。
他の記述にわかりやすい表現がありました。

「問題を科学的,つまり「筋のとおった方法」を用いて解決するための「問題解決学」であります」
これならわかりやすいです。

問題解決につながる、あらゆる科学的な手法を扱っているのがORだと考えてよいのだろうと思います。

ORの強みとしては、

  • 大規模プロジェクトなどの遂行に役立つ
  • 常識や過去の経験では判断が難しい問題に対する解の提供をしてくれる
  • 経営、工学、医学、公共政策など幅広い分野での適用可能性がある

などがあげられています。

ORの手法

問題解決につながるなら何でもありということで、手法も幅広いようです。私の持っている参考書やOR学会のサイトの情報から判断すると少なくとも以下の手法が扱われているようです。

・数理最適化
・組合せ最適化
・シミュレーション
・待ち行列
・AHP(階層的意思決定法)
・DEA(包絡分析法)
・スケジューリング
・ゲーム理論
・ネットワーク理論
・データマイニング

データマイニングはもはや機械学習ブームなので特筆したものではないですが、最近データ分析を始めた人などは知らないことも多いのではないでしょうか。

仕事やプライベートでの使い所

  • 社員のシフトを決めるタスク
  • ミーティングに参加する社員の移動距離が一番少ない、空いている会議室を見つけるタスク
  • コンビニのレジを何台おけば客の待ち時間が想定内になるのか
  • 工程A、B、Cの全てを経る必要のある作業で、生産ラインを安定させるには今日はAとBとCのどの工程をどれだけすすめるべきか決めるタスク
  • 旅行をする際に、予算を所与のもとで、どのスポットに立ち寄ることが効用が高いかを見つけるタスク

Rでの実践例

ようやく、この記事の本題です。Rでオペレーションズリサーチなどという書籍に出会えていないので、Rを用いてオペレーションズリサーチを行っている事例を集めてみようと思います。ブログを漁ればいろいろとありますね。

なお、答え合わせを兼ねて、登場する問題は『Excelで学ぶOR』の例題で表現を変えています。網羅性はないものの、出来るだけ取り上げてみようと思います。

数理最適化

シュークリーム専門店、pseudoカモノハシはシュークリームとパンケーキの生産をしているが、厨房があまりに狭すぎてシュークリームとパンケーキの同時生産ができない。また、労働基準法の観点から厨房の利用時間は40時間以内となる。

シュークリームの生産に関しての情報は以下の通りとします。

pseudoカモノハシはシュークリームとパンケーキをどれだけ生産すれば利益を最大にすることができるだろうかという問題です。

定式化すると以下のようになります。

RのlpSolveパッケージを用いてこの線形計画問題を解いてみます。

実行するとこんな結果です。

pseudoカモノハシは25個のシュークリーム、15個のパンケーキを生産することで615万円の売上をあげることができるということになります。

係数がどれくらい変わっても最適解が変化しないのかを知るための感度分析も簡単にできるようです。
係数がちょっと変わっただけで崩れる最適化とかだと実務で使う際に怖いので、大事な工程ですね。

輸送問題(ネットワーク型の線形計画法)

シュークリーム専門店、pseudoカモノハシは事業拡大につき、3つの食料庫と4つの工房を持つに至った。シュークリームやパンケーキを生産するためには食料庫から工房までトラックで輸送をする必要がある。各々の食料庫から工房までの輸送コストは以下の表の1行1列目〜3行4列目までで表される。食料庫には置ける在庫が決まっており、表の5列目で与えられている。工房には客の注文ベースの生産ノルマが課されており、表の4行目で与えられている。

以上、pseudoカモノハシは輸送コストを最も下げて生産するにはどの食料庫からどの工房に材料を輸送すればよいか、と言う問題となります。

定式化すると以下のようになります。

このコードを実行すると以下のようになり、輸送コストの最小値は3610であることが示されている。

ナップサック問題

シュークリーム専門店、pseudoカモノハシは創業3周年記念としてご当地グルメとのコラボレーションを計画している。グルメ系コンサルが提案したプランA〜Jの10件のご当地グルメとのコラボには、それぞれ費用と便益がある。この施策への予算が2000万円だとして、pseudoカモノハシは総便益が最大になるようにどのコラボ企画を採用するべきか。

定式化すると以下のようになります。

こちらにRのコードがあったので拝借しております。

これを実行すると、以下の結果が得られます。最大便益が31であること、それを実現するプランとしてE、G、Hが選ばれることがそれぞれ示しています。ただし、31を実現する解は他にもあります。

混合整数計画

シュークリーム専門店、pseudoカモノハシは傘下のパンケーキチェーン店、ヒッグス・シングスの3店舗から、「繁忙につき代替生産をお願いしたい」と言われた。ヒッグス・シングスのパンケーキは3店舗それぞれ味が異なり、それによってパンケーキの製造コストなども異なる。pseudoカモノハシの工房でパンケーキは生産可能ではあるが、そこそこに忙しいので生産できて15個が限度だと考えられる。そこで、以下の表が与えられたもとで、pseudoカモノハシはヒッグス・シングスの3店舗それぞれのパンケーキをどれだけ生産すれば利益が最大化されるか。

定式化すると以下のようになります。

今回は調査の結果、Rglpkというパッケージがあることがわかったので、そちらを用います。

このコードを実行すると、以下の結果が得られます。工房1と工房2で8個と6個の生産にコミットすることで、10万円の最大利益を実現できることが示されています。

ウェーバー問題(非線形計画)

ウェーバー問題というのは、ORWikiによると以下の定義とされています。

施設・顧客間の距離に需要量を乗じたものの総和を最小化するような単一の施設の配置を、 平面上の任意の地点の中から決定する問題。

具体例で取り組んでみましょう。

シュークリーム専門店、pseudoカモノハシは新しく食料庫を設置したいと考えている。店舗の位置が(x,y)座標で与えられているものとし、各店舗からの距離が最も小さくなるような位置に食料庫を設けるとしたらどこになるか、という問題。

定式化すると、以下のようになります。

orlocaというパッケージを使えばウェーバー問題のような非線形計画問題を簡単に解くことができます。

このコードを実行すると以下の結果が得られます。x座標が4.468528、y座標が3.843755の時に、距離の最小値が22.80259となることが示されています。

せっかくなので、いらすとやの画像を背景に結果などをプロットしてみます。赤い点が解となった点です。

待ち行列

シュークリーム専門店、pseudoカモノハシ本店における来店客の行列に関してシミュレーションするものとする。

  • pseudoカモノハシ本店は平均して1時間に50人来店し、それはポワソン分布に従うとされている。
  • pseudoカモノハシ本店の商品は商品ごとに提供するまでにかかる時間が異なり、客はどれか1品を選択するが、その選択確率がおおむね決まっている。(以下の表)

進め方としては、まず、一様乱数を生成させ、ポアソン分布の分布関数の逆関数を求め、その逆関数に乱数を入力することで顧客の来店時間の間隔を生成します。

次に、顧客の選択する商品をシミュレーションするために別で一様乱数を生成し、その乱数の取る値に応じて商品を割り当てます。
そして、商品の提供に時間がかかることから、開始時間に商品の提供時間を足して、終了時間を求めます。ただし、来店時間に前の顧客がまだ商品を受け取れていないと、待ち時間が発生するので、来店時間の時点で終了していない場合はその分だけ開始時間が遅れます。

以下のRコードを作成してみました。

まず、来店間隔のシミュレーションですが、以下のようになります。

続いて、来店時間に応じた、待ち時間の推移です。

来店者の到来とともに、ぐんぐんと伸びているのがわかります。私はシュークリームに50分も待てないですね。

続いて、任意の顧客が開始した時に、すでに待っている客の数の推移です。

最後に、空き時間です。

顧客が忍耐強く待つのであれば、ほとんどレジは休めていないという劣悪な労働環境になりそうですね。厨房の能力やレジの能力などを高める必要がありそうです。

最短路問題

シュークリーム専門店、pseudoカモノハシのオーナーSKUE氏は、いま店舗1にいるが、店舗6の店長との1on1があるため、向かおうと考えている。ついでに他の店舗の店長にも顔を出したいと考えており、他の店を経由して一番距離が短い経路を見つけたい。店舗と店舗の距離は以下のグラフのノード間のラベルで与えられているものとする。

Rのigraphパッケージに「Shortest (directed or undirected) paths between vertices」というグラフ間の最短経路を見つける関数があるので、そちらを使います。この関数は最短経路問題を解く際に効率的とされる、ダイクストラ(Dijkstra)法というアルゴリズムを採用しています。具体的に使った例が載っているブログとしてこちらを参考にしています。

このグラフから、以下のようなデータを作っておきます。

このコードを実行すると、以下のように、店舗1→店舗3→店舗5→店舗6の順番で店舗を巡ると最短経路である13が達成されることが示されます。

今回は店舗が少ないので、人間の目でも見つけることができますね。

巡回セールスマン問題

カワウソ急便の配達員が今、バレンタインデーの集配のためにpseudoカモノハシの店舗1(本店)にいる。配達員は他の全ての店舗の集配もする必要があり、集配後は確認のために本店に立ち寄る必要がある。どのようにして店舗を一度ずつ回れば移動距離が最も小さくなるかに関心がある。店舗ごとの位置は以下の図の通り。

各地点の座標は以下の表で与えられている。

このような問題を巡回セールスマン問題と呼ぶが、 解く際は座標の点から各々店舗の距離を計算し、距離行列を作成し、その距離行列をもとに、1度しか通れないという制約条件を課しながら距離が最も小さくなる組み合わせを見つける。
RではTSP(Traveling Salesperson Problem)パッケージがあるので、proxyパッケージを使って距離行列を作ってしまえば、簡単に最適な順路を示してくれる。

このコードを実行すると、以下のようになり、1→2→3→4→5→1と巡回することで総距離が12.641で済むことがわかる。

これ以降は疲れ果てたので、実践というより紹介にとどまりますが、あしからず。

スケジューリングと集合被覆問題

Rでの事例が見つかりませんでした。Qiitaで以下のようなPythonによる実践があったので、それをもとにRで書き換えてみるのも良いと思います。(後日、追記したいと思います。)
組合せ最適化 – 典型問題 – 集合被覆問題
組合せ最適化 – 典型問題 – 勤務スケジューリング問題

NPV、IRR

ファイナンスで基本のNPV(Net Present Value)やIRR(Internal Rate of Return)の計算ですが、Package ‘FinCal’パッケージで簡単に計算できます。これくらい自分で書いてもいい気もしますが。このパッケージに割引率とキャッシュフローのベクトルを入力したらNPVを返してくれるnpv関数やirr関数があります。

データ包絡分析法

先日、ブレインパッドさんがこちらの記事で公開していたデータ包絡分析法(Data Envelopment Analysis:DEA)ですが、Rのコードが付いてなかったので漁ったところ、早稲田大学の逆瀬川先生がこちらでコードを公開しているようです。

DEAは同質な複数の事業体の相対的な効率性評価のための方法と定義されています。ブレインパッドさんの例ではコストを入力として、出力としての売上が効率的かどうかを見るために使われています。

パッケージに関しては、Data Envelopment Analysisでググったら、Benchmarkingというパッケージを見つけました。このパッケージは、少なくともブレインパッドさんのブログで紹介されている、DRS(Decreasing Returns to Scale)とFDH(Free Disposal Hull)は引数で選択可能のようです。

このコードを実行すると、ブログと似たような図が描けるようです。

詳しくはドキュメントについている先行研究とかを見たいところです。

ポートフォリオ選択問題

これもファイナンスで基本となる、ポートフォリオ選択問題なのですが、投資家のリスク選好から無リスク資産(国債とか)とリスク資産(株式とか)をどれくらいの配分で持つのが効率的かを解くという問題になります。

どうやらtidyquantというパッケージを使うことで、効率的ポートフォリオを探索できるようです。こちらのドキュメント(Tidy_Portfoliomanagement_in_R)をまだ実践できていないですが、最終的に以下のような効率的ポートフォリオを描けるようです。

編集を終えて

まだまだ原理を深く理解できていないですが、いざ仕事の依頼が来た際の取っ掛かりとしては良いものが手に入った気がします。加えて、人間だと計算が厳しそうな問題を解けるというのは非常に面白いです。
「問題設定→定式化→コードに書き落とす」という一連の訓練を続けるとかなり力が付きそうな気がします。Rって統計学・機械学習以外にも本当に幅広く取り揃っていて飽きがこなくてよいですね。

参考情報

Excelで学ぶOR
サルでもわかる待ち行列
RでLinear Programming
問題解決の数理(’17)
Sensitivity Analysis 感度分析 もし○○○だったどうする?
Rでデータ解析を始めよう020 Rでナップサック問題を解いてみよう
CRAN Task View: Optimization and Mathematical Programming
Rで数理計画
RでOR:待ち行列モデル
[R]ggplot2によるグラフィックスで、図にPNG形式の画像を貼る
Rでクラスター分析〜距離行列の生成からクラスタリングまで
RでTSPの練習

R Advent Calendar 2018 一発屋芸人の検索トレンドの分析

はじめに

昨年のR Advent Calendarはポケモンのデータをrvestでスクレイピングして、レアポケモンがどのような種族値における特徴があるのかを探ったり、経験値が必要な割に種族値が低い「コスパの悪いポケモン」などを見つけました。
今年のR Advent Calendarでは、年末年始といえば一発屋芸人のテレビなどでの露出が多くなることから、一発屋芸人の検索トレンドのデータを手に入れて分析してみたいと思います。

分析工程

・データの収集
・データの整形
・可視化
・分析

データの収集

こちらのサイト(流行した一発屋芸人一覧/年代流行)に一発屋の芸人さんが列挙されていました。私は普段テレビを見ないので大体の芸人さんがわからないです。

Googleトレンドから、芸人名に関するGoogle検索の時系列データを収集します。

非常に残酷なデータだなと思いました。

ただ、一つ弁護すると、Googleトレンドはレベルではなくピークを1として標準化した数値をデータとして提供してくれていますので、
ピークが著しく高ければ、今の水準が低くてもそこそこ検索されている可能性はあるとだけ言っておきます。

本当の検索回数が必要な場合は、Google Adwords(検索連動型広告)のアカウントの開設とともに検索ボリューム取得APIなどの申請が必要なので、正確なデータが必要な場合は会社として取り組んだほうが良いと思います。個人では厳しいです。

データの整形

各芸人さん(総勢21名)の検索トレンドデータのピークの6ヶ月前までのデータと6ヶ月後のデータまでの合計1年間の検索トレンドを各々抽出してみようと思います。
GoogleトレンドのデータはCSVでダウンロードできますので、そのCSVを読み込み、トレンドのデータを文字列から数値にし、ピークの前後12ヶ月ずつのデータを抽出します。
そうすることで、一発屋芸人のピークの前後に関するデータを作ります。(ただし、今朝、Google Trendのデータを取得できるgtrendsRというパッケージがR bloggerで紹介されていました。APIないはずなんですが、URLの工夫か裏でSelenium動かしていたりするんですかね。)

可視化

作成したデータを実際にプロットしてみます。

一発屋にも盛り上がり方に違いがあるようですね。

時系列クラスタリングの適用

多様な盛り上がり方があることから、TSclustというライブラリを使って時系列クラスタリングを行い、トレンドに関しての分類的なものを得たいと思います。
今回初めて使うのですが、参考文献によると様々な類似性指標を指定して、時系列ごとの類似性を計算するようです。ピアソン相関係数のようなシンプルなものもあれば、ユークリッド距離のものやFrechet距離とかいう聞いたことないものまで幅広く用意されています。今回はシンプルにピアソン相関係数にしてみます。そして、類似性指標を出してから、そのまま階層クラスタリングを行います。

こちらが、TSclustのdiss関数を用いて計算した時系列データごとの距離を、階層クラスタリングにより描いたデンドログラムです。

この分類だけ見ても、芸人さんを知らない私からすると何も共感がありませんので、先程のクラスタリング結果をもとに可視化をしてみます。
そこで、Tokyo.Rで知らない人はいないであろう、yutaniさんの作られたgghighlightを使ってみようと思います。

ただ、日本語のラベルの表示がうまくいかなかったので、芸人さんの名前をGoogleSpreadSheetのGoogle翻訳関数(GOOGLETRANSLATE)で英訳しておきます。

(Anyway bright YasumuraやThick slice Jasonは結構キャッチーなのでは?)

まずはクラスター1

比較的短期でピークに達し、すぐに検索されなくなる、一発屋の名に相違ない傾向を持ったクラスターのように思われます。「日本エレキテル連合」とか「楽しんご」とか「8.6秒バズーカ」とかです。

続いてクラスター2

急激にピークに達するものの、ややしぶとく残り続けるような一発屋のクラスターなのかなと思います。「レイザーラモンHG」とか「厚切りジェイソン」とか「ピコ太郎」とか「世界のナベアツ」です。

そしてクラスター3

3人の芸人さんしか属していないですね。クラスターの数は2個でもよかったかもしれない。段階的にピークに達し、一気に落とされるという一発屋のクラスターのようです。「とにかく明るい安村」とか「藤崎マーケット」とか「すぎちゃん」とかです。

様々な傾向の一発屋さんがいるのがわかりました。

トレンドの推定

今回扱っているデータは芸人さんの数×時点のデータの多変量時系列となります。都合の良いものはないかと考えていましたが、古典的なVARではサンプルサイズ的にかなり苦しいと思い、Stanによるダイナミックパネルデータ分析などの事例はないか漁っていましたが、なかなかありませんでした。

松浦さんの『StanとRでベイズ統計モデリング (Wonderful R)』の241pに書かれている、モデル式12-8や12-9が今回のものに適しているなと思いましたが、コードを上げている方は見当たらなかったです。よしそれならば作ろうかと思った矢先、logics-of-blueさんのStan Advent Calendarの投稿、「Stanで推定する多変量時系列モデル」がかなりどんぴしゃな内容でしたので、コードを拝借してこの一発屋データの推定をしてみようと思います。

まずは、stanのコード

そしてキックして結果を可視化するためのRコード

そのまんま実行して、一発屋の時系列の中央値を可視化したらこんな感じになりました。一発屋のトレンドをうまく抽出できているのかなと思います。

今後の改良としては、階層性を持たせ、芸人さんごとのハイパーパラメータを持たせるとかなのですが、正月にでも取り組みたいと思います。(芸人さん以外のデータでやりたい。)

一方で、他にも多変量時系列で何かないか漁っていたのですが、Applied Time Series Analysis for Fisheries and Environmental Sciences : Dynamic factor analysisで紹介されている、Dynamic Factor Analysisというものが面白そうだなと思いました。
bayesdfaというパッケージを用いて、多変量時系列データに存在するであろうトレンドをStanを用いて推定することができるようです。元となった論文には各エリアごとのノルウェーロブスターの個体数のトレンドを推定し、3つのトレンドが発見されたとしています。ただ、同時点間のデータではないという点から今回のデータへの適用は不適切です。

同時点間に観測されていないデータであるという問題を認識した上で、このパッケージを使ってどんなトレンドを抽出できるのか試してみようと思います。

徐々に増えてから一気に落ちるトレンドや、一気に増えてから徐々に落ちるトレンドなどがうまく捉えれている気がします。
さらなる試行として、AICのような情報量基準である、Leave One Out Information Criterion (LOOIC)が最も低くなるトレンドの数を探索してみます。

トレンド数を1から5まで指定して実行した結果、5の時が一番LOOICが低くなりました。

まぁ、適切な使い方ではないのですが、徐々に増えてから一気に落ちるトレンドや、一気に増えてから徐々に落ちるトレンドなどが引き続き捉えれているようです。

今後の課題

・Stanによる多変量時系列のモデリングをしてみる。(Dynamic Panel分析とかもできると良い。少なくともStanのドキュメントにはない。)
・Dynamic Factor Analysisの適切な事例での適用をしてみる。

それでは、どうか良い年末をお過ごし下さい!
メリークリスマス!

参考情報

Introduction to gghighlight: Highlight ggplot’s Lines and Points with Predicates
{TSclust} ではじめる時系列クラスタリング
Applied Time Series Analysis for Fisheries and Environmental Sciences 9.7 Dynamic factor analysis
読了:Montero & Vilar (2014) RのTSclustパッケージで時系列クラスタリング

参加できなかった第74回TokyoRのキャッチアップと結婚式について

自身の結婚式があったため、参加できなかった第74回目のTokyoRについてキャッチアップするために、公開資料を読んだ際に感じたメモをここに記します。
そして、世のデータサイエンティストが結婚式の際に苦しまないように参考になりそうな情報を少し書きました。

今回はTogetterも初めて作ってみました。(第74回R勉強会@東京(#TokyoR)のタグが付いたものたち

初心者セッション

初心者セッション1 – Data Import & Export –

資料なし

初心者セッション2 – Data Handling –

https://ymattu.github.io/TokyoR74/slide.html#/

いつもながら良い資料です。まだ社内でアクティブにRを広めてはいないですが、広めるならこの資料が良いですね。
上から下に直感的に書けるのは初学者には大事なので。
ただ、SQLの知識がないメンバーとかには補助教材が必要ですね。

登場するパッケージはdplyrとlubridateとstringrとforcatsとpurrr。

初心者セッション3 – Plot & Visualization –

資料なし

応用セッション

How LINE Corp Use R to Compete in a Data-Driven World

資料なし

LINEでのRの活用最前線の話なのでしょうが、資料がないのが寂しいです。
きっと、タイムラインで共有されていたこの記事をベースに話されているのではないか。

LINE の全社員が必要に応じて担当サービスのデータを分析できる環境の構築

  • 総勢50名の機械学習エンジニア・データサイエンティスト・データプランナー・データエンジニアが活躍できる分析基盤を作っていますよと。
  • Hadoop クラスタのデータの全社公開&活用のために、エンドユーザ向けのWebインターフェース「OASIS」をゼロから新規に開発。
    Spark, Spark SQL, PySpark, SparkRおよびPrestoクエリを叩くことができる。
    約20のサービス・部署で利用され、月間利用者数は約200人(データ関連の人以外で150人は触っていることになる。すごく層が厚そう。)

あと、ブリスベンの写真がTwitterで写り込んでいたので、UseR!2018のお話をされているのだろうと思われます。
linerパッケージについても語られている模様。

UseR!2018に参加し、社内Rパッケージ「liner」の活用事例を紹介しました

これのことでしょう。写真も合致している。
「いらすとや」の画像がワールドワイドに使われているのがシュールでいいですね。

「データの取得、分析、レポーティング、そして結果の共有にいたるまで、様々な便利機能を提供」と書かれています。
図を見る限りは、

  • RStudioでの解析結果をDBにカジュアルに保存したり、通知したり、ドキュメント化したりできる
  • PrestoやらHiveQLなどもRStudio上で実行できる(Tab補完とかもしてくれるんでしょうか?そこは聞いてみないとわからない。)
  • コーポレートカラーを適用したggplot2を利用できる
  • A/Bテストの結果をShinyのアプリで確認できる

などの機能があるようです。50人もいる高単価な人達の時間を節約できるという点でも、すごくインパクトのある取り組みですね。

tidyeval入門以前

(speakerdeckのembedに若干苦戦しましたw)

湯谷さんの考えるtidyevalについて英語で書かれています。
環境に応じた値をRは自動で引っ張ってくるけど、たまに干渉してしまうことがあり、実際dplyrとstatsは共にfilter()関数を持っていると。
どの優先順位をもたせるかのコントロールって難しそうですね。
干渉を避けるためにquosureやunquoteというのが説明されています。私の理解が追いついていないので、正直めちゃわかった感はないのですが、
様々な環境下で動くパッケージを作ろうと考える際は不可避な領域なのではないでしょうか。

LT

「うまい飯が作りたい」

recipesパッケージの紹介です。面白い方なんですね。
tidymodelsパッケージの中に内包されているそうな。

recipesパッケージの使い方を丁寧に説明してくださっています。

目的変数と説明変数を最初に明示的に指定し、前処理の手法をパイプ演算子でつないでいくだけ。
前処理のステップは使い回せるとのことで、似たようなデータを扱う場合は使いまわして楽をできるとのこと。
これは試してみる価値がありそうですね。k-nn法を用いた欠損値補完なども関数として用意されているようです。

マジレスすると、モテるかどうかは容姿や性格によるところが大きいと思うので、
身なりを清潔に保つとか、連れて行くと喜ばれそうな場所や体験を提供するとかが近道な気がしますよね。

不連続回帰とrdrobustパッケージの紹介

計量経済学系のバックグラウンドをお持ちのkagglerの方のLTです。
イスラム教の社会が、女性の高校修了率に与える影響について、統計的因果推論をされています。
rdrobustというパッケージを用いて、不連続回帰(RDD)という手法を実践され、イスラム教政治がランダム割当になりやすい状況を作っています。
全データで推定すると負の影響が推定されていましたが、RDDを使うことで正の効果があるという結果となりました。これは他の学部領域での先行研究とも整合的とのこと。

PCAや対応分析で補完要素を使う

資料なし

CiNii API その2

資料なし

Rで健康体

資料は後ほど公開される模様

データのみは公開されている


https://github.com/weda-654/my_health_log
0 forks.
0 stars.
0 open issues.
Recent commits:

Soccer × Attribution Analysis

u++さんのLTですね!アクセスログや第三者配信データを用いた分析でよく扱われるアトリビューション分析をサッカーの貢献度に利用したという話です。
ChannelAttributionパッケージを使われています。

マーケティングの業務で使う際は、これまで見逃していた意外なページや意外な参照元などがこの分析で見えてきたりします。
ただ、価値はあってもコストがかかっては元も子もないので、コストに関する記述があったのも実務で使われている方の視点だなと思いました。

サッカーはルールくらいしか知らないので誰が意外なのかはわからないですw

地理空間データの交差検証、正しくできていますか?

地理空間データにおける交差検証する際の手法として、Spacial Cross-ValidationとTarget-oriented cross-validationが挙げられています。
空間データ向けのパッケージとしては、sfやCASTが、学習周りではmlrやcaretが扱われていました。
Referenceがあるのが嬉しいですね。

reticulateパッケージとデータサイエンスフロー

資料なし

結婚式について

なぜ開くのか

大事な思い出づくりのため。

工程

  • Willing to Payの決定(全てはここ)
    • 持ち込みし放題の式場
    • 大学OB割引などがある式場
    • 料理が美味しい(国賓を迎えたりしているか?)
    • 荘厳な雰囲気
    • アクセスが良い
  • 教会式か人前式か神前式かの選定(教会式だとめっちゃ高かったりするところもある)
  • 誘う対象の選定
    • 共に勉強を頑張った仲
    • 共に仕事を頑張った仲
    • 親族
  • 住所の聞き出し
    • Googleフォームを活用
  • 大量の切手の購入(送付用と返信用も)
    • 郵便局は21時まで空いていたりするので助かった。(どこもそうなのかな?)
  • 上司や友人へのスピーチの依頼、乾杯の依頼
    • どういうオーディエンスなのかを事前に伝えておく。スピーチ作成者の負担を軽減する。
  • 招待状の作成、発送
  • 料理の試食
    • いろんな種類を食べるのでお腹ぱんぱん。
    • 料理は1000~2000円程度の予算アップはした方が面白い。
  • ウェディングケーキのデザインのすり合わせ
  • テーブルクロス、花、引き出物、ネームプレートなどのすり合わせ
  • 司会者との打ち合わせ
  • 自分のスピーチの作成(新郎新婦ともに)
    • 意外性と感動と笑いを織り交ぜるのが良い。
  • メッセージカードの作成
    • スプレッドシートに書きなぐり、それを手書きでひたすら書ききるのみ。
  • 式場音楽の選定
  • イベントの進行のすり合わせ
  • DVDの作成(オープニング・プロフィール・エンディング)
  • 大量のピン札の調達(交通費は全額支給しろという親の教え)
  • ウェルカムボードの作成(ダイソーで4~500円くらいで材料は手に入る)

工夫した点

  • DVDを3枚自作した
    • お願いすると結構お金がかかるし、理想形に近づけるための試行錯誤の回数も限られる。
    • ハイスペックPCを持っているデータサイエンティストなら動画の編集に耐えられるはず。
    • マックのiMovieを使えばGUIで簡単に動画作成ができる。
  • 十分な大きさのメッセージカード
    • わざわざ来てくれた友人との思い出をとにかく書きまくる。

結果としてよかったこと

  • DVDの自作
    • 3枚とも笑いを提供でき、上映後は拍手している人もいたようです。(裏手にいたので会場の様子は直接見れなかったですが)
  • 定型文を避けたスピーチ
    • 書く前によくある定型文を見たんですが、面白くないと思い、独自のスピーチにしました。結果として大爆笑を提供できました。
  • 普段から写真や動画を撮りまくっていたことで、DVD制作の素材が潤沢にあった。
    • Googleフォトに写真をアップしまくれば容量に悩まずに済むので、とにかく日常的に撮り続けましょう。

やはり、結婚式はある意味でエンターテインメントなので、自分が工夫できるところ(スピーチとムービー)は少し頑張ってみるといいのかな、と思いました。

反省点

  • 結婚式前日は有給休暇を取るべき
    • 仕事を19時であがったけど、その後の準備で疲れ果てた。目にクマが若干できた。

結婚を頑張るエンジニアやデータサイエンティストの皆さんへ

結婚式の準備は色々と時間がかかります。
世の優秀なデータサイエンティストが、そのようなことに時間を割きすぎるのは社会的な損失なので、この参考情報を元に少しでも楽に準備をしていただけると幸いですね。

rstanarmパッケージを使って簡単にベイズモデリングを実行する

はじめに

今回は、rstanarmというパッケージを用いて赤ワインデータを色々といじってみようと思います。
マーケティングの意思決定のための分析などでベイズ統計を使う場面が多々あるのですが、似たような属性のデータがあるのであれば、
一つ一つstanコードを書くのではなく、Rの関数でサクッと実行して試行錯誤していくという形に持っていけたらいいなぁと感じていました。
本気を出すところではstanを、ルーティンワーク的なタスクではrstanarmをみたいな形で使い分けれると良いのではないでしょうか。

rstanarmとは

バックエンドの計算をStanに実行させて、統計モデルの推定を行うためのパッケージ。R上でlm関数のように簡単にベイズ推定を行うことができる。対象ユーザーはベイズ推定に慣れ親しんでいない頻度主義系のソフトウェアユーザー。
詳しくはこちら

インストールする

まずはrstanarmのインストールするのですが、コケまくりました。そのため、バージョンを下げてみることにします。

ここに過去のバージョンがありますが、2.17.4だと動かなかったものの、2.17.3なら動きました。

rstanarmのサンプルを回してみる

今回は、以下の文献を参考にして、大人のIrisとも言える、ワインデータを扱い、質の高いワインかどうかを決める要素を探ります。
How to Use the rstanarm Package | Jonah Gabry and Ben Goodrich

こちらの文献には、ベイズ分析の4つのステップとして以下があげられています。

  • 1.同時分布の特定(同時分布は事前分布と条件付きの尤度をかけ合わせたもの。)
  • 2.MCMCで事後分布を描く
  • 3.モデルがフィットしているか評価する
  • 4.事後予測分布を描き、結果に影響を与える予測項を確認する。

これらのステップをできるだけ素早くできると良いですね。

まずはデータを読み込んで、スケーリングしておきます。(可視化結果は前回と同じなので、載せません。)
加えて、6点以上の評価であれば1を取る二項変数を作成しておきます。

 

GLMでロジスティックモデルを推定し、rstanarmで推定した結果と比較します。rstanarmでは傾きや切片の事前分布にスチューデントのt分布を、尤度にロジスティック分布を設定しています。