Bayesian Statistics and Marketing – 混合ガウス×階層モデルのマーガリン購買データへの適用

前回の分析では、価格への反応係数の事前分布が正規分布を仮定したモデルを用いていましたが、事後分布から多峰性が観察されました。そこで今回は、各個人の価格への反応係数の事前分布が混合ガウス分布に従うとした場合の事例を扱いたいと思います。

データのおさらい

データ自体は前回のブログと同じですが、先日のTokyo.Rで松浦さんがオススメしていたGGallyパッケージのggpairs関数を用いて、今回扱うデータを可視化してみます。

まず、購買したマーガリンのブランド選択(6ブランド分)と、購入した価格(ドル)からなるデータセットの可視化をすると以下のようになります。

1つ目のブランドが最も多く選択されているようです。購入価格の分布は、ブランドによって多峰性がありそうです。

続いて、家計ごとの属性(家族構成、学歴、職位、退職の有無など)からなるデータセットの可視化をすると以下のようになります。

ほとんどダミー変数なので面白みには欠いていますね。ホワイトカラーの家計が多く、退職していない家計が多く、学歴が低い人が多いようです。年収に関しては対数正規分布に従ってそうです。

モデル

前回と同様に、「Bayesian Statistics and Marketing」の5章に載っている混合ガウスを想定した階層ベイズモデルを扱います。

yi ∼ MNL(Xi, βi) は購買レコードごとの意思決定が多項ロジスティック回帰モデルに従うということを意味し、βiは(説明変数の数×1)のベクトルとなります。
βi = Z∆[i,] + ui は購買レコードごとの価格にかかってくる係数で、その係数が家計の属性データに係数Δ(家計の数×属性データからなる説明変数の数)をかけ合わせたものと潜在的な項の和となります。なお、ここでの属性データには定数項を含めていません。係数Δは平均deltabar、分散A_delta^-1の多変量正規分布に従います。一方で、潜在的な項は平均µind、分散Σindの多変量正規分布に従います。この平均や分散に振られているindが、混合正規分布のパラメータの通し番号となり、多項分布に従います。この多項分布の割当確率がディリクレ分布に従います。最後に混合正規分布の各パラメータは、平均は正規分布に分散は逆ウィシャート分布に従います。

前回との大きな違いは、価格の係数の一部である、潜在的な項において、混合正規分布が仮定されているところになります。
ちなみに、モデルに関しての詳細はbayesmパッケージのマニュアルの61ページ目に記載されていました。

今回のモデルのDAG

Pythonのdaftというモジュールを使うことで、非常に簡単に今回のモデルのDAG(有向非巡回グラフ)を描くことができます。

今回はこちらのPythonコードで描けました。

stanコード

今回扱ったstanコードとなります。誤りがある場合はお知らせしていただけると幸いです。

stanをキックするためのコードです。先人が書かれた混合ガウスのスクリプトをHMCで実行した際に26時間ほどかかったので、今回はより複雑なモデルであることから、変分ベイズ法による推定を行ってみることにしました。松浦さんの教科書にあるように、vb関数を用いて変分ベイズ推論を行っています。

実行結果

今回はK={1,3,5,10}の混合要素数で推定を行いました。

K=1のケース

K=3のケース

K=5のケース

K=10のケース

Kが小さいと散らばりが比較的小さそうに見えます。逆に、Kが大きくなると散らばりが出てくるようです。

おわりに

 2回に渡って、マーケティングデータを用いたベイズ統計モデリングを学んできましたが、数式からstanコードに落とし込む作業の際に、stanの関数をある程度知らないとやりにくいという至極当然なことを実感しました。これまでに使ってきたstanの関数は限られたものしか扱っていなかったと言えます。特に今回は離散パラメータをstanで扱うパーツがあったので、先行研究や松浦さんの本を読みながらの手探りが多かったです。
 あと、複雑な階層ベイズモデルを扱う際に、頭の中を整理しないと手が止まってしまう感じがあったので、数式と対応するコードを横に並べながら進めました。
 マーケティングにおいては、顧客の属性ごとに多峰性のあるような事例を扱うことが多く、かつ各々のサンプルサイズも期待できないことが多いので、今回の学びを分析業務で試してみたいと思います。

参考情報

Bayesian Statistics and Marketing (Wiley Series in Probability and Statistics)
Package ‘bayesm’
Multivariate Gaussian Mixture Model done properly
StanとRでベイズ統計モデリング (Wonderful R)
機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)
daftでグラフィカルモデル

Bayesian Statistics and Marketingの5章 – 家計の異質性を考慮した階層ベイズモデル

はじめに

ゴールデンウィークで実家に持ち込む本としてチョイスしたのが、2005年出版の「Bayesian Statistics and Marketing」です。大学院のときに購入して、ちょっとしか読んでませんでした。

この本は、字面の通りマーケティング関連の分析に関してベイズ統計を使ってアプローチするというもので、この書籍のために作られた、Rのbayesmというパッケージの紹介もあり、理論だけでなくRで実践することもできます。1章から7章までの全ての分析事例に対して実行可能な関数が用意されています。(CRANにあるdocumentも120p程度と割と大きめのパッケージです。)

和書で言うと、東北大学の照井先生の「ベイズモデリングによるマーケティング分析」などがありますが、その82pでもBayesian Statistics and Marketingとbayesmパッケージが紹介されています。

今回は、5章に載っている階層ベイズモデルを用いた、家計の異質性を考慮したブランド選択モデルの分析を紹介します。加えて、GitHubでstanによる再現を試みている方がいらっしゃったので、その方のコードの紹介も行います。

最近はこれまで以上にベイズ統計が流行ってきていますが、マーケティング×ベイズの書籍は限られている印象なので、少しでもリサーチのお役に立てれば幸いです。

目的

マーガリンの購買データから、ブランドごと、家計ごとのマーガリン価格に対しての反応の違いを明らかにしたい。

データ

bayesmパッケージにある、margarineデータ。data(margarine)で呼び出せ、詳細はcranのドキュメントに載っています。

Household Panel Data on Margarine Purchasesには、516家計の購買データと、家計ごとのデモグラフィック情報が収められています。1991年の論文のデータとなるので、かなり昔のデータです。

  • 購買データは価格(USドル)と選択したブランドのID(10種類)
  • デモグラフィック情報はfamily size(家族構成)、学歴、職位、退職の有無などのダミー変数

今回の事例では、5回以上購買した家計に限定して分析しているため、
313家計・3405の購買レコードからなるデータセットとなります。

モデル

  • 家計ごとに異なる、マーガリン価格に対する反応を想定。各マーガリンのブランドの価格に対するパラメータの数は家計の数だけある。
  • 価格に対する反応は家計の属性によっても決まる。
    という前提に立ち、以下のセッティングで推論していきます。

  • 6つのブランド選択に関する多項ロジスティックモデル(カテゴリカル分布とsoftmax関数の適用)

  • 1階層目はブランドごとの価格を説明変数とし、価格に対する反応係数をかけ合わせたものを多項ロジスティックモデルの入力とする。
  • 2階層目はブランドの価格に対する反応係数が家計ごとの定数項と属性データに属性ごとの係数をかけ合わせたものからなる。
  • 家計ごとの定数項は平均0、分散V_betaの正規分布に従う。
  • 属性ごとの係数は平均vec(delta_bar)、分散V_betaクロネッカーのデルタA^(-1)の正規分布に従う。
  • 分散V_betaは平均υ、分散Vの逆ウィシャート分布に従う。
  • A = 0.01、υ = 6 + 3 = 9、V = υI(Iは単位行列)

コード

kefitsさんがいくつかの章に登場するbayesmでの実践例をstanに書き直されているようですので、そちらのコードで学ばせていただこうと思います。
https://github.com/kefits/Bayesian-Statistics-and-Marketing

以下が、stanのコードとなっています。ここでは、Hierarchical_MNL.stanとして保存します。

以下はstanをキックするためのRコードです。

実行結果

Core i5、8GBメモリのMacBook Proで40分ほどかかりました。

traceplot(d.fit)で以下のように4回の試行結果が描かれますが、収束しているようです。

summary関数を使えばわかりますが、3913行ものパラメータたちのサマリーが得られます。

  • 313家計の家計ごとのブランドに対するパラメータ(1878個)
  • 313家計の家計ごとのブランドに対する潜在パラメータ(1878個)
  • 6ブランドの係数の共分散行列(36個)
  • 6ブランドの係数の分散のハイパーパラメータの行列(36個)
  • 6ブランドの属性データ(8つ)に対する係数(48個)
  • 6ブランドの属性データに対する係数の共分散行列(36個)
  • lp(log posterior(確率密度の和でモデル比較で扱う。))(1個)

64番目の家計の各ブランドの価格に対する係数の分布を確認すると、4番目・5番目のブランドの係数が他のブランドに比べて小さいことがわかります。

続いて、家計ごとの係数に関して集計し、係数ごとの相関係数を見てみると、各ブランドごとに正の相関、負の相関がありそうです。

最後に、家計ごとに集計した、ブランドに対する価格反応係数の事後分布を描きます。

多峰性などはなく、正規分布に従っているようです。他のブランドと比較して、5番目の係数が小さいようです。

というのは誤りで、一週間後に気づいたのですが、家計ごとのブランドごとの係数の事後分布の平均値をプロットするべきでした。
正しくはこちらです。

事前情報として正規分布を仮定していましたが、係数に関して正規分布に従っていません。
そのため、事前情報として対称性のあるような正規分布を扱うのは適切ではなさそうです。

おわりに

2005年の本とは言え、十分に使いみちのある本だと思いました。まだまだ扱いきれていないですが、引き続き勉強していきます。
この本にはケーススタディが5つほどあるのですが、それのstanコード化などをしていけばかなり力がつくような気がします。

マーケティングの部署で働くデータアナリストにとって、マーケティング×ベイズの話は非常にモチベーションの上がるところなので、こういう文献を今後も見つけていきたい。

参考文献

Bayesian Statistics and Marketing (Wiley Series in Probability and Statistics)
Bayesian Statistics and Marketingのサポートサイト
ベイズモデリングによるマーケティング分析
StanとRでベイズ統計モデリング (Wonderful R)
RStanのおさらいをしながら読む 岩波DS 1 Shinya Uryu
Stanのlp__とは何なのか うなどん
‘LP__’ IN STAN OUTPUT
Package ‘bayesm’

蒙古タンメン中本コーパスに対してのLDAの適用とトピック数の探索

モチベーション

前回の記事では、Webスクレイピングにより入手した、蒙古タンメン中本の口コミデータに関して、Word2Vecを適用した特徴量エンジニアリングの事例を紹介しました。
今回はせっかく興味深いデータがあるので、どのようなトピックがあるのかをLDAを適用したいと思います。加えて、これまで記事で扱ってきたLDAの事例では評価指標であるPerplexityやCoherenceを扱ってこなかったことから、トピック数がどれくらいであるべきなのか、考察も含めて行いたいと思います。以前扱った階層ディリクレ過程であれば、トピック数を事前に決める必要が無いのですが、今回は扱わないものとします。

環境

・MacBook Pro
・Python3.5
・R version 3.4.4

Gensimで行うLDA

今回もPythonのGensimライブラリを用いて行います。

  • パープレキシティ
    • テストデータに対して計算
    • 負の対数尤度で、低いほどよい。
      • パープレキシティが低いと、高い精度で予測できるよい確率モデルと見なされる。汎化能力を表す指標。
      • トピックの数をいくらでも増やせばパープレキシティは下がる傾向が出ている。
      • 教科書でのパープレキシティの事例に関しては、トピック数を増やせば低くなるという傾向が出ている。

以下のコードでパープレキシティを計算します。

実際に、中本コーパスで計算したトピック数に対してのパープレキシティは以下のように推移しました。

Ldaのモデル選択におけるperplexityの評価によると、
”複数のトピック数で比べて、Perplexityが最も低いものを選択する。」という手法は人間にとって有益なモデルを選択するのに全く役に立たない可能性がある。”と記されています。

『トピックモデルによる統計的潜在意味解析』には、”識別問題の特徴量として使う場合は識別問題の評価方法で決定すればよい”とあるので、目的によってはパープレキシティにこだわらなくても良いと思われます。

今回のケースだと、パープレキシティだけだと、決めかねてしまいますね。

  • コヒーレンス
    • トピックごとの単語間類似度の平均
    • トピック全体のコヒーレンスが高ければ、良い学習アルゴリズムとみなす。

以下のコードでコヒーレンスを計算します。

実際に推定してみたところ、トピック数が20を超えたあたりからコヒーレンスが下がる傾向があるので、
それ以上のトピック数は追い求めない方が良いのかもしれません。

Rでもやってみる

Rでトピック数を決める良い方法がないか調べてみたところ、ldatuningとかいうパッケージがあることがわかりました。複数の論文(Griffiths2004, CaoJuan2009, Arun2010,Deveaud2014)で扱われている手法を元に、適切なトピック数を探れるようです。このパッケージを紹介しているブログの事例では、90から140の範囲で最適なトピック数となることが示されています。詳しくはこちらを見てください。
Select number of topics for LDA model

以下のコードで実行しました。一部、驚異のアニヲタさんのコードを拝借しております。なお、ldaパッケージのlexicalize関数を用いることで、ldatuningに入力するデータを作成することができます。

これを見る限りは、60〜70個の辺りに落ち着くのでしょうか。

トピックの吐き出し

Rでの結果から、60個程度のトピックで推定し、各記事に割り当てが最大のトピックを付与して、トピック別の口コミ評価をみてみようと思います。

以下のコードではトピック別の口コミ評価のしやすさからtopicmodelsパッケージを用いた推定となっています。

口コミ評価の点数が上位のトピックはこんな感じです。

口コミ評価の点数が下位のトピックはこんな感じです。

中本は社会人2〜3年目で新規メディアの立ち上げのストレス解消で数回行きましたが、北極の赤さは異常だと思います。北極を食べたり、トッピングする余裕のある人、ましてや辛さを倍にするという時点で口コミ評価も高くなると考えるのは自然なのかもしれません。

参考情報

トピックモデル (機械学習プロフェッショナルシリーズ)
トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ)
models.ldamodel – Latent Dirichlet Allocation
Ldaのモデル選択におけるperplexityの評価
pythonでgensimを使ってトピックモデル(LDA)を行う
gensim0.8.6のチュートリアルをやってみた【コーパスとベクトル空間】
LDA 実装の比較
Jupyter notebookにMatplotlibでリアルタイムにチャートを書く
Inferring the number of topics for gensim’s LDA – perplexity, CM, AIC, and BIC
Select number of topics for LDA model
47の心得シリーズをトピックモデルで分類する。 – 驚異のアニヲタ社会復帰への道

R Advent Calendar 2017 rvestを用いてポケモンデータをスクレイピング&分析してみた

R Advent Calendar 2017の11日目を担当するMr_Sakaueです。
今回はrvestパッケージを用いて、友人がハマっているポケモンの情報を集めてみようと思います。
もっとも、業務でWebスクレイピングする際はPythonでBeautifulSoupやSeleniumを使うことがほとんどなのですが、たまにはRでやってみようと思います。

目次
・やりたいこと
・rvestについて
・データの取得と集計と可視化と分析
・まとめ
・参考情報

やりたいこと

今回はポケモンたちのデータを集めた上で、以下の内容を行いたいと思います。

  • ポケモンのサイトから種族値を取得
  • ポケモンの種族値を標準化して再度ランキング
  • ポケモンのレア度や経験値に関する情報を取得
  • レア度や経験値と相関しそうな種族値を探る

今回扱った全てのコードはこちらに載せております。
https://github.com/KamonohashiPerry/r_advent_calendar_2017/tree/master

※種族値はゲームにおける隠しパラメータとして設定されている、ポケモンの能力値とされている。

rvestについて

rvestはRでWebスクレイピングを簡単に行えるパッケージです。ここでの説明は不要に思われますが、今回はread_html()、html_nodes()、html_text()、html_attr()の4つ関数を用いました。

基本的に以下の3ステップでWebの情報を取得することができます。

  • STEP1
    read_html()でHTMLからソースコードを取得する。(Pythonでいう、requestとBeautifulSoup)
  • STEP2
    html_nodes()でソースコードから指定した要素を抽出する。(PythonでいうところのfindAll)
  • STEP3
    html_text()やhtml_attr()で抽出した要素からテキストやリンクを抽出する。(Pythonでいうところのget(‘href’)など)

データの取得と集計と可視化

検索エンジンで検索してだいたい1位のサイトがあったので、そちらのWebサイトに載っているポケモンの種族値の一覧をスクレイピング対象とさせていただきます。

  • ポケモンのサイトから種族値を取得

以上のコードを実行すれば、こんな感じでポケモンの種族値一覧を得る事ができます。

とりあえず、種族値合計(Total Tribal Value 以下、TTV)のランキングの上位を確認してみます。知らないんですが、メガミュウツーとかいうイカつそうなポケモンが上位にいるようです。昭和の世代には縁のなさそうなポケモンばかりですねぇ。

■TTVランキング

取得した種族値を項目別に集計したり、Boxプロットを描いてみます。どうやら、攻撃の平均が高く、ヒットポイントや素早さの平均は低いようです。

■種族値のサマリー

■種族値のBoxプロット

  • ポケモンの種族値を標準化して再度ランキング

さて、攻撃の平均が高かったり、ヒットポイントと素早さの平均が低かったりしたので、各々の項目を標準化した上で、再度ランキングを作ってみたいと思います。

■標準化した種族値のサマリー

平均0、分散1にできているようです。

■標準化した種族値のBoxプロット

他よりも低かったヒットポイントと、高かった攻撃がならされていることが確認できます。

■標準化前後でのTTVランキングのギャップが大きかったものをピックアップ

ラッキーが144位ほど出世しています。攻撃が低く、ヒットポイントの高いラッキーが標準化により優遇されるようになったと考える事ができます。ポケモン大会の上位ランカーである後輩社員もラッキーは手強いですと言っていたのでまんざらでもないのでしょう。

  • ポケモンのレア度や経験値に関する情報を取得

今回のサイトには、個別にポケモン別のページが用意されており、そちらから、ゲットしやすさや経験値に関する情報を抽出します。

以上のコードを実行すれば、やや時間がかかりますが、全ポケモンのゲットしやすさや経験値のデータを抽出する事ができます。それらの情報がゲットできたら、まずは可視化します。

■ゲットしやすさのヒストグラム

ゲットのしやすさは、小さいほど捕まえる難易度が高くなっています。難易度の高いポケモンである0が多過ぎるので、このデータは欠損値が0になっているのではないかと疑われます。

■経験値のヒストグラム

経験値は、レベル100になるまでに要する経験値をさしています。ほとんどが100万程度となっているようです。

■ゲットしやすさと標準化TTVの散布図

やはり、ゲットしやすさに関してはデータに不備があるようで、コラッタ(アローラの姿)のような雑魚ポケのゲットのしやすさが0だったり、伝説のポケモンであるネクロズマが255だったりします。ただ、上限と下限のデータを間引けば右下がりの傾向が見られそうです。

■経験値と標準化TTVの散布図

経験値が多く必要にも関わらず、TTVが低い集団があります。どうやらこの集団に属するのは、「キノガッサ」・「マクノシタ」・「イルミーゼ」・「ゴクリン」・「シザリガー」などで、一回しか進化しないポケモンのようです。これらのポケモンは育てにくく、TTVの低い、コスパの悪そうなポケモンと考えることができるのではないでしょうか。(技や特性によってはバリューあるかもしれませんが。)

  • レア度や経験値と相関しそうな種族値を探る

先ほどのレア度に関しては、データがおかしそうだったので、レア度0と255に関しては除外してみます。

■ゲットしやすさと標準化TTVの散布図

やはり除外する事で、理想的な右下がりの傾向を示す散布図が得られたと思います。
さて、各種族値がレア度にどれだけ相関しているのかを分析したいのですが、その前にレア度を表す二項変数を作成します。

■ゲットしやすさが50以下であれば1、それ以外を0にする変数を作成

続いて、各種族値を説明変数として、レア度を目的変数としたロジスティック回帰モデルの推定をrstanで実行させます。

■stanコード

■rstanでロジスティック回帰を行い、推定結果を可視化するコード

■MCMCのシミュレーション結果のトレースプロット

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

■ロジスティック回帰の推定結果

見にくいので、推定結果を松浦さんの「StanとRでベイズ統計モデリング」にあるコードを用いて可視化します。

■推定結果の可視化

どうやら、0を含まない係数について見てみると、b3(攻撃)、b5(特殊攻撃)、b6(特殊防御)が高いほど、レア度が増す傾向があるようです。珍しいポケモンは攻撃が強いという傾向があると言えるのではないでしょうか。

まとめ

  • rvestは簡単にスクレイピングできて便利。
  • ポケモンデータは色々整備されてそうで今後も分析したら面白そう。
  • 珍しいポケモンは「攻撃」、「特殊攻撃」、「特殊防御」が高い傾向がある。
  • 経験値が必要なのにTTVの低い、コスパの悪そうなポケモンたちがいる。

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

参考情報

データサイエンティストのための最新知識と実践 Rではじめよう! [モダン]なデータ分析
StanとRでベイズ統計モデリング (Wonderful R)
【R言語】rvestパッケージによるウェブスクレイピング その2
Receiving NAs when scraping links (href) with rvest

階層ベイズモデルの直帰率分析への適用 with rstan

松浦さんの『StanとRでベイズ統計モデリング』の8章の階層ベイズがすごくわかりやすいなぁと思いつつも、自分の持っているデータで試していなかったので、これを機に実践してみようと思います。
やや変数を追加しているくらいで大した変更点はありませんが、題材としては当ブログのアクセスログにおける直帰率に関するデータで、どのような要素が直帰率に影響を与えるのかを分析します。

目次
・モデル概要
・前処理
・推定
・結果(非階層モデルとの比較)
・参考文献

モデル概要

モデルは8章のロジスティック回帰の階層モデルに一部変数を追加していますが、ほぼそのままです。記事ごとのパラメータやリファラーごとのパラメータを想定しています。
Nは記事数でnはそのインデックス、Cはリファラーの数でcはそのインデックス、Iはログとして残っているセッションの数でiはそのインデックスとなっています。hatebuは記事のはてぶ数、stringlineは記事の行数、holidayは休日or祝日ダミー変数、daytimeは12:00~18:00なら1をとるダミー変数、revisitedは再訪問ユーザーなら1を取るダミー変数となっています。記事ごと・リファラーごとに直帰のしやすさが違う(パラメータが従う正規分布のパラメータがそれぞれ異なる)という仮定のもとに立ったモデルとなります。

$$x[i] = b_{1} + x_{記事}[記事ID[i]] \\ + x_{リファラー}[リファラーID[i]] + x_{セッション}[i] $$

$$q[i] = inverselogit(x[i]) $$

$$Y[i] \sim Bernoulli(q[i]) $$

$$x_{記事}[n] = b_{2}hatebu + b_{3}stringline[n] \\ + b_{記事間の差}[n] $$

$$b_{記事間の差}[n] \sim Normal(0, \sigma_{記事番号}) $$

$$x_{リファラー}[c] = b_{リファラー間の差}[c] $$

$$b_{リファラー間の差}[c] \sim Normal(0, \sigma_{リファラー番号}) $$

$$x_{セッション}[i] = b_{4}holiday[i]  + b_{5}divice[i] \\ + b_{6}daytime[i] + b_{7}revisited[i] $$

前処理

GAのAPIからデータを取得して1セッション1記事になるようにデータを作成しています。数ヶ月で25000件ほどデータがあったのですが、計算に時間がかかるので、データ数を2400件くらいにサンプリングしています。

推定

stanコードはこちらになります。

rstanを用いたstan実行用のRコードです。ヴァイオリンプロットで主要な係数の分布を見る処理も書かれています。

結果

係数を見る限りは、符号の向きが確かなのはb5(PCダミー)とb7(再訪問ユーザーダミー)なので、PCの方が直帰しにくく、再訪問ユーザーの方が直帰しにくいという傾向があると考えることができます。

教科書ではAUCを非階層モデルと比較していましたので、比較してみようと思います。
AUCの計算を行うためのコードもGithubに載っていましたのでそちらを使います。

80%が良いとされているAUCには程遠いですが、記事やリファラーごとの差を考慮しない非階層のものよりもAUCが高いと言えます。
ちなみに、教科書の例のAUCは80%ほどでした。
Webマーケのデータ分析においてロジットは汎用性が高いで、今回のコードを土台に色々と業務で試していこうと思います。

参考文献

StanとRでベイズ統計モデリング (Wonderful R)
ベイズ統計モデリング: R,JAGS, Stanによるチュートリアル 原著第2版

RのFactoRizationMachinesパッケージを用いたFMのページビューデータへの適用

Googleタグマネージャーで集めたアクセスログデータを用いて、前回と同様に記事のレコメンドにチャレンジしてみようと思います。FactoRizationMachinesパッケージという便利そうなパッケージの存在も知れたことから、今回は以前から気になっていたFactorization Machineを扱います。

【目次】
・Factorization Machine(FM)の概要
・パッケージ紹介とインストール
・サンプルデータの構造把握と前処理
・FMの実行
・結果
・参考文献

Factorization Machine(FM)の概要

  • 組み合わせ特徴量を扱う教師あり学習モデル。行列分解とSVMを合体させた手法。
  • スパースになりやすいデータの予測問題で扱う。
  • 1ユーザーのある商品に対しての評価を、1評価1行として表して、ユーザーとアイテムの交互作用の特徴ベクトルを扱う。
  • 相互作用項に関して、時間や文脈などを自由に入れられる。
  • 相互作用項を次元圧縮する際の要素数を事前に決める必要がある。
  • Matrix Factorizationよりも精度が良いとされている。特徴量エンジニアリングなどで使われているようです。(Click-Through Rate Prediction

パッケージ紹介とインストール

FactoRizationMachinesパッケージは線形SVMと2次のFMと高次のFMを実行することができ、引数で正則化項も加えることができます。現段階においては回帰のみで分類問題への適用は今後の開発となるようです。CRANから普通に
install.packages(‘FactoRizationMachines’)
でインストールします。libFMexeパッケージの場合は、libFMをインストールしてパスを指定しておく必要がありますが、このパッケージに関しては不要となります。

サンプルデータの構造把握と前処理

FactoRizationMachinesパッケージのサンプルコードにおいては、MovieLensのデータがサンプルデータとして載せられていました。ユーザーのID(整数)、映画のID(整数)、評価(整数、5段階)、日時(整数)からなるデータに対して、sparseMatrixに変換していました。

今回は、前回の投稿非負値行列因子分解(NMF)でブログ記事のレコメンドをしてみると同じデータを使って、アクセスログデータに適用しようと思います。FactoRizationMachinesの形式に合わせるために、このブログのアクセスログも、クッキーのIDを整数に、記事のIDを整数に、閲覧回数を5段階(5以上を5に変換)に、日時を整数に変更しています。

FMの実行

デフォルトの設定c(1, 10)では線形のウェイトが有効で、2次の項の要素数が10で正則化項なしのFMを実行することになります。引数に関する詳しい情報はPackage ‘FactoRizationMachines’に書かれています。今回はサンプルを参考に正則化項ありでモデルを実行します。まず、アクセスログデータに対して、ユーザーのIDからなる整数ベクトル、記事のIDからなる整数ベクトル、セッションのあった日時のデータからなる整数ベクトルを作成し、sparseMatrix関数を用いて元データを変形し、80%のデータをトレーニングに、20%のデータをテストに割り当てます。さらに、テストデータに関して、予測値との平均二乗誤差を計算します。

結果

各モデルについての平均二乗誤差を計算しています。
線形モデルや高次元モデルよりも、2次の項を持つFMが精度が高いようです。

こちらは、この中で性能の良かった2次の項を持つFMの予測結果とテストデータの結果をプロットしたものです。4点を超える値をあまり予測できていないようです。今回はサンプルを回しただけなので、本来であれば次元の数kや正則化のセッティングをいろいろいじったり、相互作用項を新しく追加するなどして精度を高めることが必要です。

結果の比較だけでは仕事で使えないので、実際に予測した結果を取り出したいと思います。
実際に運用するとなると、ページIDを所与として、ページビュー数を0とおいて(型をそろえるため。NULLだとエラーになった)、任意のタイミング(date)を想定して、モデルにデータを適用し、評価の高いものをサジェストするスタイルになるのではないでしょうか。

この結果だと、ユーザー98に記事1を見せることに対して4.02点が与えられています。

参考文献

Factorization Machinesを今更読みました
Factorization Machines
High-order factorization machines with R #tokyor 61
Factorization Machinesのおはなし。
libFMexeを動かすまで (R Wrapper for the libFM Executable参照記事)
一歩Matrix Factorization、二歩Factorization Machines、三歩Field-aware Factorization Machines…『分解、三段突き!!』
[論文] Factorization Machines (ICDM 2010) 読んだ 22:41
Factorization machines with r
Factorization Machinesについて調べてみた

References   [ + ]

1, 2, 3, 4, 5. predict(model,data.test)-target.test)^2

非負値行列因子分解(NMF)でブログ記事のレコメンドをしてみる

概要

2017年7月より当ブログにおいて、Google AnalyticsのカスタムディメンジョンにUser IDを取得して残せるようにしました。これにより、ユーザーの一人一人が閲覧した記事の分析を行うことが可能となっています。そこで今回は、記事の数も50近くなってきたことから、レコメンド手法の適用を行ってみたいと思います。Exampleデータでレコメンドをしても、ドメイン知識がないと、どれくらい良さそうな推薦なのかの判断が付きにくいと考えていたので、自分で運営しているサイトのデータであれば評価をつけやすいのではないかと思われます。
手法としてはNMF(non-Negative Matrix Factorization)を適用します。

【目次】
・前処理
・データ確認
・行列因子分解(Matrix Factorization)
・レコメンド結果の確認
・感想
・参考文献

前処理

ユーザー単位の閲覧した記事データを作るためにGAのAPIで取得したデータを整形し、集計します。このようなデータを手に入れるには自身のブログにおいて、GTM(Google Tag Maneger)を設置して、GTMのコンテナにカスタムJS変数などを格納する必要があります。詳しくはこちらのブログを参照ください。( 2016年の新定番!ユーザーエクスプローラーをもっと活用するための簡単な方法

このような形で行にユーザーが、列に記事が、各要素にアクセスの有無が入るようになっています。本来は評価データを用いますが、評価データがないので、アクセスしたかどうかの1-0データとして扱います。

データ確認

以下より、行列の要素のうち、値が入っている要素は全体の2.9%であることがわかります。参考文献は2.6%をスパースと言っていたことから、今回のデータもスパースであると考えることができます。

行列因子分解(Matrix Factorization)

行列因子分解はユーザーごとのアイテム購買履歴などのデータからなる行列の因子分解を行うもので、レコメンドにおいては因子分解の結果として、ゼロ以上の要素を持つ場合に推薦対象とするなどの活用がされています。非負値行列因子分解は要素が非負となるように制約を付けて行列を分解するもので、正の値を取ることから推薦に使いやすいとされています。

今回はNMFパッケージを用いて、ほぼデフォルトな設定で結果を眺めてみますが。本来であれば以下のような考慮が必要となります。最適化などをがっつり進めようとするとクロスバリデーションなどを行ったり、計算リソースが結構かかったりすると思われます。

  • アルゴリズムの選択
    • brunet,KL,lee,Frobenius,offset,nsNMF,ls-nmf,pe-nmf,siNMF,snmf/r,snmf/lがある。
    • デフォルトではbrunetとなっている。
  • シーディングメソッドの選択(初期値の決め方)
    • 固定値(none),ランダム(random),ica,nndsvdがある。
    • デフォルトではランダムとなっている。
  • 実行回数の設定
    • 100〜200回を選ぶことが多い。
  • 並列処理
    • 処理量が多いので最適化などを考えると並列処理も重要。
  • 因子分解ランクの推定
    • ランクの数が多いと残差は減るが、過学習になりやすい。クロスバリデーションを行う必要がある。

レコメンド結果の確認

各ユーザーに対して推薦されたアイテムのうち、比較的値の大きなものに関してコメントをしたいと思います。

ユーザー1
「XGBoostやパラメータチューニングの仕方に関する調査」を閲覧したユーザー1に対して、「XGBoostのパラメータチューニング実践 with Python」を推薦できている一方で、「OpenCV&Pythonで画像の類似度を計算させる〜イケメンの顔比較」など関係性の薄そうな記事も推薦してしまっています。

ユーザー2
「人工知能学会全国大会2017のWebマーケティングで参考になりそうな研究9選」を閲覧したユーザー2に対して、「Kaggleで使われた特徴量エンジニアリングとアルゴリズムまとめ」を推薦できている。

ユーザー3
「洋楽の歌詞データでDoc2vecを実行してみる」を閲覧したユーザー3に対して、「LDA(潜在的ディリクレ配分法)まとめ 手法の概要と試行まで」を推薦できている一方で、「ベイジアンネットワークをRのbnlearnパッケージで推定して予測してみる」と関係性の薄そうな記事も推薦している。本来、近いと思われる、「Word2Vecでクラシックの楽曲情報をコーパスとして類似度を出してみる」を推薦できていない。

感想

  • 結果に関して、納得のいくものといかないもので半々な印象でデフォルトの設定では満足できない。
  • NMFはクロスバリデーションなどにより過学習を避けた形での最適なランク数やアルゴリズムの選定、初期値の設定などが必要と思われるので、そこらへんの取り組みも今後行いたい。
  • FM(Factorization Machine)は回帰項を追加できるので、今後挑戦してみたい。
  • 1-0のデータではなく評価のデータが欲しい。

参考文献

岩波データサイエンス Vol.5
非負値行列因子分解(NMF)によるレコメンドのちょっとした例
An introduction to NMF package
非負ッ値行列・ファクトリゼーション ~ 半群じゃないから難しくないもん!

ベイジアンネットワークをRのbnlearnパッケージで推定して予測してみる

ベイジアンネットワークの知見が無かったので、調べた情報をまとめています。一応、載せているスクリプトでRを用いて予測するということができます。

【目次】
・ベイジアンネットワークとは
・ベイジアンネットワークの用途
・ベイジアンネットワークの推定のステップ
・Rでの実行例
・おまけ:Webサービスへの応用例
・参考文献

ベイジアンネットワークとは

・複数の確率変数の間の定性的な依存関係をグラフ構造によって表し、個々の変数の間の定量的な関係を条件付確率で表した確率モデル。前提として、有向非循環(グルグルと回らないグラフ)となっているグラフの構造を持つものに限定している。
・入力となる変数と出力となる変数はモデルの中では区別されない。
・時間という明確な因果関係などをモデルに組み込みやすいので、系列データなどを扱うケースが多い。モデル設計者がデータが生成されるプロセスを考慮しながらモデルを組んでいける。
・非循環性とd分離の仮定のみによって導かれる現在考えられる最も自然な離散モデルであり、現在の様々なモデルの中でも、最も表現力と予測力を持つモデルとされる。

ベイジアンネットワークの用途

・物体追跡
・ジェスチャ認識
・Webサイトなどのレコメンデーションサービス
・広告配信
・メルマガなどの配信最適化
など画像処理問題から消費者行動問題までいろいろな分野で活用されているようです。

ベイジアンネットワークの推定のステップ

モデルの選択や推定は学習 (learning)と呼ばれ、以下の2段階のステップを踏みます。

Step1:構造学習(structure learning)…データからネットワーク構造を学習する。
Step2:パラメータ学習(parameter learning)…Step1で学習した構造によって意味付けした分布のパラメータを学習する。

ベイジアンネットワークの推定は条件付き確率やMAP(Maximum a posteriori)推定量などを用います。構造学習の際に、事前知識を導入することができます。

Rでの実行例

・データ
何千もの主要な免疫系細胞から取り出した11種類のリン酸タンパク質とリン脂質の測定値からなるデータです。以下のサイトよりダウンロードできます。
Supporting Online Material Causal Protein-Signaling Networks Derived from Multiparameter Single-Cell Data

・パッケージ
ベイジアンネットワークを実行できるパッケージはいろいろありますが、参考文献に従いbnlearnを用います。加えて、ネットワーク図の可視化のためにRgraphvizを用います。Rgraphvizがcranでは対応したものがなかったので、こちらを参考にして、インストールしました。
Provides plotting capabilities for R graph objects

・ネットワークの図示

・データ確認

・構造学習
データセットから、ネットワーク構造を推定します。推定するための手法は様々あるようですが、今回はヒルクライムアルゴリズムを用います。

構造学習より得られたネットワーク構造を可視化します。因果関係がダメな感じのネットワークが出来ていないかチェックする必要がありますね。機械的に作るだけでなく、背景的知識も考慮するという進め方が推奨されています。

・パラメータ学習
先ほど、データセットから求めた構造について、以下のコードでパラメータ推定します。パラメータの推定や構造推定は試行錯誤するところなので、こんな簡単に済む訳では無いようです。

・予測
テストデータを用意して、データセットから求めた構造・パラメータを用いて、任意の変数の予測を行います。予測したい変数の値は欠損していたらエラーで予測してくれないですが、全部0にしておけば回ります。運用上は予測したい変数の値はわからないはずなので。

全然45度線上に乗っていないので、あまり精度は高くないようです。

cpquery関数と言って、条件付き確率を予測できる関数もあるようです。
Perform conditional probability queries

おまけ:Webサービスへの応用例

参考文献には、Webサイトの閲覧データをもとに、ユーザーの行動を予測するというモデルの紹介がなされていました。ただ、Webサイトの場合、ページ数や商品の数・ユーザーの数も膨大なことから、ネットワーク構造の推定が難しいようです。そこで、ネットワークを作るに際して、変数を確率的潜在意味解析(pLSA)で扱いやすい数に絞るなどの工夫をされていました。自社のサイトに関しても適用する際に、変数の多さは待った無しだと思うので、いざ適用する際は情報圧縮技術を駆使したいですね。

参考文献

ベイジアンネットワーク技術 ユーザ・顧客のモデル化と不確実性推論

確率的グラフィカルモデル

Learning Bayesian Networks in R an Example in Systems Biology

グラフィカルモデル入門

ベイジアンネットワークを用いたWeb レコメンデーションシステムの開発

PRML8章

Billboard100位以内の楽曲の歌詞情報にLDAを適用してみた

目次

・はじめに
・データ収集
・Rによる分析
・LDAの結果
・参考文献

はじめに

前回の投稿でBillboardの週次洋楽ランキングデータをWebスクレイピングで取得し、楽曲の消費サイクルのような順位の挙動を確かめることができました。(某洋楽ヒットチャートの週次ランキングデータをBeautiful Soupで集めてみた)今回は、歌詞の情報を用いて順位データとつなぐことにより、どのような単語の入っている洋楽がBillboardにおいてTop10に入る傾向があるのかをLDAを行うことで確かめたいと思います。

データ収集

残念なことに、Billboardのサイトに歌詞の情報は載っていません。そこで、洋楽の歌詞が取り上げられている某サイトをPython(3系)でWebスクレイピングし、名寄せを頑張って順位データと歌詞データを繋ぎます。

幸いなことに某サイトのURLに規則性があったので、アーティスト名からなるURLを生成し、そのURLをWebスクレイピングして楽曲のリストを集め、今回のBillboardのランキングに入った楽曲のみに絞ります。

楽曲をランキングに含まれるもののみに絞ったら、歌詞詳細ページを取得します。

うまいこと歌詞情報を手に入れることができました。ざっと947曲です。

Rによる分析

ここから、Rにてテキストマイニングを行いたいと思います。まず、tmパッケージを用いて、不要語(stop word)を除去します。具体的にはtheとかyouとかを除外しています。

続いて、LDAを実行できるtopicmodelsパッケージで扱えるようにするために、テキストデータに以下の処理を施します。

あとは以下のコードでLDAを実行するだけです。トピック数はアドホックに20としています。研究者の方、いい加減ですみません。

LDAの結果

まずは推定されたトピックごとの上位10単語をみてみます。トピック1はラブソングとかでしょうか。トピック17にパリピっぽい単語が、トピック18にスラングが含まれていますね。

見ずらいので、行を一つにまとめて、トピックにidを割り振ります。

最後に、BillboardでTop10に入ったかどうかのデータを作っておき、そのデータと各歌詞を繋ぎ、各歌詞ごとに割りふられた確率が最大のトピックで集計をします。

BillboardのTop10ランクイン割合の高いトピックTop3
「one,ooh,call,cause,gettin,born,day,makes,came,stand」
「better,world,whoa,run,light,things,find,show,see,waiting」・・・明るい感じ?
「stop,just,hands,put,party,crazy,live,lights,play,see」・・・パリピぽい

BillboardのTop10ランクイン割合の低いトピックTop3
「wanna,want,take,rock,see,kiss,come,make,body,tonight」・・・欲求系?
「feel,heart,life,away,just,break,real,enough,every,find」・・・癒し系?
「hey,said,old,every,woo,left,told,nothing,daddy,sweet」

あまり洋楽を聴かないので、得られたトピックの解釈が中々できないのがもどかしいです。ただ、スラングの歌詞を含む歌詞はそんなにランクイン割合が悪いわけではなさそうですね。洋楽をもっと聴いて、前処理などもう少し工夫してリベンジしたいですね。

参考文献

トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ)

Pythonクローリング&スクレイピング -データ収集・解析のための実践開発ガイド-

モダンなRによるテキスト解析topicmodels: An R Package for Fitting Topic Models

某洋楽ヒットチャートの週次ランキングデータをBeautiful Soupで集めてみた

はじめに

知人より、洋楽の流行りに疎いのでキャッチアップしたいという要望があり、某洋楽ヒットチャートの週次ランキングとTop100のデータを大量に集めてみようと思うに至りました。今回は深い考察を行うには至っていませんが、簡単にRにて集計・可視化を行います。

データ収集

Webスクレイピング対象の某洋楽ヒットチャートの週次ランキングは今週の順位・先週の順位・アーティスト名・曲名・詳細ページへのリンクなどが載せられおり、毎週土曜日更新されています。サイト内から導線はありませんが、URLのパラメータに法則があるため、うまく収集できます。今回は2010年8月〜2017年6月の約7年分のデータを集めます。

URLのリストをCSVで読み込み、BeautifulSoupでタグを指定して抽出します。