データサイエンティストを超えた!Excelの条件付き書式設ティストに俺はなる!(2)〜作成したマクロをボタンで呼び出せるようにする。

こんにちは、マクロマンです。


前回の続きとして、前回作成した関数をエクセルにてボタンで呼び出す手続きをまとめたのでメモ


前回作成したマクロは業務内容によっては、結構使う事が多い*1です。
わざわざ関数を直接呼び出すのが面倒なのでボタン一つで呼び出せるようにしてしまいます。


1.関数を含むマクロをAddInsフォルダに配置する。
前回作成したRepeatConditionlFormatを含むエクセルマクロアドインファイル(Macro.xlam)を作成します。
以下のフォルダにそのファイルを追加します。

C:\Users\ユーザー名\AppData\Roaming\Microsoft\AddIns



2.Excelのオプションより、アドインを確認する。
「ファイル」⇒「Excelのオプション」⇒「アドイン」に画面遷移します。
以下のようにMacro.xlamが登録されていることを確認し、Macroを選択している状態で設定ボタンを押します。



3.アドインを有効化する。
Macroアドインを有効化します。



終わったらExcelを再起動します。


4.リボンメニューにマクロを追加する。
「ファイル」⇒「Excelのオプション」⇒「リボンのユーザー設定」に画面遷移します。

  1. コマンドの選択を「マクロ」に変更します。
  2. リボンのユーザー設定を「メイン タブ」に変更します。
  3. 「新しいタブ」ボタンでマクロを登録するタブを作成します。
  4. 「新しいグループ」ボタンでグループをマクロを登録するグループを作成します。
  5. マクロ関数を選択した状態で「追加」ボタンをクリックし、該当するグループにマクロを追加します。



5.快適なExcelライフを楽しみます。

たのしい!

*1:クラスター分析後のクラスター特徴把握、アンケートデータの概観を把握するなど

BUGS/Stan勉強会を開催します/Intoroduction of Stan

こんにちは、三十路手前のアナリストです。掲題の通り、9月29日にBUGS/Stan勉強会を開催します。


ATNDリンク BUGS/stan勉強会 #1


内容や趣旨についてはATNDを見てもらうとわかると思いますが、
独学だと理解が危ういのでいろいろな方にご指摘頂ける勉強会のほうがいいだろうというのと、
勉強会形式でアウトプットが蓄積されれば、初めて挑戦される方の敷居も低くなるだろうと考えたためです。


BUGS/stanって何?
BUGS/stanって、何かというとベイズモデリングを行うためのプログラム言語でありその実行環境です。
具体的にはBUGSではこんなことが出来ます(beroberoさんのblog)。


階層ベイズモデルで勝敗データからプロ棋士の強さを推定してみました


stanは比較的新しい同系統のツールです。


開催のきっかけ
私が、BUGS/stanをはじめて知ったのは久保拓弥先生の「データ解析のための統計モデリング入門」です。
この本は単純な回帰分析から階層ベイズまで網羅的に解説している良書であり、なんといっても説明の平易さが光っています。
この本を読んだ当時は、MCMCを使った複雑なモデリングっていらないかなーと思っていたのですが、


以下の2つのきっかけでこの技術を深めたいと考えました。


1.個体差に着目した分析がしたい
仕事柄、いつも知りたいと思うのはデータの個体差だったりするのですが個体差を変数として扱うにはデータが足りずモデル化できず、
集計に基づく分析だとデータの背後にある要因を反映できず、当たり前の結果と外れ値ばかりでつまらない。


といったことに悩まされていました。


先ほどberoberoさんのblogを読んだのが、このタイミングでこの技術は深く掘るべきだと考えました。


2.複雑なモデリングをしたい
同時期に樋口知之先生の「予測にいかす統計モデリングの基本」を読み、
自分の仕事でもこのような複雑なモデリングが出来るとメリットがある。
今後の分析の鍵になる手法だと考えるようになりました。

予測にいかす統計モデリングの基本―ベイズ統計入門から応用まで (KS理工学専門書)

予測にいかす統計モデリングの基本―ベイズ統計入門から応用まで (KS理工学専門書)


なんといっても書籍に出てくる例が飲食店の売上予測で自分が仕事で扱う問題と同じ問題を扱っていたのが大きいです。


というわけで、BUGSを始めたのですが、なんといっても日本語の情報が少なすぎます。
おそらくBUGSについて深彫りするならBUGSBook、stanについては公式ドキュメントがベストのような気がします。
というわけで勉強会を開催しようという経緯になりました。


最後に
開催に当たってはドリコム様に会場提供していただけることになりました。ありがとうございます。
また、私がBUGSに興味を持つきっかけを作ってくださり、会場確保、発表といつも親切に助けてくださっているberoberoさんに感謝です。
私のような技術に疎いな人間でも気軽に質問できることを目指して開催している勉強会です。
皆様、奮ってご参加くださいませ。後一人くらい発表者も募集中です。よろしくお願いします。


Intoroduction of Stan
最後に、勉強会の私の発表スライドを以下に掲載しておきます。stanの入門的解説です。
スライドが英語で大変申し訳ないのですが諸事情により英語で作成しました。
説明はもちろん、日本語で行う予定です。


Windows8にrstanをインストールする。

こんにちわ。BUGSにはまってしまって大変なアナリストです。

どうにもこうにもWINBUGSが遅くていや遅すぎて使えないので、今日はRからstanを叩くrstanをインストールしてみたいと思います。

環境はWindows8 64bitです。

1.Rtoolsをインストール
以下のサイトからRtoolsをインストールしてc++コンパイルができるようにしておきます。
http://cran.r-project.org/bin/windows/Rtools/


適当に次へ次へをクリックしておけば問題ないです。
今回はRtools31.exeをインストールします。


途中環境変数の設定があったりしますが、そのままで問題ないです。

【追記】この時、Windowsでは、PATHは先頭にあるものを優先して実行するのですでにコンパイラーなどをインストールしている場合はPATHの先頭にRtoolsの実行ファイルが含まれているフォルダが存在するのを確認してください。


2.rstanのインストール

次のコードを実行します。

#必要なパッケージのインストール
install.packages('inline')
install.packages('Rcpp')

#hello worldと表示されたら、インストールに問題ない。
library(inline) 
library(Rcpp)
src <- ' 
  std::vector<std::string> s; 
  s.push_back("hello");
  s.push_back("world");
  return Rcpp::wrap(s);
'
hellofun <- cxxfunction(body = src, includes = '', plugin = 'Rcpp', verbose = FALSE)
cat(hellofun(), '\n') 

#rstan自体のインストール
Sys.setenv(R_MAKEVARS_USER='')
options(repos = c(getOption("repos"), rstan = "http://wiki.stan.googlecode.com/git/R"))
install.packages('rstan', type = 'source')

library(rstan) 

Rのflexmixパッケージで混合分布モデルによるクラスタ分析を行う。

Rで混合分布クラスタリングを行うときに有名なパッケージとしてflexmixが存在します。この記事ではflexmixの簡単な使い方を解説します。
flexmix自体は潜在クラス回帰を行うパッケージなのですが、混合分布クラスタリングを行うことも出来ます。
flexmixはRのglmクラスを用いてモデルを表現出来るため、他のパッケージに比べて柔軟なモデリングが可能というメリットがあります。


そもそも、混合分布クラスタリングとはなんぞやという人は以下の本文を参考にしてください。


1.モデルベースのクラスタリングとは
クラスタリングは代表的なものとして、以下の3つの方法が存在します。
おそらくk-meansと階層的クラスタ分析はみなさんご存知でしょう。

分類 メリット・デメリット 手法
階層的手法 +データを樹形図として表現可能
‐データ数が多いと、樹形図として表現できないのでデータ数が絞られる。
階層的クラスタ分析
分割的手法 +代表点を入手可能
‐代表点以外のデータの関係性やモデルに関する知識が得られない。
k-means
PAM
モデルベース手法 +データの発生メカニズムを知ることが出来る。
‐事前にある程度のデータ発生メカニズムについて知る必要がある。
混合分布モデル


混合分布を用いたクラスタリングクラスタリングをする際にモデル構築するため、モデルベースのクラスタリングに属します。
階層的クラスタリングのようにデータ間距離の設定などがいらず、すべて確率という尺度に直した上で、クラスタリングできる。
またデータ発生メカニズムを知ることができるため成功すれば、非常に強力な手法であると言えます。


2.混合分布を用いたクラスタリング
混合分布を用いたクラスタリングについては@sleipnir002さんが、TokyoRで発表された以下のスライドを参考にしてください。



ちなみにこれは以下の本を参考にしているとのことなので、こちらの本の第4章も参考にしてください。


パターン認識 (Rで学ぶデータサイエンス 5)

パターン認識 (Rで学ぶデータサイエンス 5)


3.混合分布一般によるクラスタリング

上記で紹介している手法はパッケージMClustを用いているため、データ発生メカニズムとしては正規分布しか用いることができません。
私は小売向けのアナリストですが、回数のデータなんかはポワソン分布だし、性別比とかは二項分布だし、正規分布だけだと色々と厳しいのが実情です。


そこで、もっと柔軟にモデルを表現する必要があります。


4.Rでの実行例

flexmixは潜在クラス回帰を行うためのパッケージですが、入力するglmのモデル式の形を工夫することで、混合分布クラスタリングを行うことが出来ます。
flexmixを実行して、クラスタ分析を行う例を用いて解説しましょう。

library(flexmix)

#A.ポワソン分布×2と正規分布変数のデータ(y1,y2,y3)を作成
#クラスタ数は3(100, 100, 20)
y1<-c(rpois(100, 3), rpois(100, 8), rpois(20, 6))
y2<-c(rpois(100, 10), rpois(100, 0.5), rpois(20, 6))
y3<-c(rnorm(100, mean=200, 15), rnorm(100, mean=200, 10),rnorm(20, mean=200, 10))
class<-c(rep(1, 100), rep(2, 100), rep(3, 20))
data<-data.frame(y1,y2,y3, class)
rm(y1, y2, y3, class)
#y1
barplot(table(data[, c(4,1)]))
#y2
barplot(table(data[, c(4,2)]))
#y3
plot(data[, 3], col=data[,4])

#B.混合分布モデルのコンポーネントをそれぞれ作成
Model_1<-FLXMRglm(y1 ~ 1,family="poisson")
Model_2<-FLXMRglm(y2 ~ 1,family="poisson")
Model_3<-FLXMRglm(y3 ~ 1,family="gaussian")

#C.EMアルゴリズムの実施
mp<-flexmix(. ~ ., data=data, k=3, model=list(Model_1, Model_2, Model_3), control=list(verbose=1, nrep=5))
# Classification: weighted 
# 1 Log-likelihood :   -2405.1507 
# 2 Log-likelihood :   -2173.0186 
# 3 Log-likelihood :   -1936.5001 
# 4 Log-likelihood :   -1911.7593 
# 5 Log-likelihood :   -1905.8872 
# 6 Log-likelihood :   -1903.8647 
# 7 Log-likelihood :   -1902.4038 
# 8 Log-likelihood :   -1901.2357 
# 9 Log-likelihood :   -1900.3630 
# 10 Log-likelihood :   -1899.6985 
# 11 Log-likelihood :   -1899.1289 
# 12 Log-likelihood :   -1898.5728 
# 13 Log-likelihood :   -1897.9860 
# 14 Log-likelihood :   -1897.3486 
# 15 Log-likelihood :   -1896.6569 
# 16 Log-likelihood :   -1895.9206 
# 17 Log-likelihood :   -1895.1609 
# 18 Log-likelihood :   -1894.4037 
# 19 Log-likelihood :   -1893.6678 
# 20 Log-likelihood :   -1892.9590 
# 21 Log-likelihood :   -1892.2749 
# 22 Log-likelihood :   -1891.6134 
# 23 Log-likelihood :   -1890.9800 
# 24 Log-likelihood :   -1890.3919 
# 25 Log-likelihood :   -1889.8763 
# 26 Log-likelihood :   -1889.4597 
# 27 Log-likelihood :   -1889.1529 
# 28 Log-likelihood :   -1888.9460 
# 29 Log-likelihood :   -1888.8157 
# 30 Log-likelihood :   -1888.7373 
# 31 Log-likelihood :   -1888.6915 
# 32 Log-likelihood :   -1888.6652 
# 33 Log-likelihood :   -1888.6502 
# 34 Log-likelihood :   -1888.6417 
# 35 Log-likelihood :   -1888.6369 
# 36 Log-likelihood :   -1888.6342 
# 37 Log-likelihood :   -1888.6327 
# converged

#D.サマリを表示する。
summary(mp)
# Call:
#   flexmix(formula = . ~ ., data = data, k = 3, model = list(Model_1, Model_2, Model_3), control = list(verbose = 1, 
#                                                                                                        nrep = 5))
# 
# prior size post>0 ratio
# Comp.1 0.436   94    135 0.696
# Comp.2 0.445  100    117 0.855
# Comp.3 0.119   26    194 0.134
# 
# 'log Lik.' -1888.633 (df=14)
# AIC: 3805.265   BIC: 3852.776

#E.パラメータを取得する。
parameters(mp)
# [[1]]
# Comp.1.coef.(Intercept) Comp.2.coef.(Intercept) Comp.3.coef.(Intercept) 
# 0.9826601               2.1248895               1.7702220 
# 
# [[2]]
# Comp.1.coef.(Intercept) Comp.2.coef.(Intercept) Comp.3.coef.(Intercept) 
# 2.3623702              -0.6625998               1.6267241 
# 
# [[3]]
# Comp.1     Comp.2     Comp.3
# coef.(Intercept) 200.23493 202.094819 200.380162
# sigma             15.13504   9.412856   5.797705

#F.交差表の作成
table(clusters(mp), data$class)
# 
# 1  2  3
# 1 91  0  3
# 2  0 98  2
# 3  9  2 15

#G.結果のプロット
plot(mp)


A.ポワソン分布×2と正規分布変数のデータ(y1,y2,y3)を作成
今回、y1、y2、y3という3つの変数をクラスタ分析にかけます。
変数は以下の分布に従うとします。

変数名 分布
y1 ポワソン分布
y2 ポワソン分布
y3 正規分布


例えば、小売の顧客分類で例えると、y1は来店回数、y2商品の購入数、y3は利用金額(対数正規に従う)対数正規に従う利用金額の平均が200ってどうなの・・・とか言わないでください。といった例が考えうるでしょうか。


また、このデータには3つのクラス(クラスタ)がいて、それぞれ以下のようなパラメータを持っているとします。

変数名 n(サイズ) y1~λ y2~λ y3~μ y3~σ
クラス1 100 3 10 200 15
クラス2 100 8 0.5 200 10
クラス3 20 6 6 200 10


クラス別のy1のヒストグラムが以下の図です。


クラス別のy2のヒストグラムが以下の図です。


y3の値を各個体別にプロットしたが以下の図です。
(左からクラス1、2、3と続いています。)


B.混合分布モデルのコンポーネントをそれぞれ作成
ここでは、flexmixに入力するモデル式をglm形式で、作成しています。
ポイントは回帰式の説明変数が定数で、被説明変数が各変数yに対応しているということです。


C.EMアルゴリズムの実施
flexmix関数でEMアルゴリズムを実施して、混合分布モデルによるクラスタリングを行います。
kは想定しているコンポーネントの数を、modelには先ほど作成したglmモデルをリスト化して代入します。
controlはアルゴリズムのチューニング用です。
convergedと表示されていれば収束しています。


D.サマリを表示する。
summaryで実行結果のサマリを閲覧することが出来ます。


E.パラメータを取得する。
parametersで推定された各クラスタ毎のモデルのパラメータを見ることができます。
リストの1から順に、y1(Model_1)、y2(Model_2)、y3(Model_3)のパラメータを返しています。
(ポワソン分布のパラメータは対数です。)

実際のクラス別のパラメータ

変数名 n(サイズ) y1~λ y2~λ y3~μ y3~σ
クラス1 100 3 10 200 15
クラス2 100 8 0.5 200 10
クラス3 20 6 6 200 10


推定されたクラス別のパラメータ

変数名 n(サイズ) y1~λ y2~λ y3~μ y3~σ
クラス1 94 2.671553 10.61608 200.23493 15.13504
クラス2 100 8.371972 0.5155094 202.094819 9.412856
クラス3 26 5.872157 5.087182 200.380162 5.797705


F.交差表の作成
clusters(mp)でデータの属するクラスの値を取得出来ます。
これを実際のdataの属するクラスと比較しています、クラス3に多少の誤差はあるものの、分類できているようです。


G.結果のプロット
結果をplot関数に投げると、各クラスのクラス事後確率のヒストグラムが表示されます。


というわけで、flexmixの簡単な使い方を紹介ました。
細かなパラメータ設定についてはマニュアルを参照してください。
また、ガイドが結構丁寧なので、この記事を読まれて興味を持たれた方はぜひこちらを参照してください。

参考


FlexMix: A General Framework for Finite Mixture Models and Latent Class Regression in R
FlexMix Version 2: Finite Mixtures with Concomitant Variables and Varying and Constant Parameters

ドラマの視聴率を予測する(1)〜初回視聴率を予測して、キャストの視聴率を算出する。

こんにちは、ソーシャルビッグデータサイエンティストです。


最近のデータ解析ブームの流れにのって何か面白いデータ解析でもやってみようと思い、データサイエンティストこじらせた結果、
「今まで機械学習的な方法論使ってないところに機械学習的な方法論使うと面白いんじゃないか。」と思って
今回ドラマの視聴率で予測モデルの構築を試してみたのですが、結論を先に言うと、ニューラルネットワークやRFとか、小難しい方法尽く試した結果、失敗しました。


残念なクオリティですが、放送年度と時間帯と言った基本的な数値だけで予測モデルを構築、メモとして、一旦ブログに残しておきます。


1.データの取得
データについてはこちらのドラマ視聴率を参考に、足りないデータはWikipediaで補完して取得しました。

取得データ:ドラマの各話視聴率、キャスト、タイトル、また、Wikipediaを参考にドラマのジャンルを取得しました。
対象期間:2003年1月期〜2013年1月期
対象放送局:日本テレビ、TBS、テレビ朝日、フジテレビ
対象放送時間帯:20〜23時台

ドラマのジャンルは1ドラマが1ジャンルに所属するように以下のジャンルにWikipeddiaのページの内容を
元に機械的に推論したものを手作業で修正をかけて取得しました。
ジャンルの振り方については恣意性は否めないですが、どうせ有意差なかったので、いまやどうでもいいです。



2.ドラマの視聴率に与える要因を考えてみる。
ドラマの視聴率に影響を与える要因って色々なものが考えられると思うのですが、多分以下の様な要因かなと思っています。


ドラマの視聴率に影響を与えそうな要因
年度(長期トレンド)、4半期(周期的要因)、放送曜日、局、時間、キャスト、脚本家、同時間帯の前ドラマの最終回視聴率、ジャンル、面白さ、裏番組の強さ


ドラマの面白さや脚本家の効果が少なくなるように、連続ドラマの第1話の視聴率を予測する変数としています(理由は後述)。
今回は、最後の3つを除いてデータを取得していることになります。


3.機械学習モデルの構築と失敗
最初に作った予測モデルは、ランダムフォレストの回帰木を用いました。
時間帯やある程度主要なキャストにのみフォーカスして、各キャストの出演を表すダミー変数を作成しました。
時間帯や年度ジャンルを表すダミー変数も作成して、m=1000程度の予測変数でデータ数n=500のドラマの第1話視聴率を予測してみましたが、色々無理でした。

その後、PLSやリッジ回帰等も試して見ましたが、決定係数は0.3くらいで全然ダメ。結局、別の方法を考えてることに。



4.分析の基本、散布図と箱ひげ図に立ち返る。
やっぱりキャストの数含めて圧倒的にダミー変数が多すぎるのが問題であるので、モデル構築は一旦忘れて、
分析の基本、散布図と箱ひげ図に立ち返って、各要因を眺めてみたいと思います。
次の表はドラマの視聴率のヒストグラムです。



正規分布というよりはポワソン分布っぽい形に見える気がします。


・時系列的な要因(年度と4半期)
Rのdecompose関数を用いて、ドラマの初回視聴率を加法モデルでトレンド要因と4半期ごとの季節要因と誤差に分けたのが以下の図です。



以下がR上で出力したトレンド成分の元データです。

>dets$trend
         Qtr1     Qtr2     Qtr3     Qtr4
2003       NA       NA 14.14831 14.22333
2004 14.33354 14.70851 14.82425 15.05177
2005 15.44164 15.47893 15.55655 15.66856
2006 15.43191 15.16733 15.01629 14.67898
2007 14.38083 14.31211 14.02714 13.68280
2008 13.57845 13.47857 13.34959 13.15587
2009 12.93400 12.67909 12.36387 12.50027
2010 12.75232 12.79947 13.08385 13.01115
2011 12.73606 12.83365 12.86635 12.69154
2012 12.55217 12.29805       NA       NA


よくドラマの視聴率は下がっていると言われますが実際に下がっているようです。


2005年第4四半期の視聴率15.6685%

2012年第2四半期の視聴率12.2980%


と最高時の視聴率から最新の視聴率が3.3%も下がっています。


総務省の情報通信白書によると、時間帯別TV全体の視聴率が

2005年11月43.2%(21:00)

2012年11月40.2%(21:00)

なので、ドラマの下がり具合がTV全体の視聴率の低下同程度となっています。


これは、23時のドラマが2007年頃から増えたため、ドラマの四半期別視聴率が
最新のものほど下がっているのではないかと思ったのですが、


以下の20時代から22時台の年度別平均(ドラマの初回視聴率)でも、同様の傾向が見られたので、
どうやらそうでもないようで、ドラマ自体が厳しい状況があるようです。

> aggregate(x$rating.1, list(x$year), FUN=mean)
   Group.1        x
1     2003 14.00392
2     2004 14.93478
3     2005 15.46279
4     2006 15.10213
5     2007 14.25455
6     2008 13.30339
7     2009 12.08305
8     2010 12.38966
9     2011 12.44000
10    2012 11.58033
11    2013 11.50667


・放送時間帯(放送曜日+放送局+時間)
放送曜日と放送局と時間帯をあわせて放送時間帯と定義しました。平たく言うと「フジの月9」みたいなアレです。


以下は放送時間帯別初回視聴率の箱ひげ図です(日テレ木曜24時は横軸にありますがデータは存在しません)。



以下は各時間帯の平均初回視聴率からドラマ全体の平均視聴率を引いたものです。
こうして見ると各時間帯が持っているブランド価値みたいなのがよく見えて生々しいですね。
下位にランクインしているものはことごとく放送を終了しています。

> sort(timezoneeffect, decreasing=TRUE)
  フジ 月 21    TBS 日 21   フジ 木 22   フジ 火 21 日テレ 土 21 テレ朝 水 21    TBS 金 22   フジ 火 22 
   5.3181673    3.3284878    1.8840884    1.2284878    0.4233596    0.2490006    0.2104042    0.1981673 
日テレ 水 22    TBS 日 22 テレ朝 木 21    TBS 木 21   フジ 水 21 日テレ 火 22    TBS 木 22 テレ朝 金 21 
   0.1310519   -0.2843327   -0.6600084   -1.0109994   -1.3843327   -1.6718327   -1.9748089   -2.1474906 
   TBS 土 20    TBS 月 20    TBS 水 21   フジ 日 21 テレ朝 金 23   フジ 土 23    TBS 水 22 日テレ 月 22 
  -2.5570600   -2.9954438   -3.1843327   -3.3065549   -3.3787771   -3.7986184   -4.1843327   -5.5843327 
テレ朝 日 23 
  -6.2843327 


以下はおまけで放送曜日、放送局、時間のそれぞれの箱ひげ図です。





流石は俺達のフジテレビさんだ!


・前回のドラマの最終話視聴率
前回のドラマの最終話視聴率を横軸に、新ドラマの視聴率の初回視聴率を縦軸に取ったのが以下のプロットです。
ある程度の相関が見られるようですが、外れ値なんかも結構有りそうなので、
前処理を丁寧に行えば予測精度を上げられたんじゃないかと今は後悔しています。



・ジャンル
ドラマのジャンル別の視聴率が以下の通りです。


下記で結構強引な予測モデルを構築していながら、分散分析で有意差なんかも見ているのですが、残念ながらジャンルの効果だけありませんでした。
グッバイ!俺のジャンル集め作業、4時間。


6.上記の4つの要因を数字に置き換えて予測モデルを構築する。
ここでちょっと強引なことをします。先ほどの、各4半期での季節要因を除去したトレンドの視聴率*1、放送時間帯別視聴率の偏差、ジャンル別の平均視聴率偏差、前回ドラマの最終話の視聴率*2を予測変数として、
ドラマの初回視聴率を被説明変数とする、線形回帰モデル*3を構築します。モデルのマリーはこんなところです。

> summary(m)

Call:
lm(formula = rating.1 ~ trendp + timezonep + genrep + rating.lastdrama, 
    data = bcdf.model)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.6682 -1.9878 -0.1363  1.4224 14.2938 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      12.90161    0.39895  32.339  < 2e-16 ***
trendp            0.72435    0.12642   5.730 1.72e-08 ***
timezonep         0.91070    0.05926  15.368  < 2e-16 ***
genrep            0.03528    0.19275   0.183   0.8548    
rating.lastdrama  0.06130    0.02954   2.075   0.0385 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 2.939 on 512 degrees of freedom
Multiple R-squared: 0.4611,	Adjusted R-squared: 0.4569 
F-statistic: 109.5 on 4 and 512 DF,  p-value: < 2.2e-16 


うーん、イマイチ。まぁ、要らない変数もあったみたいですが、細かいことは気にしません(え)。
これを元にモデル構築に用いなかった2013年度1月期のドラマの初回視聴率を予測した結果が以下の図です。



グラフにしてしまうと、予測モデルとしてイケている感じしかしないのは僕だけでしょうか?
ちなみにこの場合、このグラフの見方としては斜線の右下にあるのは、予測値に対して実績値が下回っているので、
2013年1月期にいちばんイケてないのは・・・おや誰か来たよう(ry


ここで、予測された初回視聴率はドラマのキャストや脚本家によらない、ドラマの初回視聴率ということになると思います。


7.キャストによらない初回視聴率からキャストの視聴率を算出する。
さらにここで強引なことをします。ここで予測された視聴率を実際の視聴率から引いたもの(偏差)をキャストが持っている視聴率と考えます。
(実は初回視聴率を使っているのは、初回視聴率であれば脚本家による脚本の面白さみたいな効果が出てこないと考えたためです。このキャスト視聴率を他の分析に活用したいと考え初回視聴率を予測対象としました。)



数式として意味のあるものではないですが、脳内推論として


初回視聴率実績 = 時間帯 + 放送年度 + ジャンル + 前回ドラマの最終回視聴率 + キャスト + 誤差

初回視聴率実績 - キャストによらない視聴率予測 = キャスト + 誤差


こんな感じの展開をしています。


この偏差を主要5キャストで分配します。これを各俳優で過去にこれまでに主要5キャストとして出演したドラマによる平均を取ります。
例えば、俳優XがドラマA(主演)、ドラマB(2番手)、ドラマC(4番手)で出演している場合、俳優Xの視聴率は


俳優Xの視聴率=(
(ドラマAの初回視聴率実績-ドラマAの初回視聴率予測)×[主演の分配重みつけ]+
(ドラマBの初回視聴率実績-ドラマBの初回視聴率予測)×[2番手の分配重みつけ]+
(ドラマCの初回視聴率実績-ドラマCの初回視聴率予測)×[4番手の分配重みつけ]
)/3


となります。


まぁ、ややこしく話しましたが単純にキャストの過去実績から時間帯や主演年度など考えうる効果を排除した上での視聴率を按分、集計してキャスト視聴率としているだけです。


以下はキャスト視聴率の高い俳優BEST30です。

  no.木村拓哉 no.松嶋菜々子   no.福山雅治   no.吉岡秀隆 no.矢田亜希子   no.柴咲コウ   no.深津絵里 
         3.08          2.68          2.09          1.55          1.51          1.40          1.17 
no.京野ことみ   no.勝村政信     no.平泉成     no.赤西仁   no.生瀬勝久   no.岡田准一     no.阿部寛 
         1.08          0.99          0.99          0.94          0.86          0.84          0.83 
no.仲間由紀恵   no.時任三郎       no.小雪 no.綾瀬はるか   no.八嶋智人 no.黄川田将也   no.高木雄也 
         0.82          0.78          0.77          0.75          0.73          0.73          0.72 
    no.堤真一   no.江口洋介     no.渡辺大   no.賀集利樹     no.松尾諭     no.筧利夫   no.市原隼人 
         0.70          0.66          0.65          0.64          0.64          0.63          0.62 
  no.内野聖陽   no.鹿賀丈史 
         0.60          0.60 

まぁ、知らない人もいますが・・・1位がキムタク兄さんで安心しました*4


全部のキャスト視聴率を見てみたい放送関係者の方がいればデータをあげるのは公開しなければ無問題なので、知り合いの方を介してください。


8.最後に
このキャスト視聴率を用いて今後の分析をしようと思っているのですが、ちょっとこの予測モデルを作っている間に無理ゲーっぽいようで心が折れたので、
尻切れトンボですが、今回はここまで。ちなみにキャスト視聴率を予測モデルに入れた際の効果をk交差検証法で検証してみたところ、平均して、
決定係数5.5〜6程度で推移していたので、結構いい線いったのかもしれませんが、失敗ですね。


飽きてきたので、もう少し扱いやすいデータをしばらく触ってみて、気が向いたらこのデータをもう一度触ってみようと思います。
俺達の戦いはまだハジマッタばかりだ!

*1:値のつかない箇所は端の値で補完します。

*2:前回ドラマがないものは当該時間帯の平均視聴率で補完しています

*3:glm関数でもモデル色々変えて試したけど、大して結果変わらなかったので、線形回帰の結果を使用しています。

*4:個人的にはキムタク兄さんダミー変数は初期の機械学習モデルから残して予測モデルを作っても良かったような気もします。キムタク兄さんだけで残りの誤差、全部説明出来るような気がします。

PostgreSQLの集約関数で来店日付の間隔をらくらく集計

ある論文によると、顧客が優良顧客であるかそうでないかの定義に、顧客の来店間隔の標準偏差を用いることができるそうです。


また、RFM分析をやっていても、Rは現在からのスナップショットなので、一ヶ月に一回しか来ない人に対して「もう、一ヶ月も来ていないので離反している!」なんて言っても、もともと一ヶ月に一回しか来ていないので、そんなことは言えません。顧客の平均来店間隔ととの兼ね合いで見たいなんてことになります。


で、この平均来店間隔、来店間隔標準偏差(分散)、顧客の売上金額や顧客の来店回数なんかに比べるとSQLで実装するととても面倒です。
(平均来店間隔は来店回数の集計日数合計による割り算でできますが。。。)


PostgreSQLには独自の集計関数を定義する。CREATE AGGREGATE文があります。
こいつとPostgreSQLの配列機能、両方組み合わせると、相当便利で先ほどの平均来店間隔や来店間隔分散も簡単に計算できます。


そこで、顧客の平均来店間隔と来店間隔分散をPostgreSQLのファンクションとして実装しました。

--convert date to date[]
DROP FUNCTION date_aggregate(date[], date);
CREATE FUNCTION date_aggregate(date[], date) RETURNS date[] AS '
	SELECT array_append($1, $2)
' LANGUAGE SQL;

--after aggregate date, make it to interval and calculate avg
DROP FUNCTION fin_date_avg(date[]);
CREATE FUNCTION fin_date_avg(date[]) RETURNS numeric AS '
	select AVG(a.a - b.a) from
		(SELECT row_number() over(order by a), a from unnest($1) a order by a) a,
		(SELECT row_number() over(order by a), a from unnest($1) a order by a) b
	where a.row_number = b.row_number + 1
' LANGUAGE SQL;

--after aggregate date, make it to interval and calculate variance
DROP FUNCTION fin_date_variance(date[]);
CREATE FUNCTION fin_date_variance(date[]) RETURNS numeric AS '
	select variance(a.a - b.a) from
		(SELECT row_number() over(order by a), a from unnest($1) a order by a) a,
		(SELECT row_number() over(order by a), a from unnest($1) a order by a) b
	where a.row_number = b.row_number + 1
' LANGUAGE SQL;

--make aggregate function:interval average
DROP AGGREGATE date_interval_avg(date);
CREATE AGGREGATE date_interval_avg (date)(
  SFUNC = date_aggregate,
  STYPE = date[],
  FINALFUNC = fin_date_avg
);

--make aggregate function:interval variance
DROP AGGREGATE date_interval_variance(date);
CREATE AGGREGATE date_interval_variance (date)(
  SFUNC = date_aggregate,
  STYPE = date[],
  FINALFUNC = fin_date_variance
);


実際に使うときは以下の通り、date_interval_avg、date_interval_variance、に売上表の来店日付を渡すだけ。
金額の集計なんかと同様に関数として顧客の来店間隔を計算することができます。
もっと賢いやり方があるかもですが、ご参考までに。便利!

--to get customer's interval avg and interval variance
SELECT 
 c.customerid, 
 date_interval_avg(distinct f.saledate),
 date_interval_variance(distinct f.saledate),
 sum(f.salesamount)
 from sales f, customer C
where f.customerid = c.customerid AND f.saledate between '2012/04/01' AND '2013/03/31' 
group by c.customerid;

協調フィルタリングについてまとめてみた。

A Survey of Collaborative Filtering Techniques(Xiaoyuan Su and Taghi M. Khoshgoftaar, 2009,Advances in Artificial Intelligence)

仕事で協調フィルタリングについて調べる必要が出てきたのだが、あまりよい日本語の文献を見つけられなかったため(後にしましま先生の文献を見つけた)やむなく英語の論文を検索したところ、
上記のよいサーベイ論文を見つけた。というわけでこのサーベイ論文に書かれていることに自分なりに調べたことを加えて、自分用にまとめておく。
また、一部の人達の間ではとても有名なしましま先生の論文(ドラフト版)があるので、英語が苦手な人はそちらをご覧になるとよいと思われる。


協調フィルタリングは、一言で言えばユーザとアイテムのマトリックスを用いた顧客への商品のレコメンデーションです。このマトリックスより、ユーザの相関を分析し、類似したユーザはお互いが購入している商品買うという仮定に基づきレコメンデーションする仕組みといえます。Wikipediaの英語版にはこの直感的な雰囲気を説明した以下の図があります。



Wikipediahttp://en.wikipedia.org/wiki/Collaborative_filtering)より
1.ユーザとアイテム(商品)の関係がネットワークで表現されています。
2.先ほどの関係が行列で表現されています。
3.一番下のユーザ(アクティブユーザ)のパソコンに対する評価を知りたいとします。
4.アクティブユーザと類似するユーザ(緑のユーザ)を選び出します。
5.類似するユーザの評価よりアクティブユーザはパソコンを評価しないと予測します。



1.協調フィルタリングの歴史
この箇所は論文にはそれほど記述がありませんが、自分なりに調べたものです。私が学生だった頃はNetFlixのコンテストが行われ、日本でも必要のない多くのサイトで協調フィルタリングが導入されていたような記憶があります。当時は「Web2.0」とか「ロングテール」なんて今ではめっきり聞かない単語が流行ったり、「集合知プログラミング」が出版されて周りの学生の間でも簡単に実装可能なWebビジネスのアルゴリズムにほくほくしていました。

1992年:Tapestryというコメンデーションシステムの論文において協調フィルタリングという言葉が初めて登場
1994年:GroupLens、ネットニュースのレコメンデーションシステムの論文が出版
1997年:Amazon協調フィルタリングを導入
2007年:NetFlixPrize、優勝賞金$1,000,000のオープンレコメンデーションコンペの開催


2.協調フィルタリングの特性とその抱える問題
協調フィルタリングはその実ビジネスとの関連性の強さと特徴的なデータ構造より、以下のようなこの分野に特徴的な問題点を持っています。

2.1.データのスパース性
2.2.スケーラビリティ
2.3.類似性
2.4.灰色の羊
2.5.シリングアタック


2.1.データのスパース性・・・通常のビジネスで用いる協調フィルタリングでは巨大なアイテムセットを用いるためユーザ・アイテム行列は極端に疎であり、精度の面で問題になります。データがスパースであるために以下のような問題が引き起こされます。このスパース性について実感するにはAmazonの人気商品ランキングを見て、自分が買ったことのある商品を数えてみるとよいでしょう。殆どの人は大半の商品を買ったこと無いはずです。

1.コールドスタート問題・・・新しいアイテムやユーザが追加された時に類似のアイテムを見つけるのが難しい問題。
2.少カバー率問題・・・ユーザの評価が少ないアイテムは類似するアイテム等のレコメンデーションの提示が不可能になること。
3.同類推移問題・・・スパースなデータベースの場合、類似のユーザであっても、全く同じアイテムを共に評価しないと類似であると判別されない問題。

これらの問題に対処するために以下のような方法が用いられます。

1.特異値分解による次元削減
2.ハイブリッド協調フィルタリングによる補完
3.予測モデルの構築によるデータの補完

2.2.スケーラビリティ・・・多くの場合、アイテム数やユーザ数は驚異的なスピードで増大するために、古典的な協調フィルタリングアルゴリズムを使用するのには問題があります。O(n)のアルゴリズムであっても、計算量が大きすぎます(えっ)。いくつかの解決法が存在しますが、レコメンデーションの精度と計算量は常にトレードオフの関係にあります。
2.3.シノニム性・・・シノニム性は同じあるいはとても似ているアイテムが異なるアイテムとなっている傾向のことを指します。シノニムの存在はレコメンデーションの性能を下げます。LSIなどを用いてシノニムを発見することができますが、適合率においては高い性能を示すものの、再現率は低いようです。
2.4.灰色の羊・・・複数の人々と一致するするような曖昧なユーザで、協調フィルタリングの恩恵にあずかれない人を指します。Claypoolはコンテンツベースフィルタリングと協調フィルタリングのハイブリッドアプローチでこの問題を解決する方法を提案しました。
2.5.シリングアタック・・・自己への利益誘導目的の不正なアイテムに対する評価を指します。要はステマです。アイテムベースの協調フィルタリングはユーザベースの協調フィルタリングよりシリングアタックの影響を受けにくいそうです。



協調フィルタリングの分類
協調フィルタリングは大きく、以下の3つの分類に分けることができます。


3.メモリベースの協調フィルタリング・・・ユーザとアイテムの相関関係より、類似したユーザ同士の評価を基にレコメンデーションのルールを作成します。
4.モデルベースの協調フィルタリング・・・統計的なアルゴリズムを用いてユーザの購買行動をモデル化し、モデルに基づいてレコメンデーションのルールを作成します。
5.ハイブリッド協調フィルタリング・・・コンテンツベースのフィルタリングと協調フィルタリングを混合させてレコメンデーションのルールを作成します。


下の表は論文より、上記の協調フィルタリング技術をまとめた表を翻訳したものです(Su&Khoshgoftaar, 2009)。

3.メモリベースの協調フィルタリング
メモリベースの協調フィルタリングはユーザとアイテムのデータベースを用いて、これを行列として表現しユーザ(アイテム)の類似性により、レコメンドする商品をピックアップする方法です。
 

以下の2つはさらに詳細な分類です。

ユーザベースの協調フィルタリング・・・「同じ商品を買っている類似ユーザ同士の評価は似ている」と仮定してユーザの類似性にもとづき推薦する商品を決定します。
アイテムベースの協調フィルタリング・・・アイテムの類似性を中心に計算する手法。Amazonで用いられている手法もこれであり(Wikipediaより)、スケーラビリティの問題に対処するために開発された。

おそらく協調フィルタリングといえば、一般にはモデルベースではなく、こちらを指すことが多いように思われます。

3.1.類似性の計算
3.2.予測とレコメンデーションの計算
      3.2.1.重み付き予測による評価予測
      3.3.2.Top-Nレコメンデーション


アノテーションについて
先ほどのWikipediaの図で表示されていたように、協調フィルタリングではユーザアイテム間の行列を用います。以下のようにユーザがアイテムに対して5段階評価をしている行列を例に考えます。



この時、全ユーザの集合をU、その要素をu,vとします。


全アイテムの集合をI、その要素をi,jとします。


この時、あるユーザuのアイテムiに対する評価を


ユーザuの評価の平均を

とします。


3.1.類似性の計算
先ほどのWikipediaの図で言うところの3の作業に値します。ユーザ(アイテム)の距離、類似度を定義する計算を各アイテム別(ユーザ)別に行う処理を指します。協調フィルタリングにおいては致命的な計算になります。
(他の論文読んでいるとどうも、距離の指標はアイテム−ユーザ行列のスパース性や、評価がユーザによる5段階評価か商品の購買やWebページの閲覧と言ったバイナリの値をとるかといったことに依存するみたいです。)


類似性には下のようないくつかの指標が用いられます。

1.ピアソン相関係数・・・二つの変数の互いに(線形の)関係性の強さを表す度合い(最も代表的な類似性)
2.強制ピアソン相関係数・・・ピアソン相関係数の計算において平均ではなく中央値を用いる指標
3.スピアマン順位相関係数・・・ピアソン相関係数の順位版
4.ケンドールのτ相関係数・・・スピアマンの順位相関係数を相対的な順位に直したもの
5.コサイン類似性・・・ベクトルの内積(ドキュメントの単語頻度による距離に用いられる)
6.調整済みコサイン類似性・・・ユーザーにより、異なる評価スケールを行っていることより、コサイン類似性において平均を引いて標準化したもの、実はピアソン相関係数と同じ。


これらの指標を用いて、アイテム間もしくはユーザ間の類似性を定義します。
ピアソン相関係数の場合は、ユーザuとvの類似性は以下のようになります。

この類似性を用いて、以下のステップでレコメンデーションを作成します。


3.2.1.予測(レコメンデーション)計算
ユーザにレコメンドするアイテムを得るために、類似性に基づいてアクティブユーザのアイテムに対する評価を集計し、予測します。以下のような2つの指標がが存在します。

1.類似ユーザ評価の重み付き和・・・ユーザベースの協調フィルタリングでは、特定のユーザの予測を行う際に、同じ嗜好を持っているユーザはより大きな重みを付けてこの評価を行います。
2.重み付き和・・・アイテムベースの協調フィルタリングで用いられます。


アクティブユーザaに対する、類似ユーザ評価の重み付き和としてのアイテムiの評価の予測式は以下のようになります。

これでアクティブユーザに対する商品への評価の予測ができます。


3.2.2.Top-Nレコメンデーション・・・上記の予測だけでもレコメンデーションはできますが、代表的なレコメンデーションの手法であるTop-Nレコメンデーションは類似しているユーザやアイテムを絞って、そこからN個の商品を実際にユーザに推薦する方法の1つです。


3.1.ユーザベースのTopNレコメンデーション・・・商品ベクトル空間上でのユーザの類似性を測定します。あるユーザiさんに商品をレコメンドする際は、そのうち上位のk人を抽出して、その中でユーザiさんがまだ買っていない商品のうち最も購入頻度の高い上位N件をレコメンドします。この手法は計算量的に限界があります。
3.2.アイテムベースのTop-Nレコメンデーション・・・アイテムベースのTop-Nレコメンデーションはアイテムベースのレコメンデーションは計算のスケーラビリティ、早さにおいてユーザベースのものよりメリットがあります。このアルゴリズムでは事前にアイテムの類似性に基づいて全てのアイテムに対して、レコメンデーションの候補となる上位k個の類似アイテムを計算しておきます。そこからユーザuがすでに購入している商品を取り除き、その集合の中のアイテムをiとすると、uの購入しているアイテムjに対する評価とiとjの類似性で重み付けをした重み付き和を計算して、uへのjの評価を予測します。


4.モデルベースの協調フィルタリング
モデルベースの協調フィルタリングについては理解があやふやなので、羅列するだけにとどめておきます。


4.1:ベイズネット協調フィルタリング・・・ベイジアンネットワークを用いた協調フィルタリング
 4.1.1:単純ベイズ協調フィルタリング・・・ナイーブベイズを用いた協調フィルタリング
 4.1.2:NB-ELRとTAN-ELR・・・単純ベイズより、欠損の多いデータに強いアルゴリズムのようであるが、詳細は不明。性能は、ナイーブベイズやピアソン相関によるメモリベースのモデルよりも強力である。
4.2:クラスタリング協調フィルタリング・・・クラスター分析を協調フィルタリングに用いる手法。多くの場合、クラスタリングはさらなる協調フィルタリングのプロセスの中間ステップとして実行される。例えば、クラスタ間でのメモリベースの協調フィルタリングを行う等。多く場合、クラスタリングは計算量でメリットがあるが精度は低い。
4.3:回帰による協調フィルタリング・・・回帰モデルを構築を行い、特定ユーザの特定のアイテムへの評価を予測する方法。
4.4:MDP(マルコフ決定プロセル)モデルによる協調フィルタリング・・・シミュレーションのような手法・・・らしい。
4.5:潜在意味協調フィルタリング・・・潜在意味解析を用いた手法・・・らしい。


5.ハイブリッド協調フィルタリング
ハイブリッド協調フィルタリングとは、異なる2つの協調フィルタリングの仕組みを混ぜることです。
典型的には、コンテンツベースのレコメンデーションとメモリベースの協調フィルタリングを混合します。
コンテンツベースのレコメンデーションとは、テキストによる文脈情報やユーザの設定やプロフィール等の情報を用いてそこから見つけた規則性によるレコメンデーションの仕組みのことです、


例えば、コンテンツブーステッド協調フィルタリングはナイーブベイズ手法を用いて、ユーザの属性情報(年齢、性別等)とメモリベースで用いてきたユーザアイテムマトリックスを用いて欠損値を補い、擬似評価マトリックスを作成します。
その後にピアソン相関係数に基づき協調フィルタリングを行い欠損値の予測を上書きして、元々存在していた評価は上書きせずに予測を行います。この手法はコールドスタート問題を回避するのにも有効なようです。


他にも、その他利用可能なすべてのレコメンデーション手法を重み付きで混合する手法や、メモリベースの協調フィルタリングとモデルベースの協調フィルタリングを混合する手法等が存在します。
これらの方法の混合により、高い精度を維持しながら、先の2で記載した協調フィルタリング特有の問題に対処するようです。


6.協調フィルタリングの評価
協調フィルタリングの評価には一般的な機会学習の予測評価と同じ指標が用いられます。平均誤差や、再生と再認率、k交差法によるROC曲線、これについては省略します。



最後に
Tescoなんかは結構前からレコメンデーションのDMとかをやっているみたいなので、チラシやDMが大好きな日本の小売でもレコメンデーションを使ったDMが流行る日も近いと思われる。
TokyoR#30ではこれらの解説に加えてRで協調フィルタリングを実行するrecomenderlabの実行例を加えてもう少し具体的な計算手法に突っ込んで解説を行う予定です。
詳細な計算についてはまたBlogに残しておきたい。