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の水準が違っているというような情報と、その会話におけるラベル(会話の結果としての商談の成否)について分析することで、商談などがうまくいく際のヒントが得られるかもしれないですね。

参考情報

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

Word2Vecを用いて蒙古タンメン中本の口コミ評価を予測してみる

はじめに

word2vecを用いた分類は以前からやってみたいと思っていたのですが、関心を持てるテキストデータがなかったのでなかなか手を出していませんでした。
ある時、ふとしたことから某グルメ系口コミサイトから蒙古タンメン中本の口コミと評価点を抽出して、その評価をword2vecでやってみるのは面白いだろうと思いついたので、さっそくやってみます。
こういう時にはじめて、データ分析だけでなくクローリング屋としても業務をやっていて良かったなと思うところですね。
コードは以前見つけて紹介した「分散表現を特徴量として文書分類するための方法について調べてみた」のものを再利用します。

目次

・目的
・データ収集
・形態素解析
・集計
・分散表現とは
・word2vecについて
・gensimのword2vecの引数
・word2vecによる文書分類の適用
・終わりに
・参考情報

目的

某グルメ系口コミサイトの口コミを収集し、個々人の口コミの内容から個々人の店に対する評価が高いか低いかを予測する。

データ収集

BeautifulSoupで収集しており、各店舗あわせて数千件ほど集めました。(実行コードはこちらでは紹介しません。)

このようなデータが手に入っている前提で以下の分析を進めていきます。

形態素解析

文書を形態素解析して、名詞のみを抽出するためのコードを用意します。

先ほどのデータフレームに対して以下のように実行すれば、名詞のみの分かち書きを行ったカラムが手に入ります。

集計

点数のヒストグラム

3.5点から4点の間が最も評価が多いようです。1点台をつける人はほとんどいないことがわかります。

単語数のヒストグラム

大体の口コミで100単語未満のようです。

単語数と点数の散布図

どうやら口コミにおいて500語を超える記述をしている人は評価が3点を下回ることはないようですが、文字数と点数でキレイに傾向が出ているわけではないですね。

形態素解析結果の集計、単語ランキング

名詞の抽出に関して非常に便利なMeCab Neologdを用いています。蒙古タンメンもきちんと捉えることができています。

味噌よりも北極の方が出現しているようです。北極は言わずもがな、極端に辛い罰ゲームレベルの一品。味噌タンメンは辛さが抑えめのラーメンで、知人の間では最もおいしいのがこのレベルだという合意があったりしますね。

分散表現とは

  • 単語の意味を低次元の密な実数値ベクトルで表現したもの。
  • 入力層から中間層への重み自体が各単語の分散表現となっている。
  • 2017年9月のテキストアナリティクスシンポジウムにてメルカリとGunosyが特徴量として分散表現を活用しており性能が出ているとの発言があった。

word2vecについて

単語の分散表現を作ることを目的としている。

  • CBOW(Continuous Bag-of-Words)
    注目している単語の前後N単語を文脈と呼び、その文脈をBag-of-Words表現として入力し、注目している単語を出力するというニューラルネットワークを学習する。入力層から隠れ層への結合は単語の位置を問わず同じとし、隠れ層の活性化関数をただの恒等関数としている。
  • Skip-gram
    文脈のBOWを突っ込むCBOWとは異なり、入力層に1単語だけを入れる。1単語を入力し、正解データとして他の単語を入れることを繰り返して学習し、ある単語の入力に対して、どの単語の出現確率が高いかどうかを計算する。正解確率が上がるようにニューラルネットワークの重みを調整する。深層学習で使われる自己符号化器と似たような構造とされている。

gensimのword2vecの引数

gensimのword2vecには数多くの引数が存在します。gensimのドキュメントに英語で書かれていますが、せっかくなのでこちらで紹介します。

  • sentences
    解析に使う1行1センテンスで書かれた文書。日本語の場合はLineSentenceフォーマットを使えばうまくいった。単語が空白文字で区切られていて、文章は改行で区切られていれば問題ない。
  • sg
    {1,0}の整数で訓練アルゴリズムを設定できる。 1を選べばskip-gramで、0ならばCBOWを使う。
  • size
    特徴ベクトルの次元を設定する。
  • window
    文書内における現在の単語と予測した単語の間の距離の最大値を設定する。言い換えると、文脈の最大単語数を設定する。
  • alpha
    学習率の初期値を設定する。
  • min_alpha
    訓練の過程で徐々に落ちていく学習率の最小値を設定する。
  • seed
    乱数を生成する際のシード番号を設定する。
  • min_count
    一定の頻度以下の単語を除外する際の値を設定する。
  • max_vocab_size
    語彙ベクトルを構築している際のメモリ制限を設定する。
  • sample
    (0, 1e-5)の範囲で、頻度語がランダムに削除される閾値を設定する。高速化と精度向上を狙っており、自然言語処理においても高頻度語はストップワードとして除去するなどの対応が取られている。
  • workers
    モデルを訓練するために多くのワーカースレッドを利用するかどうか設定する。(並列化関連)
  • hs
    {1,0}の整数で、1であれば階層的ソフトマックスがモデルの訓練で用いられ、0であり引数negativeがnon-zeroであればネガティヴサンプリングが設定できる。全部計算することが大変なので、階層的なグループに分けて各グループごとに学習するというのがモチベーション。
  • negative
    0よりも大きければネガティブサンプリングが用いられる。5〜20などを選び、どれだけノイズワードが描かれているかを識別する。0であればネガティブサンプリングが適用されない。ネガティブサンプリングは計算高速化を目的に出力層で正解ニューロン以外のニューロンを更新しないように学習する手法。
  • cbow_mean
    {1,0}の整数で、0であれば単語ベクトルの合計を用い、1であればCBOWが用いられた際の平均が用いられる。
  • hashfxn
    訓練の再現性のためにランダムに初期値のウエイト付けできる。
  • iter
    コーパスにおける繰り返し回数(エポック数)を設定できる。
  • trim_rule
    ある単語を語彙に含めるべきかどうかを識別する、語彙のトリミングルールを設定する。
  • sorted_vocab
    {1,0}の整数で、1であれば頻度の降順で語彙を並べ替える。
  • batch_words
    ワーカースレッドにわたすバッチの大きさを指定する。(並列化関連)
  • compute_loss
    Trueであれば損失関数の計算と蓄積を行う。
  • callbacks
    訓練時の特定の段階で実行する際に必要なコールバックのリストを指定できる。

word2vecによる文書分類の適用

口コミの点数が4点以上であれば1、そうでなければ0を取る変数を作成し、それをラベルとして文書分類を行います。
以前、紹介したブログ同様に、scikit-learnのExtraTreesClassifierを用いてCountVectorizerとTfidfVectorizerを特徴量としたものをベースラインとして、同様の手法に対してword2vecで作成した分散表現を特徴量として用いたものとを比較します。評価指標はクロスバリデーションスコア(5-folds)とします。

分類の前に、せっかくword2vecを使ったので、任意の単語に類似した単語を見てみます。

まずは初心者向けの味噌ラーメン

続いて、中級者向けの蒙古タンメン

そして、上級者向けの北極ラーメン

最後に、誰もが経験する翌日という単語。

どれも関連性の高いと思われる単語が抽出できているように思われます。

それでは分類モデルの学習を以下のコードで行います。
scikit-learnを使えば、データさえあれば非常に短いコードで書けてしまいます。

一応、ベースラインよりもword2vecを特徴量としたものの方がスコアが高いのですが、わずかです。TF-IDFベースで特徴量を作成したモデルは十分に性能が出ているようです。
word2vecを用いることによる旨味はそれほどなさそうですが、パラメータを試行錯誤していけばよくなるかもしれません。

終わりに

蒙古タンメン中本のテキストをWebスクレイピングし、その口コミ情報をコーパスとして口コミ評価の二値分類に挑戦しましたが、TF-IDFよりもわずかに優秀な特徴量になりうるという結果になりました。もっと劇的な向上を夢見ていたのですが、パラメータの試行錯誤を今後の宿題としようと思います。

参考情報

Chainer v2による実践深層学習
word2vecによる自然言語処理
models.word2vec – Deep learning with word2vec
Python で「老人と海」を word2vec する
Python3 – MeCabで日本語文字列の名詞出現数の出し方について
Transform a Counter object into a Pandas DataFrame

「NOSQLの基礎知識」を読んで基礎知識を養う

同じ職場のエンジニアに『NOSQLの基礎知識』はいいぞ〜と勧められて読もうとしているのですが、マーケターレベルだと、2章と4章を読めば良いとのことです。それ以外は実際に手を動かす人が必要な知識らしい。なので、今回は1、2、4章に書かれている内容を簡単にまとめてみようと思います。

NOSQLとは(1章)

  • Not Only SQLの略で、SQLではないと言う意味。
  • RDBは単体のハードウェア上で利用するには最適だが、複数台並べてデータを分散して管理するのが不得意だがNOSQLはそれが得意。
  • SQLは膨大かつ多種多様で激しく変化するデータの塊を制限なしに保存することが難しいが、NOSQLはそれを可能にする。

NOSQLの特徴(2章)

  • 機能がRDBほど豊富ではない
    Joinなどの演算子もサポートされていない。(今はMongoDBではlookupとかいう機能が付いている模様)
  • データの整合性がゆるい
  • 結果整合性で良いという考え
    一時的なデータの不整合は問題としない。
  • NOSQLは正規化、更新書き込み量の最小化、メモリキャッシュの利用などの専門性が不要。
  • SQLとNOSQLは補完関係にある。
    SQLは豊富な機能とデータの整合性に強いが、ビッグデータに弱い。NOSQLはその逆。

NOSQLの種類(2章)

  • 世界に100種類以上存在
    ビッグデータに関連するカテゴリは以下の4つ
    ・キーバリュー
    ・カラム指向
    ・ドキュメント指向
    ・グラフ型

  • 一般的なRDBなどのSQL
    将来どのようなデータが発生するか想定の範囲内で設計が行われる。

  • キーバリュー
    ・キーとバリューの対が行として格納される。(データモデルがシンプル)
    ・キーを指定すればバリューを探し出せるのでサーバが何台増えても対応でき、スケールアウトに適している。
    ・複数のサーバに複数のバリューを複製した場合、複数のバージョンが存在する可能性がある。
    ・具体的なNOSQL:Dynamo,Voldemort,Riak,Hibari,Redis,Scalaris,Tokyo Cabinet/Tyrant

  • カラム指向
    ・キーの付された行が複数のカラムを持つことを許容している。
    ・あるカラムにおいて全てにデータが格納されていなくても良い。
    ・カラム数を後から増やしても良い。
    ・具体的なNOSQL:Bigtable,Cassandra,HBase,Hypertable

  • ドキュメント指向
    ・JSONやXMLといったデータ記述形式で記述されている。
    ・各ドキュメントは階層構造を持たず、相互の関係が横並びに管理されている。
    ・具体的なNOSQL:CouchDB,MongoDB

  • グラフ型
    ・データとデータ間のつながりを管理できる。
    ・個々のデータの属性に依存することなく、相互のつながり方だけを独立した形で保存できる。
    ・ノード(頂点)、リレーションシップ(関係の有無、関係の方向)、プロパティ(属性)の3つの要素から構成。
    ・属性は、期間や有効/無効のデータなどを持たせることができる。
    ・データ間の関係性を示すグラフを作成するという目的のために開発され、最適化されている。
    ・具体的なNOSQL:InfiniteGraph,Neo4j

HadoopはNOSQLではない(4章)

  • Hadoopは複数のソフトウェアを包含するフレームワーク
  • HDFSは大量のデータを複数のサーバに分散して格納する点でNOSQLのように見えるが、MapReduceによるバッチ処理前提となっており、リアルタイムでの応答の速さを求められる対話型のNOSQLとは分けて捉える必要がある。

MapReduce処理(4章)

  • MapReduceは大規模なデータを複数のノードで並列分散処理するためのプログラミングパターン
  • データ処理をMapタスク(膨大な元データを分解し必要な情報を抽出し有用な形に変換して出力)とReduceタスク(抽出された情報を集約し一塊のデータとして出力)の2段階のフェーズで分けて行う。
  • 多くのNOSQLにもMapReduceは実装されている。

さいごに

これまでR、Python、SQLしか心得ていなかった自分の知識の狭さを、エンジニアと一緒に働くことで実感しました。一冊読むくらいでエンジニアの話に少しはついていけると思うと読書の幅は広げておきたいですね。今後、データを活用したマーケティング支援ツールの開発を先導する際に、技術選定などの意思決定に役立つかもしれません。データアナリストという職種である者からすると、マーケティング・エンジニアリング・データサイエンスを緩急つけて学んでくというスタンスが重要なのかもと思う今日この頃です。

参考文献

NOSQLの基礎知識 (ビッグデータを活かすデータベース技術)

データアナリストもLinuxについて学んでみる

はじめに

もともとデータアナリストとして働いているのですが、最近、エンジニアと一緒に仕事をすることが増えたので、良いとされている教科書を読むなどしています。その中で、『新しいLinuxの教科書』がわかりやすく書かれていたので、今回はそこで出てきたコマンドたちを箇条書きにして忘備録としておきたいと思います。なのでいつもよりショボいです。
普段、Rをメインとして機械学習・クローリング目的でPythonを扱うものとしては、エンジニア寄りの知識が薄いなぁと課題に感じていたので、これを良い機会と思ってCLIへの耐性を付けていきたいと思います。

感想

分析業務において、仮想環境でCentOSやUbuntuを扱うことはあるのですが、分析環境を作った上で、クローリング業務・RPAシステムなどのスケジュール実行をするくらいで基礎からちゃんと勉強したわけではありませんでした。そのため、教科書には知らないコマンドもありました。
この本はLinuxを使って開発する際の作法(移植性大事、GUI使うな、十字キー使うな)から、歴史的なコマンドの変遷、シェルスクリプトを用いた作業の効率化、Gitなどのチーム開発ツールの紹介まで幅広く扱われていて十分な内容だと思いました。図による解説もあるので、普段R使っているような非エンジニアでも読めるはずです。これからエンジニアと一緒に分析業務をしていくデータアナリストの方には強くおすすめできる本だと思います。Rでの当たり前の機能が、Linuxのシェルで普通にあるんだと色々と学びがありました。

コマンドたち

date

現在の日時の表示

echo

文字列の表示

pwd

現在のディレクトリを確認

cd

カレントディレクトリの変更

ls

ディレクトリ内のファイルの表示、オプションでパターン指定も可能

mkdir

ディレクトリを作成する、pオプションで深い階層も一気に作れる

touch

ファイルを作成する

rm

ファイル・ディレクトリを削除する

rmdir

空のディレクトリを削除する

cat

ファイルを表示する

less

ファイル内容をスクロール表示する

cp

ファイル・ディレクトリをコピーする

mv

ファイル・ディレクトリを移動する

ln

ファイルにリンクを張る

find

ディレクトリツリーからファイルを探す

locate

ファイル名データベースからファイルを探す(高速)、何日か前のファイルを探すの適している。

which

コマンドのフルパスを表示する

alias

エイリアスの設定(わかりやすい別名を与える)

set

bashのオプションを切り替える

shopt

bashのオプションを切り替える(種類豊富)

PATH

コマンド検索パス

printenv

環境変数の表示

export

環境変数の設定

sudo

コマンドをスーパーユーザとして実行

ps

プロセスの表示(x, ux, ax, aux, auxww)

jobs

現在のジョブ一覧を表示

fg

ジョブをフォアグラウンドにする(入力受付状態にする)

bg

ジョブをバックグラウンドにする

kill

ジョブ・プロセスを終了させる

|

パイプライン(処理と処理をつなぐ)
R使いからするとRのdplyrでパイプ処理を知ったので、Linuxに普通にあるのを知れて新鮮な気分になった。

head

先頭の部分を表示
こちらもRでよく使うものですね。

du

指定したファイルやディレクトリの使用容量を表示する

wc

入力ファイルの行数・単語数・バイト数を数える

sort

行単位でテキストを並び替える

uniq

同じ内容の行を省く

cut

入力の一部を切り出して出力する(元データから特定の部分を抽出したい際など)

tr

入力の文字を置き換える・削除もできる(-d)

tail

ファイルも末尾部分を表示する

diff

二つのファイルの差分を表示する

grep

文字列を検索する

sed

テキストの編集を元ファイルを変更せずに行える

awk

テキストの検索や抽出・加工などの編集操作を行う(sedより高機能で、スクリプトは「パターン」と「アクション」からなる)
・列選択もできる
 ・CSVファイルからのスコア集計などもできる。
 ・あんちべさんの本『データ解析の実務プロセス入門』でも紹介されていた。

xargs

標準入力として引数のリストを与える(findやgrepとの組み合わせで使うことが多いらしい)

tar

ファイルをアーカイブファイルにまとめたり、アーカイブファイルから元のファイルを取り出したりする

gzip

ファイルを圧縮・展開する(一つのファイルしか圧縮できないのでtarとの併用を行う)

bzip2

gzipよりも高い圧縮率で圧縮を行うが、gzipよりも圧縮・展開に時間がかかる

zip

アーカイブと圧縮を同時に行う(多くのディストリビューションで標準インストールされていない。unzipも必要)

git

ファイル変更履歴を保存し管理するツール
 ・git init リポジトリを作成
 ・git add コミット対象として登録
 ・git commit コミットする
 ・git status ワークツリーの状態を表示
 ・git diff 差分を表示する
 ・git log 履歴を表示する
 ・git revert コミットを取り消す
 ・git branch ブランチの一覧表示、新しくブランチを作成
 ・git checkout ブランチを切り替える
 ・git merge 指定したブランチを現在のブランチにマージする
 ・git branch -d ブランチを削除する
 ・git push 変更の履歴を送る
 ・git clone リポジトリを複製する
 ・git remote add リポジトリパスに別名を設定する
 ・git fetch 他のリポジトリの変更を取り込む

yum

パッケージ管理(CentOS)
 ・yum install パッケージのインストール
 ・yum erase/remove パッケージのアンインストール、削除
 ・yum search パッケージの検索
 ・yum info パッケージの詳細情報を表示

apt

パッケージ管理(Ubuntu)
 ・sudo apt-get install パッケージのインストール
 ・sudo apt-get remove パッケージの削除
 ・sudo apt-get purge 設定ファイルも含めた完全削除
 ・apt-cache search パッケージを検索

ssh

リモートログインする

info

オンラインマニュアルを表示する(索引や階層構造を設定できるなどmanより充実)

参考文献

新しいLinuxの教科書

LIMEで赤ワインのデータをいじってみる with Python

はじめに

2018年1月のTokyoR( TokyoR67に行ってきました )で機械学習結果の解釈可能性について多く語られていたので、Hello World的な試行を赤ワインのデータで行ってみたいと思います。コードは論文を書いた人のGitHubのものを拝借しました。

営業やマーケティングのメンバーに機械学習手法を提案する際に、ひとつひとつの予測結果において、なんでその予測結果になったのか理由が知りたいという要望を受けることが多いです。
ある講演で、実装用のモデルにランダムフォレストやSVMを使い、マネジャーに説明する用に決定木の結果を見せたとかいう話もありました。わざわざモデルを二つ作って説明のための工数を取らなくてよくなると思うと非常にありがたい技術です。事業側の解釈可能性を重視するあまり、シンプルな手法を取らざるを得ない現場には朗報ですね。

目次
・LIMEとは
・データ
・コード
・結果の解釈
・参考情報

LIMEとは

Local Interpretable Model-agnostic Explanationsの頭文字をとったもので、機械学習によって構築したモデルに関して、その予測結果を人間が解釈しやすくする技術です。
流れとしては、

  • まずランダムフォレストなりXGBoostなりで分類器を作る。
  • 任意のデータxを取り出し、解釈可能バージョンのx’(x’∈{0,1}で、xの非ゼロ要素を1としている。)を用意する。
  • x’の周辺のデータをサンプリングする。
  • サンプリングしたデータを使って、元のモデルを近似するために、モデルの距離を目的関数とした最適化問題を解く。
     (K-LASSOという線形モデルの手法を使うことで、最終的に近似するモデルの変数の数を決めている。)
  • 学習したモデルの偏回帰係数を確認して、予測結果への影響度を見る。

という流れのようです。間違っているかもなので、参考情報をご覧になってください。
紹介動画も作られているようです。そういえば、Stanも紹介動画ありましたね。

データ

ワインデータから赤ワインのデータのみを利用します。データの詳細はこちらにあります。
http://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality.names

このようなデータセットでデータサイズは1599です。

ワインの質に関するヒストグラムはこんな感じです。

今回の分析の目的

赤ワインの評価値が7以上であれば良いワイン(Y=1)、そうでなければ悪いワイン(Y=0)として、評価値が7を超えるような赤ワインの分類器を学習させ、任意のテストデータを取り出し、そのテストデータが良いワインである要因を探るものとします。

進め方

・LIMEのインストール(pip install lime で一発)
・scikit-learnによる機械学習(ランダムフォレスト)
・LIMEを用いた予測結果の解釈の提示

コード

今回のコードはこちらにもあります。
kamonohashiperry.com/lime_study/Lime_With_Wine Data.ipynb

結果を解釈する以前に、二値分類モデルとして精度が低かったら元も子もないので、精度やAUCを確認してみます。

scikit-learnの引数でclass_weight=”balanced”にしているので、少しは不均衡データに対応できているようです。AUCは8割を超えたかったですが、いったんこれで進めます。

結果の解釈

こちらのコードで、ランダムにインスタンスを選んで、その予測結果とその予測結果に影響を与えている変数を見てみます。

どうやら、このインスタンスを良いワインと予測しており、sulphates(硫酸)が0.74を超えていた、volatile acidity(酢酸とその派生物質)が0.39よりも小さいというのが理由のようです。画像では結果が消えていますが、exp.as_list()で結果を抽出できます。

著者のコードには続きがあって、このインスタンスに対して、値を足したりすることで分類確率がどのように変わるのかが記されていました。今回はアルコールを3%ポイント増やして、総亜硫酸濃度を30増やしてみるものとします。

これを見る限り、アルコールを3%ポイント増やすと、悪いワインの確率が6.2%ポイントあがり、総亜硫酸濃度を30増やすと、悪いワインの確率が49.8%ポイントあがることが示されています。
ワインあまり飲まないので結果の解釈以前に変数の解釈ができていないのは今回のオチになるんでしょうか。このコードをもとに、仕事現場で使ってみようと思います。きっと解釈できるはず。

参考情報

“Why Should I Trust You?” Explaining the Predictions of Any Classifier
lime/doc/notebooks/Tutorial – continuous and categorical features.ipynb
3.2.4.3.1. sklearn.ensemble.RandomForestClassifier
機械学習と解釈可能性 by Sinhrks
機械学習モデルの予測結果を説明するための力が欲しいか…?
LIMEで機械学習の予測結果を解釈してみる
ワインの味(美味しさのグレード)は予測できるか?(1)