ベイジアンネットワークを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でタグを指定して抽出します。

データ取得後は簡単にpandasのstr.replaceで整形すると、以下のような結果になります。今週の順位と先週の順位が引っ付いてしまっています。

ここから横着してRで整形し、各週の順位データを作成しました。

データ確認

データ構造はこのような形です。

まず、どんな楽曲やアーティストがランキングに入っているのかを簡単に確認してみます。 ほとんど聞いたことない人の名前、曲名ですが。。

続いて、2010年8月〜2017年6月の間に100位以内に入った数を楽曲ごとにヒストグラムにしてみます。べき乗分布な形かと思いきや、20回前後で盛り上がっているのが気になりますね。

中央値が9週間なので、意外と長い期間Top100には入っています。

続いて、100位以内に入った数をアーティストごとにヒストグラムにしてみます。こちらはべき乗分布のような形になっています。

Top10入りの楽曲の実態

Top10に入っている楽曲のみに絞って、ヒストグラムを描いてみます。

Top10に入ったら、10週近くは10位以内に含まれるようです。上位はすぐに取って代わられるのかと思いきや、人気が人気を呼ぶとかなのでしょうか。確か、某ECサイトの方が、生キャラメルは売れるから売れたんだとか言っていた気がします。

順位の推移

100位以内にランクインした回数が最も多かった楽曲のTop10に関して、時系列プロットをしてみます。


初回に上位にランクインして、後は下がるだけの楽曲や、じわじわとランキングを上げていく楽曲などが観察されています。楽曲の消費のサイクルみたいなものがあるのでしょう。

今後について

せっかく面白そうなデータが手に入ったので、リンク先も辿って、どのような楽曲やアーティストの特徴が人気に繋がりうるのか見てみるのも良いですね。あと、もう少し洋楽聴いてみようと思います。私はクラシック音楽とジャズしかウォークマンに入っていないので。

参考文献

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

ラーメン二郎の某飲食店レビューサイトデータに対して共分散構造分析をしてみる

データ

ラーメン二郎に関して、某飲食店レビューサイトのデータをWebスクレイピングしたもので、「料理・味」・「サービス」・「雰囲気」・「CP」・「酒・ドリンク」の評価項目に関して、1~5の実数値が割り振られています。ラーメン二郎の店舗数40店のうち、欠損のなかった37店舗の評価データとなります。

記述統計量・ヒストグラム

まずは記述統計量とヒストグラムを見てみます。

まぁ、お酒を飲むところではないので、drinkは低いですよね。しかしながら、総じて3点台なのは驚きです。昔行ったことがあって、雰囲気は決して良くはないはずなので。

共分散構造分析

今回、検証したい仮説は「ラーメン店としての質が二郎愛につながるかどうか」です。

ラーメン店としての質に繋がりそうな評価項目
「料理・味」
「CP」
「酒・ドリンク」

二郎愛に繋がりそうな評価項目
「料理・味」
「サービス」
「雰囲気」

加えて、料理・味とCPは関係していそうな項目なので、その点もパスにおいて考慮しておきます。

以下の図は仮説のイメージです。

パスを描写する際の細かい指定は参考文献を参照されると良いと思います。
さっそく、モデルを作って実行します。

こちらが推定結果です。

lavaan (0.5-23.1097) converged normally after 100 iterationsとあるので、適切に推定されたようです。自由度が2とデータ数が少ないためかなりギリギリな推定となっています。ただ、推定すべき母数の数よりもデータ数が一応多い状態ではあります。二郎インスパイア系の店のデータも集めた方がいいかもしれないですね。

モデルの評価

モデルの評価として適合度と母数の推定に関して見ていきます。

適合度
・適合度指標であるCFI(Comparative Fit Index)が1なので、適合度に関しては良さそうです。
・同じく適合度指標であるTLI(Tucker-Lewis Index)が0.998なので、適合度に関しては良さそうです。
・0.05以下であれば当てはまりが良いとされるRMSEAは0.028なので、当てはまりは良さそうです。
・0に近いほど良いとされるSRMRは0.024となっています。

母数の推定

f2 ~ f1はラーメン店としての質と二郎愛の関係を想定したものですが、p値が0.47と全然だめでした。ラーメン店としての質が二郎愛に繋がるという仮説は正しいとは言えないです。

参考文献にもあるよう、係数の解釈を行いやすくするために標準化推定値を求めます。

これによると、f1(ラーメン店としての質)に対してはコストパフォーマンスが最も影響を与えるようです。味よりもコストパフォーマンスが勝っているという考察になりますが、それはそれで面白いですね。他方、f2(二郎愛)に対しては味よりも雰囲気・サービスが影響を与えるようです。

semPlotパッケージを用いてパス図を出力してみます。

変数名が3文字に省略されているのですが、修正方法がパッと見つからなかったので、そのまま載せています。

初歩の初歩ですが、一通りの進め方がわかったので、今後も共分散構造分析にチャレンジしてみたいと思います。

参考文献

共分散構造分析 R編―構造方程式モデリング

RStanで学部時代の研究を振り返ってみる

研究概要

大学時代に実験経済学で行った実験結果のデータがUSBに入っていたので、振り返って分析などをしてみたいと思います。

研究目的
 ピア効果に関して、競争相手が自分よりも秀でた人がいいのか劣った人がいいのかを確かめる。

実験方法
・1分間で100マス計算を2セット解いてもらう。(めちゃ速い人には3枚目も渡した)
・実験開始後、実験対象のクラスによって、途中で「平均的なクラスは○○マスまで進んでいます!」とアナウンスします。アナウンスすることで、競争相手のレベルを知り、焦るなり余裕を感じるなりしてもらおうという計画です。
なお、対照群はアナウンスをしていません。アナウンス内容は「平均告知(18秒)、上告知(15秒)、超上告知(12秒)、下告知(20秒)」と4パターンとなります。
・計算が間違っているものは加点しません。

実験対象
 某国立大学の経済学部生の1~2年の必修科目履修者217名(先生に交渉して授業の開始5分を頂いて実験を行いました。)
 内部進学やスポ専などがない分、計算能力的にある程度近い集団ではないかと思われます。

検証方法
 アナウンスごとに100マス計算の点数の水準が変わりうるのかを回帰分析などで判断。

データ可視化

以下、実験カテゴリごとの略記です。
下告知(20秒)・・・slow_20
上告知(15秒)・・・fast_15
超上告知(12秒)・・・fastest_12
平均告知(18秒)・・・average_18
対照群・・・baseline

データ構造の確認です。

平均値、中央値、標準偏差、サンプルサイズを出してみます。

中央値で見てみると、baselineに対してわずかですが点数に違いがありそうに見えます。

実験種別で点数に関するヒストグラムと確率密度関数を確認してみます。


baselineが多峰性がありそうなのが気になります。average_18は低そうに見えますね。

RStanで重回帰

『StanとRでベイズ統計モデリング』にあるコードを参考にしています。正規分布を事前分布にした線形回帰モデルです。
被説明変数が点数、説明変数が実験種別のダミー変数だけからなります。

結果

traceplot(fit)でMCMCのサンプリング結果を確認する。

収束しているように見えます。

以下推定結果ですが、残念ながらベイズ予測区間において符号の逆転が起きていないものはなかったので、アナウンスによる効果があるとは言えないようです。ただ、slow_20の係数がおしいですね。少なくとも他の実験種別よりも、アナウンス効果があるかもしれないという考察に止まりそうです。

分布でも見てみます。アナウンス効果が0を確実に超えているとは言えないですね。

結局、学部時代のレポートと結論は変わらないのですが、係数が0よりも大きい確率という観点で結果に向き合えたのは良かったと思います。

参考文献

StanとRでベイズ統計モデリング (Wonderful R)

PythonやRを用いたアルゴリズム取引・株価分析まとめ

はじめに

まわりでシステムトレードや株価の機械学習による予測などに関心が高まってきたので、私も少し調べてみようと思いPythonやRで行われた分析・実装の事例を集めてみました。自分の資産を突っ込む気にはなれないですが、事例を知っておくだけ知っておきたいですね。

調査法

Google検索において以下のクエリで上位に表示されたサイトを中心にまとめました。
「python 機械学習 株価」
「python 機械学習 為替」
「python アルゴリズム取引」
「python machine learning stock price」
(同様にRも調べました。)

アルゴリズムトレードの理論から学ぶ

システムトレードで億万長者になるぞ! coursera で Computational Investing Part I を受けた
こちらはCourseraで開かれていたシステムトレード向けの講座を受けた方のブログです。Active Portfolio Management: A Quantitative Approach for Producing Superior Returns and Selecting Superior Returns and Controlling Riskはもっとも重要なアルゴリズムトレードの本とされているようです。出版年が1999年と結構古いですが。
Active Portfolio Management: A Quantitative Approach for Producing Superior Returns and Selecting Superior Returns and Controlling Risk (McGraw-Hill Library of Investment and Finance)

Rを用いた株価データのモニタリング

Stock Analysis using R
こちらはアルゴリズム取引とかではないですが、Rで株価を簡単に取得できるquantmodライブラリについて紹介されています。

以上のコードを実行しただけで時系列データがボックスプロットで描写されました。非常に便利そうです。複数のグラフを並べてディスプレーで見てみたいですね。図ではテクニカル分析で使われるボリンジャーバンドを一発で出してくれています。傾向を掴むための分析に役立つのではないでしょうか。

MSFT_stockPrice

RSIで株価の連動性を見る

過去のデータからビッグデータ分析で株価を予測する
ここでは、主にpandasを用いて株価の分析を行っています。RSI(Relative Strength Index)という「一定の期間変動幅の中でどれ位株価が上昇しているのか、下落しているのかをはかるもの」を計算して日経平均との相関を見ています。データさえあれば、個別株の市場連動性を見るぶんにはpandasで十分に分析できそうです。

アルゴリズムトレードのシステム構成

自動トレードボット
こちらでは、仮想通貨取引のための自動トレードボットの作成のための手順などが書かれています。
工程としては
・仮想通貨価格データ取得
・バックテストの実施
・明日の価格予想
・学習パラメータの最適化
・結果のメール送信
・APIを使った成行トレード
・ジョブのスケジューリング
などとなっているようです。
完全自動化しようと思うと、作るのは大変そうですね。

WEB屋の自分が機械学習株価予想プログラムを開発した結果
こちらはよりシンプルで
「チャートを見る限り、誰がどうみても今日値上がりする銘柄」を検索して、毎朝Slackに通知してくれるシステム
をPythonで作成されています。最終的に資産を溶かす形になっているようです。

アルゴリズムトレード向けのPythonライブラリ

pythonのアルゴリズムトレードライブラリ
こちらでは
zipline( http://www.zipline.io/
PyAlgoTrade( http://gbeced.github.io/pyalgotrade/ )
pybacktest( https://github.com/ematvey/pybacktest )
backtrader( https://github.com/mementum/backtrader )
というライブラリが紹介されています。

遺伝的アルゴリズムで為替の自動売買

pythonと遺伝的アルゴリズムで作るFX自動売買システム その1
遺伝的アルゴリズムでFX自動売買 その2 進化する売買AIの実装
遺伝的アルゴリズムでFX自動売買 その3 OandaAPIで実際に取引
こちらは、外国為替の取引をPythonで自動化させた試みです。最終的に負けてしまってはいますがシステム構成についても詳しく書かれているので、勉強になりました。データはOANDAというサービスが提供しているAPIを用いて取得し、GA(遺伝的アルゴリズム)を使って為替の売買タイミングを決めているようです。実際の売買にもOandaAPIというものを利用して完全に自動化させています。

決定木による株価リターンの予測

機械学習で未来を予測する – scikit-learn の決定木で未来の株価を予測
こちらはPythonでscikit-learnを用いて、決定木による株価の予測をされています。目的変数としてはリターンインデックスをおいています。前処理においてpandasが使われています。

ランダムフォレストによる株価の予測

Pythonで機械学習を使った株価予測のコードを書こう
こちらでは、ランダムフォレストを用いた機械学習で、ETFなどのデータを特徴量にして個別の株価予測を行っています。

アカデミックサイドでの研究事例

Stock Price Forecasting Using Information from Yahoo Finance and Google Trend
UC Berkeleyの経済学部生の研究です。こちらはR言語で、ヤフーファイナンスとGoogleトレンドの情報を用いて株価を予測している研究です。これまでの時系列手法よりもパフォーマンスが良いそうです。

最後に

事例を集めてみて、どこまでのレベルのものに手を出すべきか悩ましいと思いました。完全自動化は維持するコストがかかりそうですし、ロスが増大しないか心配です。自分としては情報を自動で取得し、リターンが発生する確率の高そうな銘柄をサジェストしてくれるレベルで十分な気がします。

追記

人工知能学会の全国大会2017でAIを用いた株式運用に関する研究がなされていました。深層学習を用いたものが多そうです。
株価変動パターンの類似性を用いた株価予測
深層学習と高頻度データを用いた株式注文状況の推定
LSTM を用いた株価変動予測
深層学習を用いた株価予測の分析

Rで株価の時系列データを簡単に集計する

はじめに

かねてからYahoo!ファイナンスなどの情報に、株価について期間を指定したリターンや標準偏差や回帰分析のあてはまりなどが無いので、そのようなデータを簡単に見れるようにしたいという願望がありました。加えて、できるだけ多くの銘柄の情報をみたいという願望もありました。Rを使えばその全てが簡単に実現できるのでチャレンジしてみました。

・データ収集
・データテーブルの整形
・成長率の計算
・標準偏差の計算
・自己回帰の計算
と非常にシンプルな工程ですが、3600銘柄分のこれらの情報を一気に手に入れれるというのは魅力的だなぁと思います。

データの収集

RFinanceYJというパッケージをインストールして、以下のrfinance_module.Rのスクリプトを置いておきます。

以下のコードで2016-01-01から現在に至るまでの、指定した株価コードの株価情報を集めてきます。3600銘柄の株価コードを用意する必要があるので、集めてcodelist.csvで保存しておきます。銘柄コードの情報はここでは掲載しません。(適当なサイトにころがっています。)

以下のコードを実行して、元のデータを調整済み終値について、縦が日付で横が銘柄のテーブルに整形します。

株価テーブル

成長率

指定した期間での成長率を計算します。

return_rate

標準偏差

ここではシンプルに標準偏差の計算を行います。

標準偏差

自己回帰による決定係数

自己回帰でフィッティングの良い銘柄はどれなのかを単純に知りたいので、3600銘柄分計算させます。

AR1の係数

以上より、見てみたかった指標の3600銘柄分を一気に計算することができました。

株価指標テーブル

可視化

試しに一番成長率の高かった株価をプロットしてみます。

sdj7612

おわりに

株価データにカジュアルにアクセスできる環境を手に入れることができたので、今後は様々な指標や時系列解析・業界ごとの分析を試してみたいと思います。

顧客生涯価値(CLV)の計算 with R

顧客生涯価値(CLV:Customer Lifetime Value)を計算してくれるRのコード(Calculating Customer Lifetime Value with Recency, Frequency, and Monetary (RFM))があったので、今更感がありますが取り上げたいと思います。

目次

・顧客生涯価値の数式
・データセット
・関数
・データセットの読み込みと加工
・再購買率とRFMとの関係
・再購買率の推定
・顧客生涯価値の計算
・参考情報

顧客生涯価値の数式

まず、顧客生涯価値の数式は以下の通りです。
customer_lifetime_value

t:年や月などの期間
n:顧客が解約するまでの期間合計
r:保持率(1-解約率)
P(t):t期に顧客から得られる収益
d:割引率

rは数式上では固定ですが、実際にはデモグラ属性(年齢、地理情報、職種など)や行動(RFMなど)や在職中かどうかなどの要因により変わりうるものだと考えられます。参考文献のブログでは、このrのロジスティック回帰による推定がなされています。

データセット

データ名:CDNow
概要:1997年の第一四半期をスタート時点とした顧客の購買行動データ
期間:1997年1月〜1998年6月
顧客数:23570
取引レコード数:69659
変数:顧客ID、購入日、購入金額
入手方法:DatasetsでCDNOW dataset (full dataset)をダウンロード

関数

参考文献にはgetDataFrame関数、getPercentages関数、getCLV関数の三つの関数が出てきますが、CLVの計算に必要なのはgetDataFrame関数、getCLV関数の二つです。getPercentages関数はRecencyなどに応じて細かく分析する際に用います。

getDataFrame関数・・・生のデータセットから、指定した期間に応じたRecencyのデータを作成する関数です。

getPercentages関数・・・Recencyなどの回数に応じて、購入した顧客の割合などを計算するための関数です。

getCLV関数・・・Recency、Frequency、Monetary、購入者の数(1人と置く)、コスト(0としている)、期間、割引率、推定したモデルをもとにCLVを計算する関数です。

データセットの読み込みと加工

再購買率とRFMとの関係

まず初めにデータセットを加工します。ロジットの推定における説明変数用のデータとして19970101〜19980228のデータを用い、被説明変数にあたる購入したかどうかのデータを19980301〜19980430のデータを用いて作ります。

データを確認します。

60日以内に購入した顧客(Recency=0)のうち、45%が再び購入しているようです。

100ドル購入した顧客(Monetary=10)のうち、19%が再び購入しているようです。

10回購入したことのある顧客(Frequency=10)は49%が再び購入しているようです。

再購買率の推定

RFM(Recency、Frequency、Monetary)のデータに基づいて、再購買率をロジスティック回帰によって推定し、予測確率を用いて顧客生涯価値を計算します。

rfm_analysis_est

RecencyとFrequencyによるロジスティック回帰

顧客生涯価値の計算

推定したロジットを用いて、生涯価値を計算します。

それではさっそく、1998年5月〜6月のデータを用いて、今回推定した顧客生涯価値が妥当なのかどうかを確かめたいと思います。

予測したCLVと実際の取引金額データを散布図で描き、回帰線を引く。

life_time_value_estimation

CLVが上がれば、1998年5月1日〜6月30日の間(未来)の取引金額が増すような傾向が出ています。

参考情報

RFM Customer Analysis with R Language
Calculating Customer Lifetime Value with Recency, Frequency, and Monetary (RFM)

Japan.R 2016のスライドまとめ

まだ手に入れていないスライドもあるので随時更新しますが、Japan.R 2016(connpass)のスライドをまとめています。後日、登場したパッケージなどのサンプルコードも載せていく予定です。

目次

・石田基広さんのキーノート
・ホクソエムとは何だったのか(ホクソエムさん)
・Rと探索的データ分析で、国連での日本の立ち位置を可視化する(安田洋介さん)
・マウス操作でかんたん予測分析(鈴木了太さん)
・高速・省メモリにlibsvm形式でダンプする方法を研究してみた(@hskksk)
・Rでてんしょくかつどう(@Med_KU)
・RStudio vs Emacs(@y__mattu)
・randomforestで高次元の変数重要度見る(@siero5335)
・Rで本を作りたい(前田和寛さん)
・28歳でプログラミングを始めた話(市川太祐さん)
・LDA-Visパッケージのご紹介(@doradora09)
・【e2d3R】E2D3からDot-Bar-Chartのご紹介(楠本一哲さん)
・このIRのグラフがすごい!上場企業2016(@ito_yan)
・Rでカルマンフィルタをしたい(@tetsuroito)
・PPAP(仮)(@yutannihilation)
・スライド未公開、ユーザーの状態遷移に関する分析のお話(@sanoche16)
・私とR(高栁慎一さん)
・めくってもめくってもサンプル画像(服部恵美さん)
・木と電話と選挙(causalTree)(安井翔太さん)
・スライド未公開、dplyrの話(@tomomoto)
・てかLINEやってる?(仮)(@wonder_zone)
・心理学における「再現性」の問題とBayes Factor(@NSushi)

・石田基広さんのキーノート

スライド未公開です。

・Linux使い
・ヘブライ語の意味構造を代数学でやっていた
・S/R言語の生みの親はJohn Chambers
 以下の二つは最近書かれた本だそうです。
 Software for Data Analysis: Programming with R (Statistics and Computing)
 Extending R (Chapman & Hall/CRC The R Series)
・S→S-plus→Rの順番で発展
・purrrを最近使い始めたそうです。
・XLConnectパッケージを使って、大学教員の採点活動を効率化しているそうです。

・ホクソエムとは何だったのか(ホクソエムさん)

匿名技術者集団ホクソエムの2016年の成果
・densratio( densratio: Density Ratio Estimation
・githubinstall
githubinstall: A Helpful Way to Install R Packages Hosted on GitHub
・healthplanet( Wrapper package for healthplanet api
・RODBCDBI
RODBCDBI: Provides Access to Databases Through the ODBC Interface
・jpmesh( jpmesh: Utilities for Japanese Mesh Code

起業されたとのことです。懸命に頑張って下さい!
株式会社ホクソエム

awesomeな人材が必要とのことで、awesomeな方はアプライしてみてはいかがでしょうか。

・Rと探索的データ分析で、国連での日本の立ち位置を可視化する(安田洋介さん)

スライド未公開です。
国連のデータを使って、Exploratoryを用いた探索的データ分析の実演をされていました。

・マウス操作でかんたん予測分析(鈴木了太さん)

R AnalyticFlow
Rで実践!データサイエンス~初めの一歩から高度な応用まで~

・高速・省メモリにlibsvm形式でダンプする方法を研究してみた(@hskksk)

・Rでてんしょくかつどう(@Med_KU)

Rmd でreveal.js のhtml スライドプレゼンテーション

・RStudio vs Emacs(@y__mattu)

RStudio vs Emacs Japan.R 2016

・randomforestで高次元の変数重要度見る(@siero5335)

・Rで本を作りたい(前田和寛さん)

Rで本を作りたい

・28歳でプログラミングを始めた話(市川太祐さん)

・医療関連のアプリ開発でデータサイエンスを駆使しようとしているそうです。

スライド未公開です。
スライドがシェアされ次第載せます。

・LDA-Visパッケージのご紹介(@doradora09)

・【e2d3R】E2D3からDot-Bar-Chartのご紹介(楠本一哲さん)

スライドは未公開です。
E2D3をRで表示する試みのようです。
Experiments with e2d3 in R

・このIRのグラフがすごい!上場企業2016(@ito_yan)

スライド未公開です。後日シェアしていただけるようです。

・Rでカルマンフィルタをしたい(@tetsuroito)

・PPAP(仮)(@yutannihilation)

・スライド未公開、ユーザーの状態遷移に関する分析のお話(@sanoche16)

スライドがシェアされ次第載せます。

・私とR(高栁慎一さん)

RjpWiki
統計・データ解析
統計解析フリーソフト R の備忘録頁 ver.3.1
seekR(R限定の検索エンジン)
からだにいいもの
アブラタニブログってなんでしょう。油谷さんのブログ?

・めくってもめくってもサンプル画像(服部恵美さん)

Rのサンプルコードはあるけれども、どんな図ができるのかはわからない。そこで、サンプルコードとグラフを大量にまとめているサイトを作ったそうです。検索性は未知数ですが、暇なときに眺めておきたいですね。
R Graphical Manual

・木と電話と選挙(causalTree)(安井翔太さん)

・スライド未公開、dplyrの話(@tomomoto)

スライドがシェアされ次第載せます。

・てかLINEやってる?(仮)(@wonder_zone)

・心理学における「再現性」の問題とBayes Factor(@NSushi)

スライドは後日公開とのことです。

『マーケティング・サイエンス入門』に出てくる手法をRで実行してみる

友人に『マーケティング・サイエンス入門』がおすすめと言われて読んだんですが、やっぱり実行できないとモヤモヤしてしまいますよね。そこで、登場する手法に関連したRのコードやらを集めてみました。

・BASSモデル
・多次元尺度法
・因子分析
・ロジット&プロビット
・分散分析
・クラスター分析
・判別分析
・決定木
・コンジョイント分析
・RFM分析
・共分散構造分析

BASSモデル

市場全体の規模が動的にどのように変化するかを予測するために使われるモデル。
R を使ってバスモデルを当てはめてみた – 廿TT
こちらにRのコードや適用例がいくつか載っています。

早速、私も携帯電話の加入契約数の時系列データを用いて、コードを実行してみました。データは平成25年版の総務省の情報通信白書の表から得ました。( 第2部 情報通信の現況・政策の動向
mobile_phone_plot

当てはまりはわずかながら、BASSモデルの方が良いようです。

多次元尺度法

多次元尺度法で遊んでみる(オレ流 R入門)
こちらのブログで山手線の駅間の距離データの可視化がなされています。
各駅ごとの距離からなる行列さえ用意すれば、cmdscale()関数を実行することで可能なようです。

今回はContaminatedMixtパッケージに含まれているワインのデータセットを使って多次元尺度法を適用してみようと思います。

データはこんな感じです。
%e3%82%b9%e3%82%af%e3%83%aa%e3%83%bc%e3%83%b3%e3%82%b7%e3%83%a7%e3%83%83%e3%83%88-2016-11-03-18-48-01

以下のコードで実行しました。

wine_cmd

Barbera(バルベーラ)・・・基本的にはタンニンをあまり含まず、酸味の強い色の濃い赤ワインで庶民的。
Barolo(バローロ)・・・アルコール度数が高く、非常に重厚な味わいのワインでワインの王様と呼ばれる。
Grignolino(グリニョリーノ)・・・僅かにタンニンを感じるサッパリとした辛口の赤ワインで庶民的。

庶民と王様のワインは成分においても違いがありそうですね。

因子分析

psychパッケージというものがあるようです。こちらのサイトを参考にして進めます。( スナック菓子の食感についてRで因子分析してみた
今回は大好きなwiskyのデータセットを使ってみます。( Classification of whiskies

グレンフィディックやカリラやタリスカーがイメージ通りにプロットされています。ラガブーリンやアードベッグがはみ出しているのが残念ですが。
biplot_wiskey

ロジット・プロビット

これらの手法はビルトインの関数でできてしまいますが、せっかくウイスキーのデータがあるので、薬っぽさに繋がりそうな変数を見つけてみます。