日記20180504
日記じゃ
- Character level CNNでアーキテクチャを変える
- 回線のスピードチェックアプリを入れて見た。
- microbit
- ２、３時間のお昼寝
- ヨドバシの前でGUのセールをやってたのでTシャツを３枚買った。
- Trelloのipadアプリでチケットにチケットを添付できない。
- 今日の思い付き
- pycon miniのスライド全然できてないじゃん・・・
Character level CNNでアーキテクチャを変える
原論文（[1509.01626] Character-level Convolutional Networks for Text Classification）のアーキテクチャではどうも上手く学習できないのでアーキテクチャを全然違うものに変えた。
低次元空間に埋め込みした後に複数のカーネルサイズを通して、最後に結合して特徴を出力するやつ。
参考
[1702.08568] eXpose: A Character-Level Convolutional Neural Network with Embeddings For Detecting Malicious URLs, File Paths and Registry Keys
[1408.5882v2] Convolutional Neural Networks for Sentence Classification
character-level CNNでクリスマスを生き抜く - Qiita
CNNを利用したセンチメント分析 - にほんごのれんしゅう
VocabularyProcessor
文字を数値化するのはtensorflow.contrib.learnのpreprocessing.VocabularyProcessorを使ったので日本語も問題なくできる。
注意点としてはこれは元々は単語の処理をする為のものなので最大長を4とかに設定すると
from tensorflow.contrib.learn import preprocessing vocab_processor = preprocessing.VocabularyProcessor(4) docs = ["a aa aaa aaaa", "aaa a"] print(list(vocab_processor.fit_transform(docs)))
出力
[array([1, 2, 3, 4]), array([3, 1, 0, 0])]
なので単語の間にスペースを入れて文字を整数にマッピングするようにした。
gist.github.com
あとcontrib.learn.preprocessingは古いのでtf.dataを使えと言う警告が出るがtf.dataの使い方がわからない。
と言うかドキュメント読んでも同じことが出来るようには思えない・・・ うーむ・・・
方針
学習も上手くいってるのでこのアーキテクチャで色々実験する。
回線のスピードチェックアプリを入れて見た。
なんか回線が遅い気がしたのでせっかくなのでアプリを入れて確認して見た。
speedcheck
これ
チェックした上下の速度だけでなくSSIDや場所も履歴として残るので良い。
色んな場所のFree wifiの速度をチェックしたくなる。
結果をipadやiphoneで同期できるらしいのでアカウントも作った。
ただ広告がウザすぎるので360円くらい払って広告出さないようにした。
また、ipadの方のアプリで支払ったのにiphoneの方でログインしたのにまだ広告出てるのはなんかちょっと微妙やな。
同じアカウントなら一つ買ったら別の端末でも広告出さないようにしてよ。iphoneの方は広告出したままにした。。。
結果
wifiで二種類繋げられるんだけど、チェックしたらなぜか知らないけど遅い方に接続してることが発覚した。
だから遅く感じたのか・・・
速い方に変えた。
速度
- 遅い方
くだり 21.36Mbps
のぼり 5.63Mbps
- 早い方
くだり 59.24Mbps
のぼり 4.69Mbps
どっちも普通に速いやんけ。ってかのぼりは同じあまり変わらないのね。
microbit
twitterで見かけたのでちょっと調べた。
http://microbit.org/ja/guide/
教育用のマイコン。BBCが作っていてイギリスの小学生とかが使ってるらしい。
マイクロusbでパソコンに繋いでプログラミングできるみたいだ。pythonでもコーディング可。
ジャイロ、温度センサー、無線通信、bluetoothと機能は一通り揃ってるみたいだ。
２、３時間のお昼寝
例のごとくMONSTER ENERGYを昨日二本飲んだので夜は3時間くらいしか眠れなかったので２、３時間お昼寝した。
昼寝すると午後からの作業が捗るな。
ヨドバシの前でGUのセールをやってたのでTシャツを３枚買った。
３枚で3000円ちょっとだった。
なんてないTシャツなので写真は割愛
今日の思い付き
あ、プラグラミング言語の特定とかなら出来そうか。
— tdual@情報幾何 (@tdualdir) 2018年5月4日
ブログのネタとしてはいいかもしれないけどあまり実用的ではないかも・・・
pycon miniのスライド全然できてないじゃん・・・
やばい・・・
日記20180503
日記じゃ
- 今日はずっとcharacter level CNNの実装してた。
- 久しぶりにAWSのGPUインスタンスを起動した。
- 『前処理大全』の3-3を読んだ。
- ChromeのdevToolを開くショートカット
- 歌詞取得スクリプト
- コンソールに求人を埋め込むAbemaTV
今日はずっとcharacter level CNNの実装してた。
完全に論文通りなのになんか上手く学習されないのでどこかバグがあるんだろうか。
tensorboardでsparsityを確認したらいきなりほぼ1になっていたのでなんかおかしい・・・
久しぶりにAWSのGPUインスタンスを起動した。
Char level CNNの学習を待って居られないので時間をお金で買うことにした。
だが、Amazon reviewのデータセットを送ったらすぐにディスクが一杯になった。
前にSSDの学習のための画像データセットを保存していて、それが圧迫してたので消した。
『前処理大全』の3-3を読んだ。
極値、代表値の算出の話。
平均値
「平均値」の他の代表値とは違う利点として計算計算コストが低いから良いって言うのは、言われてみそうだがあまり意識したことなかった。
SQLの関数
SQLのMEDIAN関数やPERCENTILE_CONT関数を使ったことないな・・・
まあ、SQLでそもそもデータ分析をしたことない。手元に持ってきてPythonでやってる。そこまで巨大なデータを扱わないからか？
python
agg関数はlambdaが使えるのか。
あとは、パーセンタイルの計算としてnumpy.percentileがある。
gist.github.com
バージョン1.15からnumpy.quantileも使える。（今は1.14）
percentileは引数が1~100の整数だが、quantileは0~1の少数で指定できるようだ。
numpy.quantile — NumPy v1.15.dev0 Manual
歌詞取得スクリプト
「歌詞タイム」ってサイトで歌詞をちょっとコピーしようとしたら出来なかったので取得スクリプトを作った。
DevToolを開いても取得できたけど
タグなどが挟まっていたので「もうスクリプト作るか・・・」ってなった。
最初、HTML要素を見て歌詞が見当たらなかったので歌詞は別APIで取得してるのかと思ったけどただエンコードしてjavascriptで持ってるだけだった。
なのでHTMLをテキストで取得してデコードした後に正規表現で対象箇所を抽出するだけで歌詞を取得できた。（
HTTP HEADERなどを駆使してもっとマシな作りに出来ないか・・・）
なんかの法律に抵触すると面倒臭いのでスクリプトは公開しないようにします。
コンソールに求人を埋め込むAbemaTV
AbemaTVのページでコンソールを開くと求人が出てくるの、割と面白くないですか? pic.twitter.com/BjLqlUnILm
— keymoon (@kymn_) 2018年5月3日
ワロタwwww
ちなみに要素はアイコンが。
要素はもっとワロタ pic.twitter.com/JcvNIH9OPJ
— tdual@情報幾何 (@tdualdir) 2018年5月3日
日記20180502
まとめ
- 昨日眠れなかったので一日中頭があまり働かない感じだった。
- 『Pythonで学ぶ新しい統計学の教科書』1部6章まで読んだ。
- 『前処理大全』3章「集約」の3−1まで読んだ。
- はてなにJupyter notebookを貼る方法
- アイロン買った。
- 畳み込み層の処理は相互相関と言う記事をポストした。
- RAKERUに行った
- 明日の予定
昨日眠れなかったので一日中頭があまり働かない感じだった。
昨日飲んだMONSTER ENERGYのせいだと思う。ああ言うのを飲むとよく眠れなくなる。（が集中力はやはり上がるのでつい飲んでしまう。）
『Pythonで学ぶ新しい統計学の教科書』1部6章まで読んだ。
まだ、
- 母集団分布の推定 = 分布の決定（正規分布や二項分布など） + パラメータの推定
- 標本の統計量=母集団分布のパラメータの推定値
みたいな統計学の基礎の話
「母数」は分布のパラメータのことでサンプルサイズじゃないってことを強調して書いてた。
個人的にはもう「母数」って言葉も使うのやめてパラメータで良いんじゃないかと思う。
『前処理大全』3章「集約」の3−1まで読んだ。
groupbyした後に
aggで集計処理をまとめてできる。
gist.github.com
"May the agg be with groupby(aggがgroupbyとともにあらんことを)"
はてなにJupyter notebookを貼る方法
gistにipynbファイルをD&Dで追加して、そのリンクを貼るだけだった。
気をつけるのはtextareaに入れるとちゃんと貼れないのでこの状態の時にDropする。
アイロン買った。
最近買ったシャツのシワが目立つので、人生初のアイロンを買った。
t-falの変な形のスチームアイロン。アイロン台を買いたくなかったのでハンガーとかに掛けたままシワが伸ばせるやつにした。
まだ使ってない。
畳み込み層の処理は相互相関と言う記事をポストした。
tdual.hatenablog.com
割と反響があったみたい。
Toeplitz matrixの解説やnp.correlationでcnnを実装する話もいつか書きたい。
RAKERUに行った
オムライス専門のチェーン店。初めて行った。吉祥寺のヨドバシカメラ横のガストの下。
オムライスだけでなくパンも美味かった。後、店内がカラフルで子供受けしそうな内装だった。
明日の予定
pycon miniのスライドをいい加減完成させたい。
深層学習の畳み込み層の処理は「畳み込み」じゃなかった件
畳み込み層の処理は厳密には畳み込みではなかったのか・・・
畳み込み層
数学で言うところの畳み込み
畳み込みの定義
２次元の畳み込みは次のように定義されています。
離散値の場合は行列で表現できて
となります。
畳み込み層の処理と比較
(1)と(3)は非常に似てますが、(1)では、(3)では
となっていて、添え字が違うことに気付きます。
実は(1)の処理は相互相関(cross-correlation)と呼ばれていて畳み込みとは別物です。
相互相関
相互相関と畳み込みの相違点
まずは相互相関と畳み込みの違いを見ます。
以後、簡単のために１次元実数値連続関数を扱います。
１次元実数値連続関数の畳み込み
１次元実数値連続関数の相互相関
((5)の関数が実数値でない場合はの代わりに複素共役
を取る必要がある。）
交換律
畳み込みの重要な性質の一つとして交換律があります。
ここでとすると
\begin{align}
(f*g)(a) &= \int_{-\infty}^{\infty} g(a-x)f(x)dx\\
&= \int_{\infty}^{-\infty} g(x')f(a-x')d(a-x')\\
&= \int_{\infty}^{-\infty} g(x')f(a-x')d(-x')\\
&=\int_{-\infty}^{\infty} f(a-x')g(x')dx' = (g*f)(a)
\end{align}
つまり、
と言う交換律が成り立つことがわかります。
しかし、(5)の相互相関の形では交換律は満たしません。
フーリエ変換
畳み込みのフーリエ変換がそれぞれの関数のフーリエ変換の積になってるは有名な定理で、工学などでもよく応用されます。
\begin{align}
F[f*g] &= \int_{-\infty}^{\infty}\left(\int_{-\infty}^{\infty} g(a-x)f(x)dx\right)e^{-ika}da\\
&= \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} g(a')f(x)e^{-ik(a'+x)}dxda'\\
&=\int_{-\infty}^{\infty}f(x)e^{-ikx}dx \int_{-\infty}^{\infty}f(a')e^{-ika'}da' \\
&= F[f] F[g]
\end{align}
式変形の二行目ではとした。
これを相互相関でやると
\begin{align}
F[f\star g] &= \int_{-\infty}^{\infty}\left(\int_{-\infty}^{\infty} g(a+x)f(x)dx\right)e^{-ika}da\\
&= \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} g(a')f(x)e^{-ik(a'-x)}dxda'\\
&=\int_{-\infty}^{\infty}f(x)e^{+ikx}dx \int_{-\infty}^{\infty}f(a')e^{-ika'}da' \\
&= F^{*}[f] F[g]
\end{align}
ここではとした。
フーリエ変換の複素共役を取ったものとフーリエ変換に積になっています。（クロススペクトルと言う。）
非常に似た形をしてますが性質の相違点はいくつかあるようです。
相互相関の意味
相互相関の意味ですが、座標aでの二つの関数の類似度を表しています。似てるほど相互相関は大きくなります。
例えば、ノイズに混ざった二つの信号があった場合に相互相関を使ってそれが同じ信号かどうか判断できます。
（参考: https://qiita.com/inoory/items/3ea2d447f6f1e8c40ffa）
また、相互相関の性質を利用してフィルターなのども出来ます。
import numpy as np import matplotlib.pyplot as plt sig = np.repeat([0., 1., 1., 0., 1., 0., 0., 1.], 128) sig_noise = sig + np.random.randn(len(sig)) corr = np.correlate(sig_noise, np.ones(128), mode='same') / 128 clock = np.arange(64, len(sig), 128) fig, (ax_orig, ax_noise, ax_corr) = plt.subplots(3, 1, sharex=True) ax_orig.plot(sig) ax_orig.plot(clock, sig[clock], 'ro') ax_orig.set_title('Original signal') ax_noise.plot(sig_noise) ax_noise.set_title('Signal with noise') ax_corr.plot(corr) ax_corr.plot(clock, corr[clock], 'ro') ax_corr.axhline(0.5, ls=':') ax_corr.set_title('Cross-correlated with rectangular pulse') ax_orig.margins(0, 0.1) fig.tight_layout() fig.show()
出力
畳み込み層に「畳み込み」を適応するとどうなるのか？
畳み込み層に「畳み込み」を適応するとどうなるのか？
結論から言うと何も問題ないです。ニューラルネットワークの層で相互相関と畳み込みの相違点は関係ありません。
ただカーネルの上下左右を反転して学習してるだけです。（このカーネルをフリップカーネルと言います。）
図2の様にインプットのラベルを−1から始めてカーネルを反転すると畳み込みの式と一致してることがわかると思います。
終わりに
な ぜ Convolutional と 名 付 け た し。
まあ、ただの名前だし、数学の用語が誤用されることはよくあることなのであまり気にしても仕方ないね。
notebookはここに置いときますね。
github.com
ツイッターやっているのでフォローお願いします！
↓今すぐフォローすべきキラキラ アカウント
じゃあの。
参考文献
日記20180501
なんとなく日記書いてみるか。（いつまで続くかわからないけど）
- Finderのメニューバーにデスクトップを置いた
- 「PYTHONで学ぶあたらしい統計学の教科書」を買った
- 「PYTHONで学ぶあたらしい統計学の教科書」を１部４章まで読んだ
- 畳み込み層
- Python2.7.15は最後（last）のバージョンじゃなく最新(latest)のバージョンらしいな・・・
Finderのメニューバーにデスクトップを置いた
基本的にMacBookのデスクトップにファイルを置かないんだけど、友達に少し貸してたらなんか色々ファイルを置かれていた。
デスクトップのファイルを確認するためにFinderからユーザーディレクトリを開いてもデスクトップフォルダが見つからない・・・
ターミナルに行ってlsコマンドで確認するとちゃんとある・・・なんで隠してるんだよ・・・
仕方ないのでopenコマンドで開いてファイル確認したけど、これから毎回それするの面倒臭いのでFinderのメニューにデスクトップを置くことにした。
まずはFinderから隠してるフォルダを表示しなければと思い、隠してるフォルダの表示の仕方を調べた。
ターミナルでやる方法が見つかった。
defaults write com.apple.finder AppleShowAllFiles TRUE killAll Finder
これで全てのファイルがFinderから確認できるようになった。デスクトップフォルダも表示されたのでドラッグ&ドロップでメニューに追加した。
全てのファイルが表示されたままではみづらいので
defaults write com.apple.finder AppleShowAllFiles FALSE killAll Finder
で元に戻した。
ってか、なんでデフォルトでデスクトップフォルダ隠すんだよ・・・（それとも過去の自分が何かやったのか？ゴールデンウィークが終わって出社したら会社のMacBookを確認してみるか。）
「PYTHONで学ぶあたらしい統計学の教科書」を買った
「PYTHONで学ぶあたらしい統計学の教科書」のKindle版を買った。
会社のチャットで面白そうと言う発言があったのでちょっと興味を持っていた。
紙の本の方が先に出てて本屋で少し立ち読みして良さそうなのでKindle版が出たら買おうと思ってた。
「前処理大全」をまだ読み終わってないので並行で読んでいくか。
「PYTHONで学ぶあたらしい統計学の教科書」を１部４章まで読んだ
この本は７部構成でその中の章も沢山あるので量にビビってたけど細かく切ってるだけみたいだ。１ページに３章が詰まってるとかもざらにある。
母集団と標本を池の魚を釣る例えてその後もその例えがちょいちょい出てくる。
１部は統計学の基礎で、ここまでは用語の説明だけだった。
畳み込み層
Pycon miniの資料を作ってる途中でCNNのことを調べてたら畳み込み層の処理は厳密には畳み込みでないと言うのをGoodfellowらの『深層学習』の中で見つけた。
ひょっとしたら研究者なら当然知っているような内容かもしれないけど自分なりにまとめて別の記事として書く。
Python2.7.15は最後（last）のバージョンじゃなく最新(latest)のバージョンらしいな・・・
Twitterで見た面白いポスト
I misread the original post -- it's the *latest*, not the *last*. Maybe I was *hoping* it's the last, but the intention of the team is not to stop supporting it now. Apologies!
— Guido van Rossum (@gvanrossum) 2018年5月1日
ワロタwwwww
MacBook Proにアプリインストール
日記残すか。
MacBook Proにアプリインストール
朝にこの記事を読んで「Alfred」と「CLCL lite」と「RecordIt」をインストールした。
www.danshihack.com
Alfredは昔は使ったけどこのMacBook Proにしてからインストールしてなかったから再び入れてみた。
まあ、あまり使ってなかったので入れてなかったワケで、これからも使わないかもしれないけど定番アプリなのでとりあえずは入れてみた感じですね。
CLCL liteは結構いいな。
commandやoptionを連打することでアプリを開くことが出来る。これくらいシンプルな方が好きかも。
commandにTrelloのアプリを割り当ててすぐにメモが出来るようにした。
optionにはChromeを割り当ててすぐにブラウザに戻れるようにした。
もう一つあれば、iTermに割り当てたのにぐぬぬ・・・
CLCLは多用しそうですね。
RecordItも使える！
Recordit: Record screencasts fast & free! with GIF Support!
画面を撮影してGIFアニメーションにしてくれるアプリ。
今までLICEcapを使っていたけどRedordItの方が使いやすい。
上のメニューバーにアイコンが出て来てそこで範囲選択や録画/停止が出来る。
また、アニメーションはローカルに保存されるわけではなくrecordit.coにアップロードされているのでURLを貼るだけで共有できる。
アニメーションにはメニューバーから選択できてブラウザで開かれる。（パスワードで鍵も掛けられるみたい。）
デフォルトではmp4形式で再生されてる。右下のGIFボタンを押すことでGIFファイルが開かれる。
ローカルに保存したければ、一度ブラウザで開いてから右クリックで保存。（サイトがいつ無くなるかわからないのでローカルに保存した方がいいか・・・？）
これはめっちゃ良い。
A/Bテストの評価をベイズ統計でやってみない？
遊びでA/Bテストの評価をベイズ統計でやってみたら、思いのほか面白かったので記事に残します。
用語の定義
コンバージョン
コンバージョンとは「Webサイト上で起きた最終的な成果」のことです。
具体的に何を意味するかはサイトの種類によっては様々です。
例えば、ECサイトでは商品の購入で、SNSでは会員登録などです。
コンバージョン率
コンバージョン率は「成果に繋がる最初の行動に対して実際に成果に繋がった割合」のことです。
ECサイトではある商品が購入された数をその商品が表示された数で割ったものがコンバージョン率となります。
（例えば、ある商品の表示数が100でその内で購入数が2ならばコンバージョン率は 2/100=0.02になります。）
A/Bテスト
A/Bテストとは複数の施策に対して良否をつけるためのテストです。
例えば、ウェブサイトのユーザビリティーを上げるためにデザインAとデザインBを同じ期間試してどちらが表示数が多いのか比べたり、
広告のパターンAとパターンBを出し分けてどちらのコンバージョン率が大きいか比べたりすることです。
今回は広告のA/Bテストを考えます。
広告の表示に対して購入が何件有ったのか（コンバージョン率）をA/Bテストします。
コンバージョンの確率分布
表示に対して購入が有ったか、無かったのかの二通りです。
それはBernoulli（ベルヌーイ）分布で表現できます。
購入が有った時x=1
購入が無かった時x=0
とすれば、がx=1になる時の確率なので
がコンバージョン率ということになります。
xは確率変数,はパラメータ（母数とも呼ばれる）と呼ばれています。
なぜベイズ統計を使うのか
割合の問題点
コンバージョン率は購入数÷表示数であることを説明しました。ここでパターンAのコンバージョン率を確かめる実験を二回行ったとします。
実験①: 表示数100, 購入数2 でした。コンバージョン率は2/100=0.02になります。
実験②: 表示数10000、購入数200でした。コンバージョン率は200/10000=0.02になります。
例えば実験② は実験①に比べて長い間表示してたので表示数も多かったとします。
どちらの方が信頼できますか？
実験①はデータが少なくてたまたまコンバージョン率が0.02になった様に思えますが、実験②は実験①よりは信頼できそうです。
しかし、両方とも同じコンバージョン率になります。
コンバージョン率の様な割合はこの信頼性を表現出来ていない事になります。
尤度と最尤法
ベイズ統計の前に最尤法という方法でパラメータの推定をしてみます。
尤度
確率分布はパラメータを固定した場合の確率変数の関数でしたが、同じ関数を確率変数を固定した場合のパラメータの関数と捉えた物を尤度と言います。
この場合の確率変数には観測データが入ります。観測データの表示数がNの場合に尤度は以下の様になります。
最尤法
尤度が最大値の時のパラメータが尤も（もっとも）適した値であるという方針のもとでパラメータを推定する方法です。
最大値の時の引数の値を求めるには微分が0のところを求めればいいのですが、計算を楽にするために尤度の対数を微分をします。（対数関数は単調増加関数なので引数が最大のときに対数関数の値も最大になるので問題ない。）
\begin{align*}
\frac{d\log L(\theta)}{d\theta}
&=\frac{d}{d\theta}\log\prod_{i=1}^N\theta^{x_i}(1-\theta)^{1-x_i}\\
&=\frac{d}{d\theta} \sum_{i=1}^N \left(x_i \log\theta +(1-x_i)\log(1-\theta)\right) \\
&=\sum_{i=1}^N \left(\frac{x_i}{\theta} - \frac{1-x_i}{1-\theta}\right) =0
\end{align*}
\begin{align*}
\sum_{i=1}^N \left(x_i(1-\theta) - (1-x_i)\theta \right) &=0\\
\sum_{i=1}^N(x_i -\theta) &=0\\
\sum_{i=1}^N x_i - N\theta &= 0\\
\theta = \frac{\sum_{i=1}^N x_i}{N}
\end{align*}
Nは表示数で、はコンバージョンが有ったときに1で無かったときは0なので、
の尤も適した値はコンバージョン率だとわかります。
ですので、最尤法を使っても割合の問題が出てきて信頼性を測るすべはありません。
ベイズ統計
ベイズ統計ではパラメータを１つの決まった値ではなくて確率分布として考えます。
こうする事で割合や最尤法では測れなかった信頼性が測れる事が後にわかります。
ベイズの定理
確率分布のベイズの定理は以下の様になります。
は事後分布と呼ばれます。これをパラメータ(コンバージョン率)の分布と捉えます。
は
に依存してないので今回はあまり気にする必要がありません。
は尤度でベルヌーイ分布と同じ形でしたね。
は事前分布と呼ばれています。
共役事前分布
事前分布の形を決定しないと事後分布の計算はできませんが形は謎です。
ですので、この分布は適当に決めてしまいましょう。
ただ、事前分布を適当に決めるにしても何らかの指針は欲しいものです。また適当に決めると計算できない場合がほとんどです。
そこで共役事前分布というものを考えます。この共役事前分布は事後分布が計算できるように選んだ事前分布デス。また、事前分布と事後分布は同じ形の確率分布になります。
どういうものがあるのかはwikiなどを見て欲しいのですが、尤度がベルヌーイ分布の場合は事前分布としてベータ分布を選べば良いことがわかっています。
ベータ分布
はベータ関数ですが、積分をして1とするための規格化係数なのであまり気にしなくていいです。
形は図1の様になります。
に依存して形が大きく変わっていることがわかります。
また、が共に１の時は一様分布になっているのもわかります。
事後分布の導出
共役事前分布のベータ分布を使って事後分布を計算して見ましょう。（事後分布もベータ分布になるはずですが。）
\begin{align*}
f(\theta\mid x) &\propto f(x\mid\theta)p(\theta)\\
&\propto \left(\prod_{i=1}^N\theta^{x_i}(1-\theta)^{1-x_i}\right)\theta^{\alpha -1}(1-\theta)^{\beta-1}\\
&=\theta^{\sum^N_{i=1}x_i+\alpha-1}(1-\theta)^{\sum^N_{i=1}(1-x_i)+\beta-1}\\
&= \theta^{\alpha'-1}(1-\theta)^{\beta'-1}
\end{align*}
パラメータがのベータ分布になりました。
但し、
です。はコンバージョン数だけ増加して、
はコンバージョンされなかった数だけ増加しています。
事前分布のベータ分布のパラメータとして（一様分布になる）を選ぶことにしましょう。
事後分布のグラフ
必要なライブラリをpipで入れます。（pandasは描画には使いませんが後で使うので一緒にインストールします。）
pip install numpy scipy matplotlib pandas
事後分布を描画する関数を作ります。
import numpy as np import scipy.stats from matplotlib import pyplot as plt def show_beta(h, t): x = np.linspace(0, 1, 1000) plt.xlim(0, 1) plt.ylim(0, 5) plt.xlabel(r"$\theta$",fontsize=15) plt.ylabel(r"$f(\theta\mid x)$",fontsize=15) alpha = 1 beta = 1 alpha += h beta += t f = scipy.stats.beta(alpha, beta) m = f.mean() v = f.var() y = f.pdf(x) plt.title(r'$\alpha$=%d+1, $\beta$=%d+1 [mean=%.2f, var=%.7f]' % (h,t,m,v), fontsize=15) plt.plot(x, y) plt.show()
hにコンバージョン数、tに非コンバージョン数を入れると事後分布を描画する様にしました。
また、図のタイトルに平均(mean)と分散(var)も表示しています。
実験①：表示数100,購入数2の時を描画してみましょう。
show_beta(2,100-2)
出力
実験②：表示数10000,購入数200の時は
show_beta(200,10000-200)
出力
図2と図３を見比べると、図２の分布には幅があり不確定であることがわかります。一方で図３は幅がほとんどなく値がほぼ確定してることがわかります。この幅は分散に依存してるので、実際に分散（var）を見ると図３の方が全然小さいことがわかります。この様にパラメータ（コンバージョン率）を分布で表現すれば、その値の信頼性が分散という形で確認できることがわかりました。
これがベイズ統計を使う利点の一つです。
ベイジアンA/Bテストの実装
準備ができたので、ベイズ統計を使ってA/Bテストの結果を評価しましょう。
これはベイジアンA/Bテストと呼ばれています。
早速実装しまし。
コード
class BayesianAB: def __init__(self, params=None, size=2): """ params are a list of Beta distribution's initial params. e.g. [[alpha of A, beta of A], [alpha of B, beta of B]] """ if params is None: self.params = [] for _ in range(size): self.params.append([1,1]) self.size = len(self.params) self.data = [] for _ in range(size): self.data.append([0,0]) print("the number of comparison: ", self.size) self.sampling() def update(self, data, sampling=True): """ data are a list of pairs of impression and conversion. e.g. [[imp of A, conv of A], [imp of B, conv of B]] """ if self.size != len(data): print("No match of the size.") for p, current, new in zip(self.params, self.data, data): imp = new[0] conv = new[1] current[0] += imp current[1] += conv p[0] += conv p[1] += (imp - conv) if sampling: self.sampling() def mean_ver(self): """ return [(mean of A, variance of A), (mean of B, variance of B), ...] """ return [(posterior.mean(), posterior.var())for posterior in self.posterior_list] def sampling(self, n_samples=50000): print("num of samples: ", n_samples) self.posterior_list = [beta(*p).rvs(n_samples) for p in self.params] def show_beta(self, title="", save=False, labels=None): plt.figure(figsize=(10, 5)) plt.title("Posterior distribution "+ title) cmap = plt.get_cmap('jet') color_list= [] for i, posterior in enumerate(self.posterior_list): color =cmap(0.25*(i+1)) color_list.append(color) plt.hist(posterior, bins=100, histtype="stepfilled", normed=True, color=color, alpha=0.5) handles = [Rectangle((0,0),1,1,color=c, ec="k", alpha=0.5) for c in color_list] if labels is None: labels = [chr(65+i) for i in range(self.size)] # create A,B,... plt.legend(handles, labels) if save: plt.savefig("{}.png".format(title)) plt.show() def diff_prob(self, index_high, index_low): prob = (self.posterior_list[index_low] < self.posterior_list[index_high]).mean() if prob < 0.5: prob = 1 - prob return prob def show_metrics(self): print(self.data) def metrics(self, labels=None): if labels is None: labels = [chr(65+i) for i in range(self.size)] # create A,B,... return pd.DataFrame(self.data, index=labels, columns=["Impressions", "Conversions"])
使用例
初期化するときに事前分布のベータ分布のパラメータを設定できる様にしましたが、ほとんどの場合は一様分布で十分なのでそのまま初期化しましょう。
またsizeで比較する数を指定できます。
abtest = BayesianAB(size=2)
出力
the number of comparison: 2 num of samples: 50000
作成したクラスの使い方を説明するのですが、ここではダミーの行動データを作りましょう。
def dummy_genrator(conversion_rate, n_impression): return np.array(bernoulli(conversion_rate).rvs(n_impression)) a = dummy_genrator(0.015, 204) b = dummy_genrator(0.021, 200) imp_a = len(a) imp_b = len(b) conv_a = np.count_nonzero(a) conv_b = np.count_nonzero(b) print("A \n表示数: {}\n購入数: {}\n".format(imp_a, conv_a)) print("B \n表示数: {}\n購入数: {}\n".format(imp_b, conv_b))
出力
A 表示数: 204 購入数: 2 B 表示数: 200 購入数: 3
このデータを使って更新します。
abtest.update([[imp_a, conv_a], [imp_b, conv_b]])
これで更新されているので、平均と分散を確認すると、
abtest.mean_ver()
出力
[(0.014610260919694517, 7.003148406424335e-05), (0.019850764074203515, 9.568870563433339e-05)]
[(平均,分散),(平均,分散)]という形式で出力してるので
パターンAの分布の平均が0.014610260919694517, 分散が7.003148406424335e-05
パターンBの分布の平均が0.019850764074203515, 分散が9.568870563433339e-05
です。
図も描画して見ると
abtest.show_beta()
出力
updateの時の順番と同じインデックスを指定してやります。Aなら0,Bなら１です。
abtest.diff_prob(0,1)
出力
0.6682600000000001
完全に一致してる場合は0.5になります。ここでは0.668...なのでやはりあまり有意差がないということです。
更新したデータを見るにはshow_metricsを使います。
abtest.show_metrics()
出力
[[204, 2], [200, 3]]
[[Aの表示数,Aの購入数],[Bの表示数],[Bの表示数]]という形式で出力してくれます。
本番っぽい使い方
大体の使い方はわかったと思うので、本番っぽい使い方をしてみましょう。
今回はパターンは2パターンではなく３パターンの比較します。それぞれのパターンをA,B,Cと呼びます。
ある程度データが溜まった想定でダミーデータの生成します。
a = dummy_genrator(0.015, 1002) b = dummy_genrator(0.021, 1200) c = dummy_genrator(0.016, 1130) imp_a = len(a) imp_b = len(b) imp_c = len(c) conv_a = np.count_nonzero(a) conv_b = np.count_nonzero(b) conv_c = np.count_nonzero(c) print("A \n表示数: {}\n購入数: {}\n".format(imp_a, conv_a)) print("B \n表示数: {}\n購入数: {}\n".format(imp_b, conv_b)) print("C \n表示数: {}\n購入数: {}\n".format(imp_c, conv_c))
出力
A 表示数: 1002 購入数: 22 B 表示数: 1200 購入数: 22 C 表示数: 1130 購入数: 18
早速グラフをみてみましょう。
abtest = BayesianAB(size=3) #今回は3パターンなのでsizeは３ abtest.update([[imp_a, conv_a], [imp_b, conv_b],[imp_c, conv_c]]) abtest.show_beta()
出力
3つとも重なっている部分が大きいですね。
数値でも確認しましょう。
print(abtest.diff_prob(1,0)) # BとAの比較 print(abtest.diff_prob(1,2)) #BとCの比較 print(abtest.diff_prob(2,0)) #CとAの比較
出力
0.72712 0.66822 0.84508
有意差がわかりませんね。もっとデータが必要です。データが溜まるまで待ちましょう。
データが溜まった想定でダミーデータの生成します。
a = dummy_genrator(0.015, 25100) b = dummy_genrator(0.021, 22000) c = dummy_genrator(0.016, 23300) imp_a = len(a) imp_b = len(b) imp_c = len(c) conv_a = np.count_nonzero(a) conv_b = np.count_nonzero(b) conv_c = np.count_nonzero(c) print("A \n表示数: {}\n購入数: {}\n".format(imp_a, conv_a)) print("B \n表示数: {}\n購入数: {}\n".format(imp_b, conv_b)) print("C \n表示数: {}\n購入数: {}\n".format(imp_c, conv_c))
出力
A 表示数: 25100 購入数: 364 B 表示数: 22000 購入数: 469 C 表示数: 23300 購入数: 354
先ほど作ったBayesianABのインスタンスにデータの追加更新をしてやりましょう。
abtest.update([[imp_a, conv_a], [imp_b, conv_b],[imp_c, conv_c]])
abtest.show_beta()
print(abtest.mean_ver())
出力
[(0.014823136368344494, 5.586078465062482e-07), (0.0212059201285066, 8.839507467322899e-07), (0.01526854730312879, 6.223394152971891e-07)]
全ての分布が最初の時よりも幅が小さくなっていることがわかります。また、パターンBが明らかに離れているので有意差があることがわかります。
print(abtest.diff_prob(1,0)) # BとAの比較 print(abtest.diff_prob(1,2)) #BとCの比較 print(abtest.diff_prob(2,0)) #CとAの比較
出力
1.0 1.0 0.65888
A,CとBは100%違うことがわかります。そして、Bのコンバージョン率の平均0.021と一番高いのでこの広告を選ぶべきだとわかります。
しかし、一番効果がある広告を選ぶならこれでいいのですが、目的がA,B,Cのそれぞれの有意差を調べることならまだAとCの有意差ははっきりしていません。
ですので、もうちょっとデータを溜めましょう。
データが溜まった想定でダミーデータの生成します。
a = dummy_genrator(0.015, 251000) b = dummy_genrator(0.021, 220000) c = dummy_genrator(0.016, 233000) imp_a = len(a) imp_b = len(b) imp_c = len(c) conv_a = np.count_nonzero(a) conv_b = np.count_nonzero(b) conv_c = np.count_nonzero(c) print("A \n表示数: {}\n購入数: {}\n".format(imp_a, conv_a)) print("B \n表示数: {}\n購入数: {}\n".format(imp_b, conv_b)) print("C \n表示数: {}\n購入数: {}\n".format(imp_c, conv_c))
出力
A 表示数: 251000 購入数: 3808 B 表示数: 220000 購入数: 4544 C 表示数: 233000 購入数: 3811
追加更新しましょう
abtest.update([[imp_a, conv_a], [imp_b, conv_b],[imp_c, conv_c]])
abtest.show_beta()
print(abtest.mean_ver())
出力
[(0.015141034587728202, 5.332992221848848e-08), (0.0207078937798314, 8.243283715216886e-08), (0.016250559254538843, 6.199512187766625e-08)]
abtest.diff_prob(2,0)
出力
0.99964
Cの方が良いと言う事がわかりましたね。
このようにベイジアンA/Bテストでは予め有意差を測るためにデータ数がどの位必要なのかを知ってる必要はなくインクリメントに更新して評価することが出来きます。例えば、リアルタイムに更新をすれば有意差が出たと判断した時点ですぐにA/Bテストをやめて機会損失を最小限に抑える事ができます。
これもベイズ統計の強みです。
カイ二乗検定と比較
せっかくなのでカイ二乗検定と比較してみましょう。
from scipy.stats import chi2_contingency def chi2(imp_a, conv_a, imp_b, conv_b,title=""): print(title) metrics = [ [conv_a, imp_a - conv_a], [conv_b, imp_b - conv_b], ] print(metrics[0]) print(metrics[1]) x2, p, dof, expected = chi2_contingency(metrics) print("カイ二乗値: ", round(x2, 5)) print("p値: ", p) if p < 0.05: print("有意差あり。") else: print("有意差なし。") print("\n")
最初のデータのとき
imp_a_1 = 1002 conv_a_1 = 22 imp_b_1 = 1200 conv_b_1 = 22 imp_c_1 = 1130 conv_c_1 = 18 chi2(imp_a_1, conv_a_1, imp_b_1, conv_b_1,"###1nd A vs B ###") chi2(imp_a_1, conv_a_1, imp_c_1, conv_c_1,"###1nd A vs C ###") chi2(imp_c_1, conv_c_1, imp_b_1, conv_b_1,"###1nd B vs C ###")
出力
###1nd A vs B ### [22, 980] [22, 1178] カイ二乗値: 0.20435 p値: 0.6512354862810796 有意差なし。 ###1nd A vs C ### [22, 980] [18, 1112] カイ二乗値: 0.74604 p値: 0.38773258145798684 有意差なし。 ###1nd B vs C ### [18, 1112] [22, 1178] カイ二乗値: 0.08233 p値: 0.7741616454611604 有意差なし。
2回目にデータが溜まった時
imp_a_2 = 1002 + 25100 conv_a_2 = 22 + 364 imp_b_2 = 1200+ 22000 conv_b_2 = 22 + 469 imp_c_2 = 1130 + 23300 conv_c_2 = 18 + 354 chi2(imp_a_2, conv_a_2, imp_b_2,conv_b_2,"###2nd A vs B ###") chi2(imp_a_2, conv_a_2, imp_c_2,conv_c_2,"###2nd A vs C ###") chi2(imp_c_2, conv_c_2, imp_b_2,conv_b_2,"###2nd B vs C ###")
（ベイジアンA/Bテストでは追加更新してたので表示数と購入数は1回目のデータとの合計になる。）
出力
###2nd A vs B ### [386, 25716] [491, 22709] カイ二乗値: 28.2126 p値: 1.086949244640548e-07 有意差あり。 ###2nd A vs C ### [386, 25716] [372, 24058] カイ二乗値: 0.13625 p値: 0.7120340373703777 有意差なし。 ###2nd B vs C ### [372, 24058] [491, 22709] カイ二乗値: 23.24073 p値: 1.4293812107138311e-06 有意差あり。
A,BとCは有意差があって、AとCの有意差は分からない。ベイジアンA/Bテストと結論は同じ。
3回目のデータが溜まった時
imp_a_3 = 1002 + 25100 + 251000 conv_a_3 = 22 + 364 + 3803 imp_b_3 = 1200 + 22000 + 220000 conv_b_3 = 22 + 469 + 4544 imp_c_3 = 1130 + 23300 + 233000 conv_c_3 = 18 + 354 + 3811 chi2(imp_a_3, conv_a_3, imp_b_3,conv_b_3,"###3rd A vs B ###") chi2(imp_a_3, conv_a_3, imp_c_3,conv_c_3,"###3rd A vs C ###") chi2(imp_c_3, conv_c_3, imp_b_3,conv_b_3,"###3rd B vs C ###")
出力
###3rd A vs B ### [4189, 272913] [5035, 238165] カイ二乗値: 231.76377 p値: 2.4586705996516234e-52 有意差あり。 ###3rd A vs C ### [4189, 272913] [4183, 253247] カイ二乗値: 11.01697 p値: 0.0009028171800658043 有意差あり。 ###3rd B vs C ### [4183, 253247] [5035, 238165] カイ二乗値: 137.02089 p値: 1.1932289659422658e-31 有意差あり。
AとCに有意差がある事がわかりましたね。
有意差を調べると言う意味では結果は同じになりましたね。
