2018年に参加したデータ分析系の勉強会で得た知識の詰め合わせ

社内に分析チームがないことから、私は月に3~4件は刺激を求めて勉強会に足を運んでいます。新しい知見を得れることは然ることながら、社内だともらえないフィードバックをいただけたり、課題の共有などをできるのが良いと思います。

目の肥えた皆さんにとって新規性のある情報はあまりないかもしれませんが、詰め合わせた情報をお楽しみください。

統計学まわり

  • 勉強会名
    KDD論文読み会

論文読んだ「Winner’s Curse: Bias Estimation for Total Effects of Features in Online Controlled Experiments 」

機械学習まわり

  • 勉強会名
    merpay×M3 機械学習 NIGHT

    • 会社名
      M3
    • 知見
      コンテンツをレコメンドする際のテクニックとして、MFとCNNの合わせ技について紹介されていました。訓練時には、アクセスログデータをもとにMFで潜在的な表現を抽出しそれのアイテム間の類似度を計算し、推薦時には、テキストのタイトルとキーワードなどをCNNで学習し訓練時と同じ次元になるようにアイテムのベクトルを出力する。そして、訓練時のものと近いアイテムを推薦することでCold-Start問題を克服するとのことでした。
    • 発表資料
      Matrix Factorization と Text CNN による Cold Start Problem への取り組み

The Road to Machine Learning Engineer from Data Scientist

  • 勉強会名
    NetaDashi Meetup

    • 会社名
      NRI
    • 知見
      Elmoを用いた文書分類。Word2Vecなどではできなかった、文脈を考慮して類似度などを算出できる。
      Elmoの多言語対応に関しては、このGitHubを参照すると良いらしい。
      https://github.com/HIT-SCIR/ELMoForManyLangs

異常検知の評価指標って何を使えばいいの? / Metrics for one-class classification

  • 勉強会名
    グリー開発本部 Meetup #1 DataEngConf NYC報告会

    • 会社名
      GREE
    • 知見
      Contextual Banditについての紹介
      Artwork Personalization at Netflix」という記事で2017年ごろに取り上げられていたようです。

エンジニアリングまわり

  • 勉強会名
    グリー開発本部 Meetup #1 DataEngConf NYC報告会

    • 会社名
      GREE
    • 知見
      LUIJI
      ・メリットとしては、少なくとも以下のものがあるそうな。
       ・Pythonで書ける
       ・エラーの途中で処理を止めて、それを解消したら、止めたポイントから開始できる
       ・様々なツール群と連携できる柔軟性
       ・10行程度でスクリプト書ける。
       ・複雑な依存関係も描ける。
    • 発表資料
      https://www.slideshare.net/greetech/dataengconf-nyc18-1

  • 勉強会名
    bq_sushi tokyo #9 2018総集編

    • 会社名
      オープンハウス
    • 知見
      BigQueryGIS
      BigQueryからGISの情報を扱うことが可能になったらしい。顧客の希望する物件の情報をレコメンドするために地理情報を扱うらしいです。
      ただ、基準とする測地系が国によって異なり、それらを考慮しないで推薦すると1~2kmはズレてしまうとのこと。家買う際にそんだけズレるとキツイですね。こちら(BigQueryGIS: Google und PostGIS )はBigQueryGISに関連した情報を漁って見つけた記事ですが、BQで抽出した情報をそのままGoogleMapに表示できるのは面白いですね。

データ分析のツラミ系

  • 勉強会名
    merpay×M3 機械学習 NIGHT

    • 会社名
      M3
    • 知見
      メタデータの検索システムについて
      データセット名、テーブル名、カラム名、カラムのディスクリプションをキーワードで検索できる。
      日次でディスクリプションを取ってくるようにしている。どのドキュメントが一番見られているのかもモニタリングできるとのこと。似たような取り組みとして、リクルートがMetaLookingとかいう内製ツールを作っていたりしますね。私は各サービスごとのDBのテーブルの注意点などを適宜スプレッドシートに残す程度しかしていませんが、分析者がすぐにキャッチアップできる環境は重要ですね。

  • 勉強会名
    MLCT

    • 会社名
      ???
    • 知見
      事業計画書を作るリーンキャンバスの機械学習版とも言える、機械学習キャンバス0.1というものが質疑応答の際に紹介されていました。

  • 勉強会名
    グリー開発本部 Meetup #1 DataEngConf NYC報告会

    • 会社名
      GREE
    • 知見
      データリーク問題はどこも苦しんでいる?
      SalesForce社が顧客企業15万社の情報を活用して、機械学習モデルを構築しようとしたが、
      蓄積されたデータにおいては、ビジネスプロセスをやたらと予測できてしまうようなデータリーク問題が起きまくっていた。
      原因としては、データサイエンティスト不足(分析を前提としたデータ蓄積ができていない。)、手入力によるラベリングミスなどがあるらしい。
      どこの企業も苦しんでいると思うと、分析を前提にスナップショットを残し続けるという取り組みは競争優位性につながるのだろうか。
      SalesForce社は、訓練と検証の精度の差が大きいと注意したり時系列データを確認するなどして、データの信憑性に気をつけてモデルを作ったそうです。
      15万社にうまくフィットするモデルなので、精度は70~75%で満足できるものらしい。
    • 発表資料
      https://www.slideshare.net/greetech/dataengconf-nyc18-1

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時であがったけど、その後の準備で疲れ果てた。目にクマが若干できた。

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

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

学習済み分散表現を用いた文書分類に挑戦(一部再学習も)

はじめに

2018年9月のテキストアナリティクスシンポジウムに行った際に、学習済みの分散表現で事前学習したモデルを使って分類してうまくいく事例が紹介されていました。
全てのタスクにおいてうまくいくとは思えませんが、試すコストはあまりかからないので試してみます。

2017年のテキストアナリティクスシンポジウムにおいても、メルカリやGunosyでは分散表現を用いた手法が一番精度が高いと言われていましたし、今年の会ではNLP系の学会でも分散表現はデファクトスタンダードになっているという話も伺いました。
2013~14年はLDAを使った研究が多かった気がしますが、徐々にシフトしていっているんですね。

これまで(Word2Vecを用いて蒙古タンメン中本の口コミ評価を予測してみる)は4000件程度の蒙古タンメン中本の口コミの情報を元に分散表現を手に入れていましたが、学習済みの分散表現を用いたアプローチも有効かもしれないと思い、試してみようと思います。

分類タスク

某グルメ口コミサイトの蒙古タンメン中本の口コミのテキストから、3.5点以上の評価かどうかを予測するタスクを扱います。
本当は、ポケモン図鑑の説明文から水やら炎やらのタイプを予測するとかをしたいのですが、あいにく手元にデータがないので、以前集めた蒙古タンメン中本の口コミを使います。(実は後日、ポケモン図鑑のデータを集めたのですが、平仮名にまみれたデータな上に、データ数も800件しかなかったので、どのみち厳しかったです。)

学習済み分散表現

Word2Vecなどで大量の文書をもとに学習させた分散表現のことを指します。
大規模コーパスで分散表現を手に入れる際は、数十GBにも相当するテキストデータを数時間かけて推定するので、学習済みのモデルは非常にありがたいです。(4年前に会社のPCで計算した際は、12時間くらいかかったこともありました。)

無料で提供してくださっている分散表現については、すでにこちらのブログで紹介されています。そこで紹介されているものに少し付け足すと、日本語の分散表現に関しては以下のようなものがあります。

  • 白ヤギコーポレーションのモデル:Gensim
  • 東北大学 乾・岡崎研究室のモデル:Gensim
  • Facebookの学習済みFastTextモデル:Gensim
  • NWJC から取得した単語の分散表現データ (nwjc2vec):Gensim
  • NNLM embedding trained on Google News:TensorFlow

そこで、今回は各種学習済み分散表現と蒙古タンメン中本コーパスで求めた分散表現の文書分類の性能バトルをしてみたいと思います。
ただ、分散表現ではなく、単語の頻度をもとに特徴量を作ったものが一番精度が高いのですが、分散表現同士の比較でもってどの学習済み分散表現が中本の口コミ分類に役に立ちそうなのかを明らかにしようと思います。(本来は分析という観点から即でボツですが、見苦しくも比較していきます。)

前処理

前処理は以下の通りで、テキストデータを分かち書きして、数値や低頻度・高頻度語を除外しています。

処理を施すとこのようなデータになります。

特徴量は、scikit-learnのCountVectorizerやTfidfVectorizer、分散表現の合計・平均・TF-IDFを求めたものを用意します。

蒙古タンメン中本の口コミ4000件から作成した分散表現:Gensim

まず、以前のブログで紹介した蒙古タンメン中本の分散表現ですが、以下のように推定しています。

Pipelineを用いてExtraTreesClassifierによる学習をします。特徴量は先程あげた、テキストベースのCountVectorizerやTfidfVectorizer、分散表現の合計・平均・TF-IDFで、評価指標はAUCのクロスバリデーションスコアとします。

汗に関してコンテキストの似ている単語を抽出しています。

結果は、以下の通りで、分散表現を使わない方がAUCが高いです。ただ、w2v_tfidf(分散表現のTF-IDFを特徴量にしたもの)が分散表現の中でAUCが高いようです。今回はこの60.5%をベースラインに比較していこうと思います。

白ヤギコーポレーションのモデル:Gensim

こちらのリンク、「word2vecの学習済み日本語モデルを公開します」から、ダウンロードしてそのまま以下のコードでモデルを扱えます。

汗の関連語を抽出していますが、中国の歴史の何かですか?可汗とかいう単語は聞いたことあるかも。

まずは白ヤギさんの分散表現をそのまま使って予測してみます。(コードは先程のものとほぼ重複するので省略しています。)
残念ながら、ベースラインの60.5%には至りませんでした。

hogehoge.modelというフルモデル形式の場合は、再学習が可能です。詳しくはこちら(models.word2vec – Word2vec embeddings model)に書かれています。

今回は、白ヤギさんの分散表現に対して、追加で蒙古タンメン中本のテキストを食わせて再学習させます。

ベースラインの60.5%よりも下回り、さきほどの白ヤギさんのもともとの分散表現よりも下回りました。

再学習してもかえって精度が下がったりすることから、簡単に精度が出るわけではなさそうです。まぁ、理想はその適用領域での大量のテキストデータがあることで、Wikipediaを元に作成した分散表現に強く依存しても駄目なのだろうと思われます。

東北大学 乾・岡崎研究室のモデル:Gensim

日本語 Wikipedia エンティティベクトルからダウンロードした学習済み分散表現を用います。ダウンロード後は普通に「gzip -d file.txt.gz」みたいにターミナル上で解凍します。以下のコードを実行すればすぐに使うことができます。
ただし、KeyedVectors形式のものは白ヤギさんのように再学習ができません。(Why use KeyedVectors instead of a full model?

汗の類似語に関しては、難しい単語が高めに出ているようです。

残念ながら、ベースラインの60.5%には至りませんでした。

Facebookの学習済みFastTextモデル:Gensim

FastTextはGoogleにいたTomas Mikolov氏がFacebookに転職されて作られた分散表現を求めるためのモデルです。Gensimでも呼び出せます。学習済みのものはこちらのGitHub(Pre-trained word vectors)にあるのですが、NEologdで形態素解析したものをベースに学習し公開されている方がいるとのことで、こちら(fastTextの学習済みモデルを公開しました)からダウンロードしたものを使わせていただきました。

何だこれはレベルの結果が返ってきました。中国の歴史上の人物か何かなんでしょうか。

若干ですがベースラインの60.5%よりも良い結果が得られましたが、 誤差の範囲な気がします。

NWJC から取得した単語の分散表現データ (nwjc2vec):Gensim

国立国語研究所の収集されたテキストデータを元に学習した分散表現が提供されています。ただし、利用するためには申請する必要があります。申請が受理されたらこちら(NWJC から取得した単語の分散表現データ (nwjc2vec) を頒布)からダウンロードして使えます。

汗の関連語ですが、うまく関連付けれているように思われます。少なくとも中国史ぽくはありません。しかしながら、顔文字まで学習していたとは。

ベースラインの60.5%よりも1%ポイントほど高い結果となりました。

NNLM embedding trained on Google News:TensorFlow

こちら(tensorflow-hubで超簡単にテキスト分類モデルが作成できる)で紹介されているように、GoogleがTensorFlowでGoogleニュースのテキストをもとに学習した分散表現が提供されています。

こちらのGitHub(NNLM embedding trained on Google News)から、Japaneseのnnlm-ja-dim50、nnlm-ja-dim50-with-normalizationなどが使えます。分散表現の説明についてはこちらのドキュメント(Token based text embedding trained on Japanese Google News 6B corpus.)にあります。

AUCが65%となっているものの、先程のsklearnでのクロスバリデーションのものとの比較ではないので、なんとも言えないですが、Googleニュースのデータだし結構精度が出そうな可能性を感じますね。
今後、TensorFlowでクロスバリデーションによるAUCスコアの出し方を調べてみて、順当に比較できるようにしたいです。(Kerasを使って計算している事例は見つけた。)

比較

今回の分類タスクはそもそも分散表現では精度が出なかったのですが、学習済み分散表現の中で序列を作るとすると、梵天が一番良く、FastTextが少しだけ良かったです。
TensorFlowをほぼ業務で使わないので、Googleニュースの分散表現を今回の比較対象にできなかったのですが、後日比較できるようにしたいと思います。

あと、今回の口コミの点数を当てるタスクよりも、分散表現にとって相性がいいタスクがあるかもしれないので、今回の結果で諦めることなく色々と試して行きたいです。

おわりに

様々なシンポジウムなどでスタンダードとなってきた分散表現ですが、学習済み分散表現をそのまま使って分類問題で役に立つのかを見てきました。残念ながら、口コミの評価予測タスクにおいては全然効果がなさそうでした。ただ、分散表現の中でもタスクによって相性の良い学習済み分散表現がありそうです。
先程も述べたように、理想は大量のテキストデータで学習した分散表現を求め、それを予測に使うことなので大量のテキストデータを集めて再チャレンジしたいです。どれくらいのテキストデータがあれば十分なのかの規模感もわからないので、実践あるのみなんですかね。

参考情報

Word Embeddingだけで文書分類する
tensorflow-hubで超簡単にテキスト分類モデルが作成できる
Error: ”Word2vec’ object has no attribute index2word
Word2vec Tutorial Online training / Resuming training
Word Embeddingモデル再訪
Googleの事前学習済みモデルを手軽に利用出来るTensorFlow Hub
ゼロから作るDeep Learning ❷ ―自然言語処理編

pyreaperで音声データのピッチを掴むためのメモ書き(F0の抽出)

はじめに

最近、仕事で音声データを扱うことがあるのですが、アプローチの仕方を色々と調べています。

今のところ検討しているのが、

  • クラウドAPIなどによる文字起こしによるテキストアナリティクス
  • 音声の波形に関する分析
  • 音声のピッチに関する分析

の3つですが、1つ目と2つ目については比較的すぐにググってチャレンジできましたし、ビジネス上で役に立ちそうな発見もありました。

そこで今回は3つ目の「音声のピッチ」(音程)について扱おうと思います。

ピッチがわかると何が嬉しいかというと、顧客対応などの音声でピッチがどのようなものであれば顧客満足度が高いのかを分析できるようになり、担当者へのフィードバックに使える可能性があります。
今までなんとなく蓄積していた音声データが宝の山になるかどうか、それをジャッジするための一つの手段と考えることができます。

後輩からの聞き込み

大学院で動画の圧縮について研究していた元後輩社員くんがいたので、きっと知っているはずだろうと思い、その方からピッチに関する研究を進める上でのキーワードを教えてもらいました。
彼曰く、F1やF0というキーワードで探していけばピッチの抽出ができるそうなので、早速調べて実践できないかチャレンジします。

F0とF1

Wikipediaで調べて結果を以下に記します。

  • F0:基本周波数
    信号を正弦波の合成で表したときの最も低い周波数成分の周波数を基本周波数と呼ぶ。

  • F1:フォルマント周波数
    フォルマント:言葉を発している人の音声スペクトルを観察すると分かる、時間的に移動している複数のピークのこと。
    周波数の低い順に、第1フォルマント、第2フォルマントというように数字を当てて呼び、それぞれF1、F2とも表記する。
    観察の仕方としては、縦軸に周波数、横軸に時間を置くものとする。

ピッチはどうすれば求めることができるのか?

断創録さんのブログによると、ピッチとは声の高さのことで、基本周波数やF0とも呼ばれているそうです。
そのため、今回はF0を求めれば目的を果たすことができそうです。

今回扱うモジュール/今回のツール

今回は、ピッチ検出を行うためのツール、REAPERのPythonラッパーであるpyreaperを使います。

一発で入るはずです。

続いて、音声ファイルをWAVファイルに変換するためのツールとして、FFmpegをインストールします。
Macだと、

で一発で入るはずです。
使い方としては、

のようにターミナル上で使います。

データ

今回は、無料の音声素材でなにか無いかと探していたのですが、NHKの百人一首の読み上げデータにしました。
句に関してはこだわりはないので、適当に第5句にしています。

NHKの百人一首 第5句
奥山に紅葉踏み分け鳴く鹿の 声聞く時ぞ秋は悲しき/猿丸大夫

もう秋ですし、風情があって良いですね。

プログラム

実行プログラムはご丁寧にpyreaperのPyPIのページに良いものがありましたので、こちらをまるごと使います。

まず、百人一首の音声がm4a形式だったので、wav形式になるように、
先程インストールしたFFmpegを用いて、サンプルレート16000、チャネル数1のWAVファイルに変換します。

なお、電話などの話者が2名いるステレオのデータの場合、FFmpegを用いて左側成分、右側成分だけを抽出することができますので、分けて解析すると良いと思います。左右の分割に関してはFFmpegの参考文献で紹介されています。そうすれば、片方の話者の声だけを抽出することができますので、解析が捗ります。

実行結果

百人一首の音声の波形です。

ピッチマークです。ピッチマークとは、「単位波形を抽出する基準点」とされているとの説明があるのですが、あまり良くわかりません。無音のときに生成されていないので、ピッチ抽出の確信度という意味なんでしょうか。

F0です。音程とみなすことができるので、こちらを会話の分析の際に使ってみるのが良いのだと思います。

相関係数です。こちらもピッチ検出に使う指標のようです。

今後の使いみち

F0を集計してある水準のものが多い、時系列で見ると前半と後半でF0の水準が違っているというような情報と、その会話におけるラベル(会話の結果としての商談の成否)について分析することで、商談などがうまくいく際のヒントが得られるかもしれないですね。

追記

数十分の音声に関してF0を抽出しようとしたら、処理に非常に時間がかかり、結局最初の数分とかのF0の抽出にとどまりました。予感はしていましたが、結構計算リソースが必要な領域なのですね。

参考情報

IPythonデータサイエンスクックブック ―対話型コンピューティングと可視化のためのレシピ集
音声生成の基礎と音声学
Is there a way to convert audio files in Mac OS X or the command line without using iTunes?
やる夫で学ぶディジタル信号処理 東北大学 大学院情報科学研究科 鏡 慎吾
[vDSP][信号処理]オーディオ・音声分析への道7 FFT 自己相関関数 ピッチ検出
REAPER: Robust Epoch And Pitch EstimatoR
FFmpeg wiki:AudioChannelManipulation
基本周波数抽出器 REAPERによる抽出実験
SPTKの使い方 (3) ピッチ抽出

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分布を、尤度にロジスティック分布を設定しています。

ほとんど係数の大きさが同じであることが確認できます。

ベイズ推定の良いところは事後分布から関心のある係数に関しての取りうる値などをシミュレーションできるところですが、
posterior_interval関数で簡単に計算することができます。

肝心のMCMCの収束診断ですが、shinystanを使います。

やや余談ですが、他のデータセットでshinystanを用いた際に、予測結果にNAsが含まれている場合に、
shinystanが起動しないという問題があり、以下のようなエラー文が吐かれます。

調べたところ、こちらのgithubにあるように、

のように引数でppd=FALSEのように設定することで、立ち上げることができました。

3つの基準をクリアしているため、収束しています。

係数の分布についても可視化します。

rstanarmの良い点の一つとして、モデルのアップデートが容易に行える点があげられると思いますが、実際、以下のように先程のモデルに変数を追加して推定することができます。
今回は、alcoholを二乗したものを新しい変数として加えます。

次に、looパッケージを用いて、更新したモデルと元のモデルの性能の比較を行います。
looパッケージは統計モデルの予測精度の指標として扱われる、WAIC(Widely Applicable Information Criterion)を計算するためのパッケージで、WICは事後分布から得られる対数尤度の平均や分散からなる値として表されます。looはleave-one-out cross-validationのleave-one-outの頭文字。

さっそく入れようと思ったところ、

というエラーが出ました。
こちらでも議論されていましたが、

でバージョンを2.0.0から1.1.0に落としたら動きました。

ここで、事後分布が特定のサンプルデータに対して敏感であるかどうかをlooパッケージを用いて可視化します。

縦軸のshape parameter kは推定の信頼性の指標とされ、大きければ大きいほど信頼できないと見なし、横軸は今回推定したワインデータのデータの番号で、左が元のモデル、右が変数を追加したモデルのものです。
どうやらどちらも0.4未満のkに収まっているようです。参考情報の事例では0.5を超えていましたが、moderate outliersと説明されていたので、今回の推定は問題ないと思われます。

続いてモデルの比較を行います。

elpd_diffに関しては右のモデルの精度が高ければ正の値を、低ければ負の値を取るようになっています。標準誤差も返されます。
どうやら変数を追加したモデルの方が、ちょっとだけ良さそうです。

続いて、事後予測分布から、どの変数がどのように予測に影響を与えるのかを確かめます。
比較のためにデータを2つほど作成し、両者において一つだけ変数が違うという状況下での、予測される確率の比較を行います。

他の要素をある一定水準で保った際に、alcoholだけ1度下げることで、平均3%ほど高い評価が得られる確率が高まるという考察となります。

以上で、rstanarmの一連の使い方となるのですが、
一部の関数においては、階層ベイズモデルも行えるので、試してみようと思います。

ただ、階層ベイズにするにも、赤ワインのデータしかないので、グループ変数をどうにか作らないといけません。
あまりやりたくはありませんが、データがないので、説明変数を元にk-means(K=3)によるクラスタリングを行い、それをグループ変数とします。

stan_glmer関数を使えば、以下のような簡単な記述で定数項や係数がグループごとに異なるパラメータの分布に従うとする階層ベイズモデルを推定できます。

今回は、切片だけがグループごとに異なるモデル、傾きだけがグループごとに異なるモデル、切片も傾きも異なるモデルを作成します。

先程紹介した、looパッケージを使って、ベースとなるモデルとの比較を行います。

うーん、残念ながらどのモデルもベースモデルよりも圧倒的に強いものは無さそうです。

感想

まだまだrstanarmの関数やら機能やら定義を全て把握しきれていないですが、そこらへんがクリアーになれば、これまでのstanでの推定業務において生産性が高まる可能性を感じました。
簡単な階層ベイズモデルくらいなら、非常に直感的に書ける点や、変数の追加によるモデルのアップデートが容易な点などはポイント高いなぁと思います。
とはいえ、実務としてマニュアルでstanコードを作成していくのは必須なので、このパッケージを使うことによって、stanコードの改善に時間をより一層割けるようになるなら、それが一番だと思いました。
あと、「ベイズ初めてです!」という新入りの方とかには慣れ親しんでもらうには良さそうですね。lm関数レベルで実行できてしまうので。
今回、mc-stan.orgの配下にあるページなどを漁る過程で、ベイズ推定結果の可視化などで知らないことにも色々と出会えたので、今後も読み進めていきます。

追記

2018-09-10: Stanを使って変数選択したいにprojpredというパッケージが紹介されており、これを使えば、情報量基準に従った変数選択を簡単に行えるそうです。こうなると、「ベイズ推定に慣れ親しんでいない頻度主義系のソフトウェアユーザー」に限らず多くの人が幸せになれるパッケージなのかもしれませんね。

参考情報

Using the loo package (version >= 2.0.0)
Leave-one-outクロスバリデーションの2つのデメリット、からの解決方法
stan_glmer | Bayesian Generalized Linear Models With Group-Specific Terms Via Stan
WAICを計算してみる
Package ‘bayesm’
Hierarchical Partial Pooling for Repeated Binary Trials
Leave-one-out cross-validation for non-factorizable models
priors | Prior Distributions And Options
StanとRでベイズ統計モデリング (Wonderful R)

Stanで順序プロビット(Ordered Probit)の推定のためのメモ書き

最近はBayesian Statistics and Marketingという本に関心があって、そこで取り上げられているモデルをStanに落とし込めないか模索しています。
そこで順序プロビット(Ordered Probit)の推定が必要であることがわかったため、Stanでの適用事例を漁っていました。まだマーケティング事例への適用はうまくいってないですが、いったん順序プロビットを簡単にまとめて今後の備忘録としておきます。

順序プロビットとは

被説明変数yが連続潜在変数y∗に対応していると考えるとする。
潜在変数は観察できないが、被説明変数yは観察でき、これらの2つ変数の関係は次のように表される。
(今回扱うデータは3から8までの順序データのため、以下のような表記になる。)

この対応関係は閾値メカニズムと呼ばれている。
各被説明変数をとる確率は以下のように記され、プロビットでは正規分布を扱うため以下のようになる。

これらの選択確率からなる尤度関数を最大にしたものが順序プロビットの推定となる。(c0=−∞でc6=∞とする。σは1とする。)

このように、潜在的な順序関係を想定し、それを満たすように閾値とパラメータを推定する点において、潜在変数を用いたモデルの柔軟性の高さが感じられる。

なぜ順序プロビットを使うのか

マーケティング業務において扱うデータにおいて、NPSやアンケートなど順序尺度の質的変数が多いので、それらのデータを二値データに落とし込んだり、そのまま基数データとして扱うのではなく、適切に扱いたいというモチベーションがあります。加えて、順序尺度の質的変数をもとに予測する際は普通のOLSだと、今回のケースで3を下回ったり、8を超えたりする可能性があり、予測結果として使いにくいです。
アンケートの点数をそのまま被説明変数として回帰しているケースは、データアナリティクスにこだわりの無いメンバーとかであればままあることなので、順序プロビットの民主化というか、布教していきたいと思います。

今回扱うデータ

勝手ながら大人のirisだと思っているワインデータです。今回は赤ワインに絞って、品質に関する順序変数を被説明変数として、各変数との相関を見ていきます。
まずはGGallyパッケージのggpairs()関数を適用して傾向を掴みます。見にくいので是非コードを回して確かめてください。

データに関する説明はワインの味(美味しさのグレード)は予測できるか?(1)で丁寧になされていますので、ご確認ください。

モデル

データセットに含まれる全部を含めて順序プロビットで回帰してみようと思います。
つまり、「酒石酸濃度、酢酸濃度、クエン酸濃度
残留糖分濃度、塩化ナトリウム濃度、遊離亜硫酸濃度
総亜硫酸濃度、密度、pH、硫酸カリウム濃度、アルコール度数」
の全てを使って赤ワインの質への影響を見ていきます。

Stanコード

最初に、Stanのユーザーガイド2.17の138ページにあるOrdered Probitのサンプルコードを使ってみたのですが、
収束しなかったので、初期値を設定するか弱情報事前分布を導入するかの判断が必要となります。
そこで、jabranhamさんが係数が平均0で分散10の正規分布に従うとするサンプルコードを書かれていたので、そちらを使って推定します。
書き換えているところはデータの制約くらいです。

StanをキックするためのRコード

推定結果の可視化を行うためのcommon.Rは松浦さんのGitHubにあるものになります。

結果

まず、MCMCが収束したかどうかの判断ですが、ShinyStanに従うものとします。

ShinyStanによる収束診断をクリアできています。

続いて、推定したパラメータです。

係数の符号がはっきりと分かれている、赤ワインの品質に影響を与えそうな変数としては、volatile acidity(酢酸濃度)、chlorides(塩化ナトリウム濃度)、total.sulfur.dioxide(総亜硫酸濃度)、sulphates(硫酸カリウム濃度)、alcohol(アルコール度数)のようです。

最後に、推定した閾値です。

1と2の閾値が近く、2と3の開きが大きく、あとは比較的均等のようです。

比較

質的変数をそのまま重回帰した際の結果ですが、符号やその大小はあまり変わらないです。やはり予測の際にどの順序尺度の値に対応するかがわかるのが使う利点だと思います。

考察

マーケティングにおいて、順序尺度の質的変数を扱う際に順序プロビットを積極的に使っていきたいと思いますが、アンケート分析を行う際は、ユーザーごとの評点の癖が点数に影響を与えている可能性があります。
そのため、点数の付き方がユーザーごとに違うとする階層モデルへの拡張を今後行っていくのが面白いと思いますし、実際に研究されている論文があります。

参考文献

StanとRでベイズ統計モデリング (Wonderful R)
stan-examples/limited-dv/oprobit.stan
stan-dev/stan users-guide-2.18.0.pdf
Stanによる順序ロジット回帰
第9章 順序選択モデル: 年金投資選択問題
Wine Quality Data Set
ワインの味(美味しさのグレード)は予測できるか?(1)
世界一簡単な収束[シナイ]Stanコード
RStanとShinyStanによるベイズ統計モデリング入門

おまけ

数式をブログに載せる際は、こちら
Online LaTeX Equation Editor – create, integrate and download
でインタラクティブに数式を作成し、その結果を
QuickLaTex Publish Math on the Web without compromising quality
に貼り付けて画像を出力しています。

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の心得シリーズをトピックモデルで分類する。 – 驚異のアニヲタ社会復帰への道