ロジスティック回帰を使って確率を予測したいときに「評価メトリクスとして使いたいのは AUC (areaUnderROC) じゃなくて Logarithmic loss (LogLoss) なんだよ！」と常々思っているのですが、現在の MLlib には二値分類 (BinaryClassificationEvaluator) 、多クラス分類 (MulticlassClassificationEvaluator)、回帰 (RegressionEvaluator) 用の Evaluator 実装しか用意されていなかったので、ついカッとなって実装してしまった次第です。
これくらいの機能は標準で取り揃えていてもいいんじゃないかと思うんだけどなあ…
XGBoost で構築した予測モデルを Java から利用したい、それも特徴ベクトルが一つ一つ、任意のタイミングで与えられるような オンライン環境下 で リアルタイムな予測 を実現するために利用したい、という目的を叶えるためのモジュールを作りました。
(XGBoost の凄さとか XGBoost そのものの使い方とか GBDT/GBRT の解説は本エントリにはありませんので、そのような情報を求めている方は他のブログエントリを読まれることをおすすめします。)
Java から XGBoost を利用しようとすると、XGBoostをJavaのwrapperを使用して実行する - TASK NOTES にあるように、DMLC が提供している xgboost4j を利用する手段がすでに存在しています。ただ、この xgboost4j をオンライン予測に適用する場合、下記に挙げるようないくつかの懸念があります。
DMatrix で表現するのが面倒mvn install する必要がある
そういうわけで、
XGBoost 互換な予測器を作って jCenter で公開しました。
README.md にもサンプルコード込みで使い方を書いていますので、合わせてご参照ください。
なお、タイトルでも「予測器」と明言しているとおり、学習 (training) 機能については割りきって一切実装をしていません (モデルの構築をオンラインですることは流石にないと思うので)。 そのため学習データを用意し、XGBoost を CLI で直接実行、もしくは Python や R のラッパーを経由するなどして別途モデルの構築を事前に済ませておく必要があります。
学習済みのモデルが用意できたら、そのモデルを new Predictor("/path/to/model-file") としてロードします。
予測をするには Predictor#predict() のメソッドを呼び出します。
予測の際の入力となる特徴ベクトルは、FVec インタフェースのオブジェクトとして表現する必要があります。
基本はお手持ちのデータの形式に合わせて FVec インタフェースを実装したクラスを用意していただくことになりますが、配列や Map でデータが表現されている場合には、以下のユーティリティメソッドを利用することができます。
FVec.Transformer#fromArray() で FVec オブジェクトに変換できますMap オブジェクト → FVec.Transformer#fromMap() で FVec オブジェクトに変換できますまた GBDT/GBRT のモデルを特徴ベクトルの変換に利用することを目的に、(分類・回帰の結果を出力するのではなく) GBDT/GBRT の各ツリーにおいて辿り着いたリーフのノード番号を出力する Predictor#predictLeaf() メソッドも用意されています。
元々このモジュールを作った動機の一つとして、オンライン利用に耐えうる速度性能で予測をしたい、という目的があったので xgboost4j と合わせてベンチマークをとってみました。 ベンチマーク計測の詳細は こちら のプログラムにあるとおりです。
| 機能 | xgboost-predictor | xgboost4j | 性能比 |
|---|---|---|---|
| モデルの読み込み | 49017.60 ops/s | 39669.36 ops/s | 1.24 |
| 単一の特徴ベクトルでの予測 | 6016955.46 ops/s | 1018.01 ops/s | 5910.48 |
| 複数の特徴ベクトルの予測 | 44985.71 ops/s | 5.04 ops/s | 8931.47 |
| リーフノードの出力 | 11115853.34 ops/s | 1076.54 ops/s | 10325.53 |
結果は上記のとおり、
となります。Pure Java にしただけで xgboost4j との性能差がこんなに出るものなのか… とちょっと不思議ではありますが、今回のベンチマークの計測に利用したテストデータは人工的に生成されたものなので、その影響があるのかもしれません。 そのため、実世界のデータを与えた場合にはまた違った結果になる可能性がありますのでご注意ください。
現時点の xgboost-predictor は、モデルとしては “gbtree” のみをサポートしています。また目的関数は “binary” および “multi” のみのサポートとなります。 (“gblinear” や他の目的関数のサポートについては、必要そうであれば対応する予定です。)
追記: 2016-09-27 最新のビルド手順は こちら に記載しています。
なんで Java から XGBoost を扱いたいのかはさておき、概ね XGBoostをJavaのwrapperを使用して実行する - TASK NOTES こちらのサイトの解説どおりではあるのですが、OS X で素直に java/create_wrap.sh を叩いてビルドしようとすると
$ ./create_wrap.sh
build java wrapper
clang-omp++ -Wall -O3 -msse2 -Wno-unknown-pragmas -funroll-loops -fopenmp -fPIC -fPIC -shared -o java/libxgboostjavawrapper.so java/xgboost4j_wrapper.cpp wrapper/xgboost_wrapper.cpp updater.o gbm.o io.o subtree/rabit/lib/librabit.a dmlc_simple.o -pthread -lm -I/Library/Java/JavaVirtualMachines/jdk1.8.0_66.jdk/Contents/Home/include -I/Library/Java/JavaVirtualMachines/jdk1.8.0_66.jdk/Contents/Home/include/linux -I./java
In file included from java/xgboost4j_wrapper.cpp:15:
/Library/Java/JavaVirtualMachines/jdk1.8.0_66.jdk/Contents/Home/include/jni.h:45:10: fatal error: 'jni_md.h' file not found
#include "jni_md.h"
^
みたいなエラーが出てしまうはずです。
この問題の対処法は実に簡単で、xgboost リポジトリのトップディレクトリ配下にある Makefile の 8行目を
before: export JAVAINCFLAGS = -I${JAVA_HOME}/include -I${JAVA_HOME}/include/linux -I./java
after : export JAVAINCFLAGS = -I${JAVA_HOME}/include -I${JAVA_HOME}/include/darwin -I./java
と書き換えるだけ、です。
続いて、java/xgboost4j ディレクトリ配下で mvn package コマンドを実行して jar ファイルを生成しようとしてみると、今度は
-------------------------------------------------------
T E S T S
-------------------------------------------------------
Running org.dmlc.xgboost4j.BoosterTest
Oct 29, 2015 2:51:55 AM org.dmlc.xgboost4j.DMatrix <clinit>
SEVERE: load native library failed.
Oct 29, 2015 2:51:55 AM org.dmlc.xgboost4j.DMatrix <clinit>
SEVERE: java.io.FileNotFoundException: File /lib/libxgboostjavawrapper.dylib was not found inside JAR.
Tests run: 1, Failures: 0, Errors: 1, Skipped: 0, Time elapsed: 0.121 sec <<< FAILURE!
testBoosterBasic(org.dmlc.xgboost4j.BoosterTest) Time elapsed: 0.063 sec <<< ERROR!
java.lang.UnsatisfiedLinkError: org.dmlc.xgboost4j.wrapper.XgboostJNI.XGDMatrixCreateFromFile(Ljava/lang/String;I[J)I
at org.dmlc.xgboost4j.wrapper.XgboostJNI.XGDMatrixCreateFromFile(Native Method)
at org.dmlc.xgboost4j.DMatrix.<init>(DMatrix.java:62)
at org.dmlc.xgboost4j.BoosterTest.testBoosterBasic(BoosterTest.java:75)
at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
...
と、test フェーズでこけてしまうかと思います。これは native ライブラリの拡張子が .so となってしまっているのが原因なので、これは create_wrap.sh の 12 行目 の移動先ファイル名を明示的に libxgboostjavawrapper.dylib と指定するか、もしくはすでに移動先にある native ライブラリの拡張子を .so から .dylib に書き換えてやれば OK です。
Approximate membership query (AMQ) が実現できるデータ構造としてそれなりに広く使われていそうな Bloom filter ですが、機能性を高めたバリエーションはそこそこ存在する一方で、空間効率を追求した・コンパクトなバリエーションはあんまり見つからないものです。
ここ最近、ふとデータ構造熱が高まってきたこともあったので、オリジナルの Bloom filter よりもコンパクトに表現できる代替データ構造を探してメモしてみました。
空間効率を追求するといっても、転送・永続化している状態での空間効率なのか、それともルックアップ可能な状態においての空間効率なのかで全然違うわけですが、この Compressed Bloom Filters は前者の転送時・永続化時の空間効率を選択したバリエーションです。
端的に言えば、Bloom filter のビット列における 0 と 1 の出現確率はランダム (1/2 = 0.5) というわけではなく偏りが生じることがあるため、算術符号で符号化すればよりコンパクトになるじゃん、というものです。 ただ、この圧縮はあくまでも転送・永続化時のことしか考えていないようで、Bloom filter が素早くルックアップできる状態にあるためには結局復号した状態で保持する以外に他ないようです。
Cache-, Hash- and Space-Efficient Bloom Filters という論文で、Space-Efficient な手法として提案されている approximate membership query を提供するデータ構造です (Bloom filter と同じ機能を提供しているが、厳密には Bloom fitler ではない)。
要素のハッシュ値を算出するところまでは Bloom filter と同じですが、そのハッシュ値を
することで、Bloom filter よりもコンパクトな表現を実現しています (ハッシュ値は一様分布しているものとして、そのハッシュ値の gap をとると幾何分布になり、幾何分布に従う数値を符号化するには Golomb 符号が都合いい、ということ)。
この方式でのルックアップ操作は、あるハッシュ値が Golomb-Compressed Sequence に含まれるか否かだけを探索により判定すればいいので、Bloom filter のようなビット列への複数回のランダムアクセスは生じません。 とは言えど、この圧縮表現のままでは効率的なランダムアクセスが実現できないため、
ことで、復号にかかる時間効率を高められるようにしています。 $I$ を大きくすれば大きくするほど、空間効率がよくなるが探索効率は悪くなり、 $I$ を小さくすると空間効率は悪くなるが探索効率が上がる、というトレードオフの関係になっています。
Cuckoo Filter: Practically Better Than Bloom では、Cuckoo hashing を利用して、approximate membership query を提供するデータ構造を提案しています。
Cuckoo hashing の詳細説明は kumagi 先生の資料 に譲るとして、Cuckoo filter では要素から $f$ ビットのフィンガープリント (ハッシュ値の一つだと考えれば OK) と、2 つのハッシュ値を算出し、 Cuckoo hashing に利用します。より具体的には、2 つのハッシュ値を配列で表現されるバケット (複数のフィンガープリントを記録できる入れ物) のインデックス決定に利用し、フィンガープリントをそのバケットに追記する、という使い方になります。
この Cuckoo filter で使われる 2 つのハッシュ関数はちょっとだけ特殊で、$h_1(x) = hash(x), h_2(x) = h_1(x) \oplus hash(x's fingerprint)$ となっています。つまりは、フィンガープリントと一方のハッシュ値がわかっていれば $h_1(x)$ と $h_2(x)$ を算出することができるわけです。
パラメータとしてはフィンガープリントを表現するビット数や、バケットの大きさ (異なるフィンガープリントを格納できる個数)、バケット数と、ちょっと多めなのが気になりますが空間効率も参照性能も Bloom filter より良さそうです。
さて今回はこのぐらいにして、Golomb-Cmpressed Sequence あたりを実装してみようかな…
最適化問題、とりわけ線形計画問題が気になり始めるお年頃になってきたので、 SCIP というソルバーを使ってみましたよ、というメモ。
まずは手元にある MBP に、 SCIPをMacにインストール - Qiita の手順を参考にインストールを試してみる。ところが、
src/rational.h:32:10: fatal error: 'gmp.h' file not found
などの gmp.h がない旨のコンパイルエラーが発生してしまう状況に遭遇してインストールできない。
なので、ひとまず make GMP=false としてお茶を濁すことしてみる (GMP を要求するのは ZIMPL らしく、かつその機能を利用することは今のところはないので)。
ビルド後の scip のバイナリ scip-X.X.X/bin/scip を立ち上げるとインタラクティブシェルが立ち上がる。
ひとまず適当な問題をソルバーで解かせてみようと思い、講義テキストらしきもの (PDF) を参考に、最適化問題を記述した LP ファイル (sample.lp) を作ってみる。
maximize
400 x1 + 300 x2
subject to
60 x1 + 40 x2 <= 3800
20 x1 + 30 x2 <= 2100
20 x1 + 10 x2 <= 1200
end
このファイルを scip のインタラクティブシェル上で読み込んでみる。読み込みには read コマンドを利用する。
SCIP> read /path/to/sample.lp
read problem </path/to/sample.lp>
============
original problem has 2 variables (0 bin, 0 int, 0 impl, 2 cont) and 3 constraints
2 つの変数、3 つの制約条件、と出力がでて、ちゃんと読み込めたようだ。
読み込んだ問題を実際に解いてみよう。問題を解くには optimize コマンドを利用する。
SCIP> optimize
feasible solution found by trivial heuristic after 0.0 seconds, objective value 0.000000e+00
presolving:
(round 1, fast) 0 del vars, 0 del conss, 0 add conss, 4 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs
presolving (2 rounds: 2 fast, 1 medium, 1 exhaustive):
0 deleted vars, 0 deleted constraints, 0 added constraints, 4 tightened bounds, 0 added holes, 0 changed sides, 0 changed coefficients
0 implications, 0 cliques
presolved problem has 2 variables (0 bin, 0 int, 0 impl, 2 cont) and 3 constraints
3 constraints of type <linear>
Presolving Time: 0.00
transformed 1/1 original solutions to the transformed problem space
time | node | left |LP iter|LP it/n| mem |mdpt |frac |vars |cons |cols |rows |cuts |confs|strbr| dualbound | primalbound | gap
* 0.0s| 1 | 0 | 3 | - | 194k| 0 | - | 2 | 3 | 2 | 3 | 0 | 0 | 0 | 2.700000e+04 | 2.700000e+04 | 0.00%
0.0s| 1 | 0 | 3 | - | 194k| 0 | - | 2 | 3 | 2 | 3 | 0 | 0 | 0 | 2.700000e+04 | 2.700000e+04 | 0.00%
SCIP Status : problem is solved [optimal solution found]
Solving Time (sec) : 0.00
Solving Nodes : 1
Primal Bound : +2.70000000000000e+04 (2 solutions)
Dual Bound : +2.70000000000000e+04
Gap : 0.00 %
表示内容が豊富でちょっと面食らってしまうけど、 problem is solved [optimal solution found] がポイント。最適解が求まったぽい。
optimize コマンドで解が求まったら、 display のコマンドでその解を表示してみよう。
SCIP> display solution
objective value: 27000
x1 30 (obj:400)
x2 50 (obj:300)
ちゃんと最適解が求まっているのがわかる。
optimize コマンドで求まった解を他の用途で使いたい、ということは結構よくあるはず。ファイルに書き出すことができればスクリプト言語でパースして利用する事もできるだろう。
それでは write コマンドでファイルに書き出してみよう。
SCIP> write solution result.sol
written solution information to file <result.sol>
出力された result.sol を less などで覗いてみれば、display solution でコンソールに表示した内容と同等のものが出力されているのがわかるはず。
SCIP のひととおりの使い方がだいたい分かってきたところで、他の問題も SCIP で解けるように、問題の記述方法を勉強してみよう。
まずは最適解問題で定番のナップサック問題。具体的な問題は Wikipedia のページ ナップサック問題 の図より拝借。
maximize
4 d4w12 + 2 d2w2 + 2 d2w1 + 1 d1w1 + 10 d10w4
subject to
weight: 12 d4w12 + 2 d2w2 + 1 d2w1 + 1 d1w1 + 4 d10w4 <= 15
binary
d4w12 d2w2 d2w1 d1w1 d10w4
ナップサック問題に入れる・入れないをフラグ的に binary で表現し、目的関数に価値の合計を、制約条件に重さの合計に対する上限値を設定すれば OK だろう。
ソルバーで解いて嬉しいかどうか別として、 N クイーン問題 をソルバーで解かせてみよう。
問題の記述が大変だったので、盤面の大きさは 4x4 にしてみる。
maximize
a1 + a2 + a3 + a4 + b1 + b2 + b3 + b4 + c1 + c2 + c3 + c4 + d1 + d2 + d3 + d4
subject to
rowA: a1 + a2 + a3 + a4 = 1
rowB: b1 + b2 + b3 + b4 = 1
rowC: c1 + c2 + c3 + c4 = 1
rowD: d1 + d2 + d3 + d4 = 1
col1: a1 + b1 + c1 + d1 = 1
col2: a2 + b2 + c2 + d2 = 1
col3: a3 + b3 + c3 + d3 = 1
col4: a4 + b4 + c4 + d4 = 1
diagC1D2: c1 + d2 <= 1
diagB1D3: b1 + c2 + d3 <= 1
diagA1D4: a1 + b2 + c3 + d4 <= 1
diagA2C4: a2 + b3 + c4 <= 1
diagA3B4: a3 + b4 <= 1
diagC4D3: c4 + d3 <= 1
diagB4D2: b4 + c3 + d2 <= 1
diagA4D1: a4 + b3 + c2 + d1 <= 1
diagA3C1: a3 + b2 + c1 <= 1
diagA2B1: a2 + c1 <= 1
binary
a1 a2 a3 a4 b1 b2 b3 b4 c1 c2 c3 c4 d1 d2 d3 d4
すごく冗長な表現になってしまっているように思えるが、ちゃんと解くことができる。いいね！
過去二回にわたって、
と綴ってきた、CTR, CVR の区間推定をするお話の総集編的なエントリです。
今回は commmons-math3 を使って、各区間推定方法における実際のカバレッジを測定するシミュレータを作ってみました。あわせて測定結果より、それぞれの方法の特性を確認していきます。
「二項比率の区間推定 (英語 Wikipedia)」は、統計的に独立・有限回で、各々の試行において「成功」か「失敗」のどちらかの結果が得られる問題の、「成功が発生する確率 $p$ (二項比率)」の信頼区間を求めることに相当します。
これをクリック率 (CTR) の区間推定に置き換えると、「成功」はすなわちクリックされたこと、「失敗」はクリックされなかったことに相当します。そして「成功が発生する確率」というのは、クリック回数 / インプレッション回数、つまりは CTR となります。
この「二項比率の区間推定」をする方法はいくつかあって、先の英語 Wikipedia にも記載されている方法として
などがあります。ただそれぞれ、論文で指摘されているとおりの特徴（というよりも欠点）があって、サンプルサイズの大きさや $p$ の値次第で実際のカバレッジが想定外の値になってしまう、という状況が起こりえます。
そのため、自身がとりあつかっている CTR や CVR のインプレッション数、クリック数、コンバージョン数に応じて、これらのいずれかをうまく使い分ける必要がある… と考えたのが、今回の一連のエントリを書き始めたきっかけでした。
それではそろそろシミュレーションをしてみて、各区間推定方法の特性を実際に確認してみましょう。
シミュレーションの内容は次のとおりです。
シミュレーションのプログラムは次のリンク先になります。
komiya-atsushi/binomial-proportion-confidence-interval
測定結果は Google スプレッドシート Comparison of Binomial proportion confidence interval にまとめました。
なお、以下の精度面の評価では Wald confidence interval の結果を省いています (Wald は評価するまでもなく突出して精度が悪く、グラフに含めると見辛くなってしまうため)。
まずはサンプルサイズの大きさを変化させたときの結果に着目してみます。
次に、二項比率 $p$ を変化させたときの結果を見てみます。
最後に、今回のシミュレーションをするのに要した時間を、計算方法ごとにグラフにしてみます。
論文を読んだり実際にシミュレーションをしてみて、それぞれの計算方法の特性が明らかになったので表にまとめてみます。
| Method | 実装が容易 | 区間の幅 | $p$ の偏り | 区間の上限/下限 |
|---|---|---|---|---|
| Wald | ◎ | × (*1) | × (*3) | × (*4) |
| Clopper-Pearson / exact | × | △ (*2) | - | - |
| Winson / score | × | - | × (*3) | - |
| Agresti-Coull / Adjusted Wald | ◯ | - | - | × (*4) |
元々の目的であるクリック率やコンバージョン率の区間推定に立ち返って、今回の調査・シミュレーション結果を見てみると、
と言えることでしょう。
今回の TokyoWebmining はここ最近、特に注目を集めている Deep Learning と word2vec がトピックということで、参加者募集も数分で枠が埋まってしまうほどの大人気っぷりでした。
せっかくなので、（参加したくても参加できなかった方々も多いことかと思いますし）たまにはまとめエントリを、会場内での質問やディスカッションを中心に書いてみようかと思います。
(聞き間違い・勘違いなどがあるかと思いますので、気づかれたかたはツッコミ願います)
「両方とも RNN って略すのやめろ」
前回 は CTR (クリック率)、CVR (コンバージョン率) に対するいくつかの区間推定方法を、それぞれの特徴とともに列挙してみました。今回はそれらの区間推定方法による実際の信頼区間を、Java や Python, R を用いて求める方法をまとめてみます。
Java で二項比率の区間推定をするには、commons-math3 の org.apache.commons.math3.stat.interval パッケージ以下のクラスを使うのが手っ取り早いでしょう。
それぞれの区間推定方法に対応するクラスは以下のとおりです。
NormalApproximationInterval クラスClopperPearsonInterval クラスWilsonScoreInterval クラスAgrestiCoullInterval クラスこれらのクラスに定義されている #createInterval(int numberOfTrials, int numberOfSuccesses, double confidenceLevel) を呼び出すことで信頼区間を求めることができます。たとえば numberOfTrials にインプレッション数を、 numberOfSuccesses にクリック数を、 confidenceLevel に 0.95 を設定して呼び出せば、95% 信頼水準での CTR の信頼区間が得られます。
サンプルコードは以下のとおり。
注意点として、 CTR や CVR の割合が 0% or 100% になるケース (つまりは、クリックやコンバージョンがまったく発生していない or 毎回発生している状態) において NormalApproximationInterval クラス、もしくは ClopperPearsonInterval クラスの #createInterval(int, int, double) メソッドで信頼区間を求めようとすると、 MathIllegalArgumentException や NotStrictlyPositiveException などの例外が発生してしまいます。
Wald confidence interval ではそもそも区間の幅が 0 になるケースに相当するのでどうしようもないのですが、Clopper-Pearson confidence interval では区間の片側だけでも算出することはできるはずです。
なので、このようなケースでも無理矢理に算出することはできるわけで、例えば ClopperPearsonInterval クラスの実装を以下のように修正すれば例外を生じることなく信頼区間を求めることができます。
Python の場合は、statsmodels を使います。
使い方は こちらのドキュメント を参考に、statsmodels.stats.proportion.proportion_confint(count, nobs, alpha, method) を呼び出します。CTR を算出するのであれば count にクリック数を、nobs にインプレッション数を指定し、加えて信頼水準 $100(1-\alpha)\%$ の $\alpha$ を alpha に指定します。
method には、区間推定方法を文字列で指定します。
normalbetawilsonagresti_coullサンプルコードは以下のとおり。
Java 同様に注意すべきこととして、CTR や CVR の割合が 0% or 100% の場合に、Clopper-Pearson confidence interval による区間推定の結果のうち、一方の片側が NaN になってしまうことが挙げられます。この場合、下側のエンドポイントが NaN であれば 0% と、上側のエンドポイントが NaN であれば 100% と読み替えればよいかと思います。
R では、binom パッケージ の binom.confintt(x, n, conf.level, methods, ...) を使って二項比率の区間推定をします。
methods に "all" を指定すると、すべての区間推定方法の結果を一覧で出力してくれます。もしくは以下の文字列を指定することで、対応する区間推定方法での結果を出力してくれます。
prop.testexactwilsonagresti-coull今回挙げた言語の各ライブラリでは、いずれもメソッド・関数を呼び出す程度の簡単なコードで二項比率の区間推定をすることができました。
ただ Java と Python については、CTR / CVR が 0% or 100% といったコーナーケースにおいてあまり好ましくない振る舞いをするため、多少の注意が必要となります。
また、各ライブラリで区間推定方法の呼称が異なることがあるため (Python / statsmodels の beta (Clopper-Pearson) や、R / binom の prop.test (Wald) など)、この点にも注意すべきかと思われます。
(次回こそはシミュレーション結果を…)
わけあってクリック率・コンバージョン率の信頼区間を算出したくなったのだけど、そのやり方を調べてみたら結構ややこしかったので、調べた結果をメモに残しておきます。
クリック率 (Click-through rate, CTR) やらコンバージョン率 (Conversion rate, CVR) を扱う仕事をしていると、少なくとも一度ぐらいはそれらの信頼区間を求めて (区間推定して) みたくなるものかと思います。
それというのも、例えば「100 回のインプレッションのうち、1 回のクリックが得られた」という標本 (サンプル) があったとして、これから CTR を点推定すると 1% になるものの、これは「サンプルサイズを増やしたときにも同様に 1% になるのか？」と言ったらそんなことは言えないわけで、ならば「どれくらいの信頼水準のときにどれくらいの範囲に真の CTR が存在しうるのか？」ということを知りたくなるわけです。この範囲を求めることがすなわち信頼区間を求める・区間推定をすることに相当します。
さて以上のようにクリック率・コンバージョン率の区間推定をしてみたいのですが、具体的にはどのようにすればいいのか？ これは二項比率 (binomial proportion, この日本語訳で適切なのか、わからない…) の区間推定をすることに等しくなります。
二項比率の区間推定をする方法について、あいにく日本語で網羅的にまとまった解説が Web 上には存在しないのですが (だからこのブログエントリを書いているわけでして…) 英語 Wikipedia のページ がそこそこ充実しているので、わかる人はこちらのページを合わせて参照することをおすすめします。以降はこの Wikipedia ページと、同ページでリファレンスされているいくつかの論文 (後述) をもとに話を進めます。
さて二項比率の区間推定をする方法として、今回は以下の 4 つを取り上げてみます (他にもベイズ的な確信区間などがあるわけですが、こちらは僕自身まだちゃんと理解しきれていないので割愛します)。
以下、それぞれの区間推定方法について、数式とともにその特徴などを列挙していきます。
$$\hat{p} \pm z_{\alpha/2} \sqrt{ \hat{p} (1 - \hat{p}) / n}$$
($\alpha$ は信頼係数、$z_{\alpha/2}$ は標準正規分布の上側 $100(\alpha/2)$ % 点、$n$ はサンプルサイズ or 試行回数、$\hat{p}$ は二項比率の推定値)
二項分布 $B(n,p)$ は正規分布 $N(np, np(1-p))$ で近似できることから、二項比率の信頼区間も正規分布 $N(p,p(1-p)/n)$ で近似することができるため、上記の式で信頼区間を求めることができます。
$$\left[1 + \frac{n-x+1}{x F_{2x,2(n-x+1),1-\alpha/2}} \right]^{-1} < p < \left[1 + \frac{n-x}{(x+1)F_{2(x+1),2(n-x),\alpha/2}} \right]^{-1}$$
($F_{n,m,z_{\alpha/2}}$ は、自由度 $n,m$ の F 分布における右側 $100\alpha$ 点)
Wald confidence interval は正規分布で近似することで信頼区間を求めていましたが、サンプルサイズが小さい場合や $p$ が 0 もしくは 1 に偏っている場合は正規分布での近似が難しくなります。その代わりに、F 分布を用いることで正確な (?) 信頼区間を求めることができるそうです。
$$\left( \hat{p} + \frac{z_{\alpha/2}^{2}}{2n} \pm z_{\alpha/2} \sqrt{[\hat{p}(1 - \hat{p}) + z_{\alpha/2}^{2} / 4n] / n} \right) / (1+z_{\alpha/2}^{2}/n)$$
Wald confidence interval も Clopper-Pearson confidence interval も、それぞれ信頼区間の幅については狭かったり広かったりしてちょっと扱いづらいわけですが、Wilson confidence interval ではその点においてバランスがとれた幅の信頼区間が求まるようです。
$$\tilde{p} \pm z_{a/2} \sqrt{ \tilde{p} (1 - \tilde{p}) / \tilde{n}}$$
$$(\tilde{n} = n + z_{\alpha/2}^{2},\ \tilde{p} = \frac{1}{\tilde{n}} \left(x + \frac{z_{\alpha/2}^{p2}}{2} \right) )$$
Wald confidence interval はサンプルサイズが小さい場合に信頼区間の幅が狭く、結果としてカバレッジが低下する問題がありました。一方でこの Adjusted Wald confidence interval では、信頼係数 $\alpha$ から定まる $z_{\alpha/2}$ を用いて $n,x$ を調整し、 $\tilde{n}, \tilde{p}$ を算出しています。そして、この $\tilde{n}, \tilde{p}$ を用いて、Wald confidence interval の式を使い、信頼区間の近似値を算出しています。
特に $\alpha$ が 0.05 の場合は adding two "successes" and two "failures" と言っているとおり、$n$ に $2 + 2 = 4$ を加え、$x$ に $2$ を加える操作をすればだいたいいい感じになってくれます。
ここまで区間推定方法について見てきましたが、これらを適用する先のクリック率やコンバージョン率の特性についても見ておきましょう。
長々と書いてきましたが、上記をふまえると、クリック率やコンバージョン率の区間推定には Agresti & Coull の Adjusted Wald confidence interval を利用するのがよいのではないか、と考えられます。
ただ厳密には、クリック率・コンバージョン率の特性を想定したテストケースをいくつか用意して、実際のカバレッジを測定するシミュレーションをしてみないことには胸を張って「○○ がいい！」とは言えないかな… と思います (シミュレーションは次回のブログエントリに書く予定)。
本ブログエントリは、主に以下 2 つの論文で述べられている内容をまとめたものとなっています。より詳しく知りたい方はこれらの論文を読まれることをおすすめします。
