2015-12-12

ロジスティック回帰を使って確率を予測したいときに「評価メトリクスとして使いたいのは AUC (areaUnderROC) じゃなくて Logarithmic loss (LogLoss) なんだよ!」と常々思っているのですが、現在の MLlib には二値分類 (BinaryClassificationEvaluator) 、多クラス分類 (MulticlassClassificationEvaluator)、回帰 (RegressionEvaluator) 用の Evaluator 実装しか用意されていなかったので、ついカッとなって実装してしまった次第です。

これくらいの機能は標準で取り揃えていてもいいんじゃないかと思うんだけどなあ…

2015-11-08

TL;DR

XGBoost で構築した予測モデルを Java から利用したい、それも特徴ベクトルが一つ一つ、任意のタイミングで与えられるような オンライン環境下リアルタイムな予測 を実現するために利用したい、という目的を叶えるためのモジュールを作りました。

(XGBoost の凄さとか XGBoost そのものの使い方とか GBDT/GBRT の解説は本エントリにはありませんので、そのような情報を求めている方は他のブログエントリを読まれることをおすすめします。)

xgboost4j という選択肢

Java から XGBoost を利用しようとすると、XGBoostをJavaのwrapperを使用して実行する - TASK NOTES にあるように、DMLC が提供している xgboost4j を利用する手段がすでに存在しています。ただ、この xgboost4j をオンライン予測に適用する場合、下記に挙げるようないくつかの懸念があります。

  • XGBoost の Java wrapper でしかないので、オンライン予測の目的で利用をするにはインタフェースがちょっと使いづらい
    • 大量の特徴ベクトルを入力して一括予測するようなバッチ処理に適したインタフェースになっている
  • LIBSVM フォーマットじゃないデータを入力するのに手間がかかる
    • 特に疎な特徴ベクトルを DMatrix で表現するのが面倒
  • JNI 由来のオーバーヘッドが気になる
    • 予測処理などが C++ で書かれているので高速処理が期待できる一方で、特徴ベクトルを一つ一つ与えて予測させる場合、ネイティブコードの呼び出しにかかるオーバーヘッドが全体のパフォーマンスに大きな影響を与えそう
  • どこの Maven repository にもアップロードされていないので、自前で mvn install する必要がある
    • OS X で開発をしている場合、ネイティブライブラリをビルドするのも一苦労 (参考)

Pure Java での予測を実現する

そういうわけで、

  • 速度性能がオンラインでの利用に耐えうる水準で
  • そこそこ使いやすいインタフェースで
  • ネイティブコードを必要としない (= Pure Java な)

XGBoost 互換な予測器を作って jCenter で公開しました。

使い方

README.md にもサンプルコード込みで使い方を書いていますので、合わせてご参照ください。

なお、タイトルでも「予測器」と明言しているとおり、学習 (training) 機能については割りきって一切実装をしていません (モデルの構築をオンラインですることは流石にないと思うので)。 そのため学習データを用意し、XGBoost を CLI で直接実行、もしくは Python や R のラッパーを経由するなどして別途モデルの構築を事前に済ませておく必要があります。

学習済みのモデルが用意できたら、そのモデルを new Predictor("/path/to/model-file") としてロードします。 予測をするには Predictor#predict() のメソッドを呼び出します。

予測の際の入力となる特徴ベクトルは、FVec インタフェースのオブジェクトとして表現する必要があります。 基本はお手持ちのデータの形式に合わせて FVec インタフェースを実装したクラスを用意していただくことになりますが、配列や Map でデータが表現されている場合には、以下のユーティリティメソッドを利用することができます。

  • double の配列で表現された密な特徴ベクトルで → 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

結果は上記のとおり、

  • 予測処理は 1ms 以下で処理できている
  • xgboost4j と比較して約 6,000 倍近い速度性能がでている (要は xgboost4j よりも十分速い)

となります。Pure Java にしただけで xgboost4j との性能差がこんなに出るものなのか… とちょっと不思議ではありますが、今回のベンチマークの計測に利用したテストデータは人工的に生成されたものなので、その影響があるのかもしれません。 そのため、実世界のデータを与えた場合にはまた違った結果になる可能性がありますのでご注意ください。

制限

現時点の xgboost-predictor は、モデルとしては “gbtree” のみをサポートしています。また目的関数は “binary” および “multi” のみのサポートとなります。 (“gblinear” や他の目的関数のサポートについては、必要そうであれば対応する予定です。)

2015-10-29

追記: 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 です。

2015-10-17

Approximate membership query (AMQ) が実現できるデータ構造としてそれなりに広く使われていそうな Bloom filter ですが、機能性を高めたバリエーションはそこそこ存在する一方で、空間効率を追求した・コンパクトなバリエーションはあんまり見つからないものです。

ここ最近、ふとデータ構造熱が高まってきたこともあったので、オリジナルの Bloom filter よりもコンパクトに表現できる代替データ構造を探してメモしてみました。

Compressed bloom filters

空間効率を追求するといっても、転送・永続化している状態での空間効率なのか、それともルックアップ可能な状態においての空間効率なのかで全然違うわけですが、この Compressed Bloom Filters は前者の転送時・永続化時の空間効率を選択したバリエーションです。

端的に言えば、Bloom filter のビット列における 0 と 1 の出現確率はランダム (1/2 = 0.5) というわけではなく偏りが生じることがあるため、算術符号で符号化すればよりコンパクトになるじゃん、というものです。 ただ、この圧縮はあくまでも転送・永続化時のことしか考えていないようで、Bloom filter が素早くルックアップできる状態にあるためには結局復号した状態で保持する以外に他ないようです。

Golomb-Compressed Sequence (or Golomb-code sets)

Cache-, Hash- and Space-Efficient Bloom Filters という論文で、Space-Efficient な手法として提案されている approximate membership query を提供するデータ構造です (Bloom filter と同じ機能を提供しているが、厳密には Bloom fitler ではない)。

要素のハッシュ値を算出するところまでは Bloom filter と同じですが、そのハッシュ値を

  1. 昇順ソートして
  2. そのソート順における、左に隣接する値との gap を求めて
  3. その gap を Golomb 符号で符号化

することで、Bloom filter よりもコンパクトな表現を実現しています (ハッシュ値は一様分布しているものとして、そのハッシュ値の gap をとると幾何分布になり、幾何分布に従う数値を符号化するには Golomb 符号が都合いい、ということ)。

この方式でのルックアップ操作は、あるハッシュ値が Golomb-Compressed Sequence に含まれるか否かだけを探索により判定すればいいので、Bloom filter のようなビット列への複数回のランダムアクセスは生じません。 とは言えど、この圧縮表現のままでは効率的なランダムアクセスが実現できないため、

  1. ハッシュ空間を大きさ $I$ の空間に分割して
  2. その分割された空間ごとに Golomb 符号化をし
  3. ルックアップの際はハッシュ値から探索先のハッシュ空間を絞り込んでから逐次 Golomnb 符号を復号する

ことで、復号にかかる時間効率を高められるようにしています。 $I$ を大きくすれば大きくするほど、空間効率がよくなるが探索効率は悪くなり、 $I$ を小さくすると空間効率は悪くなるが探索効率が上がる、というトレードオフの関係になっています。

Cuckoo filter

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 あたりを実装してみようかな…

2015-10-01

最適化問題、とりわけ線形計画問題が気になり始めるお年頃になってきたので、 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 を立ち上げるとインタラクティブシェルが立ち上がる。

read: 問題を記述したファイルを読み込む

ひとまず適当な問題をソルバーで解かせてみようと思い、講義テキストらしきもの (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: 最適化問題を解く

読み込んだ問題を実際に解いてみよう。問題を解くには 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] がポイント。最適解が求まったぽい。

display: 解を表示する

optimize コマンドで解が求まったら、 display のコマンドでその解を表示してみよう。

SCIP> display solution

objective value:                                27000
x1                                                 30   (obj:400)
x2                                                 50   (obj:300)

ちゃんと最適解が求まっているのがわかる。

write: 解をファイルに書き出す

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 クイーン問題

ソルバーで解いて嬉しいかどうか別として、 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

すごく冗長な表現になってしまっているように思えるが、ちゃんと解くことができる。いいね!

参考文献

宮代 隆平 の web ページ(整数計画法メモ)

2015-07-25

過去二回にわたって、

と綴ってきた、CTR, CVR の区間推定をするお話の総集編的なエントリです。

今回は commmons-math3 を使って、各区間推定方法における実際のカバレッジを測定するシミュレータを作ってみました。あわせて測定結果より、それぞれの方法の特性を確認していきます。

まずは復習から

「二項比率の区間推定 (英語 Wikipedia)」は、統計的に独立・有限回で、各々の試行において「成功」か「失敗」のどちらかの結果が得られる問題の、「成功が発生する確率 $p$ (二項比率)」の信頼区間を求めることに相当します。

これをクリック率 (CTR) の区間推定に置き換えると、「成功」はすなわちクリックされたこと、「失敗」はクリックされなかったことに相当します。そして「成功が発生する確率」というのは、クリック回数 / インプレッション回数、つまりは CTR となります。

この「二項比率の区間推定」をする方法はいくつかあって、先の英語 Wikipedia にも記載されている方法として

  • Wald confidence interval
  • Clopper-Pearson (or 'Exact') confidence interval
  • Wilson (or 'Score') confidence interval
  • Adjusted Wald confidence interval / Agresti-Coull confidence interval

などがあります。ただそれぞれ、論文で指摘されているとおりの特徴(というよりも欠点)があって、サンプルサイズの大きさや $p$ の値次第で実際のカバレッジが想定外の値になってしまう、という状況が起こりえます。

そのため、自身がとりあつかっている CTR や CVR のインプレッション数、クリック数、コンバージョン数に応じて、これらのいずれかをうまく使い分ける必要がある… と考えたのが、今回の一連のエントリを書き始めたきっかけでした。

シミュレーションしてみる

それではそろそろシミュレーションをしてみて、各区間推定方法の特性を実際に確認してみましょう。

シミュレーション内容

シミュレーションの内容は次のとおりです。

  • 次のサンプルサイズおよび真の二項比率 $p$ のすべての組み合わせごとに、1,000,000 回 の試行をします
    • p : 10% から 1% までの 1 ポイント刻みの値、および 1% から 0.1% までの 0.1 ポイント刻みの値
    • サンプルサイズ : 100 から 1,000 までの 100 刻みの値、および 1,000 から 10,000 までの 1,000 刻みの値
  • 信頼度を 95% としたときのカバレッジを測定します
    • このカバレッジが大きすぎず小さすぎず、 95% の値により近いことが望ましい結果となります
  • 1,000,000 回のサンプリングのうち、求められた信頼区間に $p$ の値が含まれている割合がカバレッジとなります
  • 精度面の検証以外にも、処理時間の測定をして評価をします

プログラム

シミュレーションのプログラムは次のリンク先になります。

komiya-atsushi/binomial-proportion-confidence-interval

結果と評価

測定結果は Google スプレッドシート Comparison of Binomial proportion confidence interval にまとめました。

なお、以下の精度面の評価では Wald confidence interval の結果を省いています (Wald は評価するまでもなく突出して精度が悪く、グラフに含めると見辛くなってしまうため)。

サンプルサイズを変化させたとき

まずはサンプルサイズの大きさを変化させたときの結果に着目してみます。

  • サンプルサイズを大きくすれば、いずれも 95% の信頼度に近いカバレッジが得られる
    • カバレッジのぶれもサンプルサイズが増えるごとに落ち着く
  • 全体的なカバレッジの傾向は、Clopper-Pearson > Agresti-Coull > Wilson となる
    • Clopper-Pearson のカバレッジが一番広く、95% を下回ることはまれ
    • Wilson は常に一番狭いカバレッジとなる
    • Agresti-Coull は、サンプルサイズが小さいときは Clopper-Pearson のように広いカバレッジをとり、サンプルサイズが増えるに従って Wilson のカバレッジに近づく

p を変化させたとき

次に、二項比率 $p$ を変化させたときの結果を見てみます。

  • サンプルサイズがそこそこある状態であれば、 $p$ が大きくなるにつれていずれも 95% の信頼度に近づいていく
    • こちらも、カバレッジのぶれは $p$ が大きくなることで落ち着く傾向となる
  • 全体的なカバレッジの傾向は、Clopper-Pearson > Agresti-Coull > Wilson となる
    • これはサンプルサイズを変化させたときと同じ傾向である
  • サンプルサイズあまりない状態で $p$ が小さい場合、Wilson は極端に狭いカバレッジ (90% 前後) となることがある

計算時間

最後に、今回のシミュレーションをするのに要した時間を、計算方法ごとにグラフにしてみます。

  • Wilson と Agresti-Coull はどちらも大差なし
  • Clopper-Pearson は他に比べて遅くなりがち
    • 他より複雑な計算が必要なため
    • サンプルサイズに比例して実行時間がかかってしまうのは、利用したライブラリ (commons-math3) の実装に起因している?
  • Normal Approximation (Wald) は、区間の幅が 0 になってしまうようなケースの計算時間が含まれていないため、Wilson, Agresti-Coull よりも計算時間が不当に短くなっている

まとめ

計算方法それぞれの特徴

論文を読んだり実際にシミュレーションをしてみて、それぞれの計算方法の特性が明らかになったので表にまとめてみます。

Method 実装が容易 区間の幅 $p$ の偏り 区間の上限/下限
Wald × (*1) × (*3) × (*4)
Clopper-Pearson / exact × △ (*2) - -
Winson / score × - × (*3) -
Agresti-Coull / Adjusted Wald - - × (*4)
  • *1 : 狭くなりがち
  • *2 : 広くなりがち
  • *3 : 狭くなる
  • *4 : 下限が 0 を下回る / 上限が 1 を超える

結論

元々の目的であるクリック率やコンバージョン率の区間推定に立ち返って、今回の調査・シミュレーション結果を見てみると、

  • Wald は欠点が多すぎてクリック率などの区間推定には適していない
  • クリック率のように、$p$ が極端に小さい値になりうることを考えると、Winson を利用するのは注意が必要
    • サンプルサイズが十分でない場合は使わない方がいい
  • 適切な幅の区間推定をする上では、Clopper-Pearson もしくは Agresti-Coull がよさそう
    • Clopper-Pearson は、意図して信頼区間を保守的に (慎重に) 広くとりたい場合に向いている
    • Agresti-Coull は下限が 0 より小さくなったり上限が 1 より大きくなったりする点に注意すれば (計算は楽だし) 使い勝手がいい

と言えることでしょう。

2015-01-24

今回の TokyoWebmining はここ最近、特に注目を集めている Deep Learning と word2vec がトピックということで、参加者募集も数分で枠が埋まってしまうほどの大人気っぷりでした。

せっかくなので、(参加したくても参加できなかった方々も多いことかと思いますし)たまにはまとめエントリを、会場内での質問やディスカッションを中心に書いてみようかと思います。

(聞き間違い・勘違いなどがあるかと思いますので、気づかれたかたはツッコミ願います)

深層学習時代の自然言語処理 by @unnonouno さん

Recurrent Neural Network

  • p.11 Recurrent Neural Network
    • 1 個の単語に対して、1 個の dense ベクトルがある
      • 次元数は 100〜1,000 ぐらい
      • RNN への入力ベクトルはこれになる
      • この単語ベクトルは各層で共有される
      • 正則化して sparse にしてみたけど、別によくはならなかった
  • p.14
    • 時間方向に層を重ねることに相当する
      • 層が文長に相当する、ここが他の NN と異なる
  • p.16
    • RNN の学習は、Back Propagation (Through Time) で学習する
      • 時間をさかのぼって学習しているように見える
      • 図的には、赤が Back propagation になる
  • p.17
    • 誤差が最初の方に伝搬しないという問題がある
    • これを解決したのが Long Short-Term Memory という技術
  • p.18 Long Short-Term Memory
    • 隠れ層のベクトル間でされる演算
    • 影響の与えるタイミングと与えないタイミングがあるはず、という考え
    • これによりパラメータが増えることになるが、これは学習対象となっている
    • 昨年後半からこれが流行ってる
    • 機械翻訳や Wikipedia 文章ぽいものの再生成、構文解析などに適用されてる by Google
    • 1 個の文がベクトルになって、そこから文が湧き出てくる

Recursive Neural Network

「両方とも RNN って略すのやめろ」

  • p.22
    • こちらは木構造を学習してつくる
    • 文章の構文解析だけでなく、画像の構造推定に使える
      • 木構造を作る
      • 二つの要素をくっつける順番を学習する
  • p.23 評判分析にも使われてる
    • こちらは木構造が事前に与えらて、ポジネガを判定している
    • 前は feature engineering で頑張って精度を高めていたけど、この RNN を使ったら良くなった
    • 文章全体ではなく、一つの文に対して適用される

本題:構文解析

  • p.32 Shift-Reduce 法
    • ガーデンパス文に弱い
    • これが Recurrent Neural Network に近いんじゃないか?
  • p.41
    • 品詞情報だけだと、構文解析は全然できない
  • p.46
    • 構文解析はまだルール作りはなんとかなるが、意味解析は膨大過ぎてルール作りは辛い

まとめ

  • p.48
    • Recurrent が流行っている
    • 音声認識、特に G とか MS とかの大手では DNN が使われているんじゃなかろうか?

ディスカッション

  • 単語分割、構文解析は精度が出ているが、意味解析とか談話解析はまだまだなので、そのあたりで DNN 使って精度出せるといいね
    • そもそもの問題設定があいまいだったりするけど…
  • テトリスブロックを回転させたもの同士が同じかどうかを判定するタスクで、DNN はうまくいかないというツッコミを入れた論文があった
  • 「言語学者をクビにすればするほど精度が上がる」

ディープラーニング徹底活用 画像認識編 by @atelierhide さん

この発表で一番言いたいこと

  • 学習済みモデルを徹底活用しよう!
    • Convolutional Neural Networks (CNNs) のモデル
  • 世界一のモデルを使うことができる
  • 1000 次元のベクトルが出力として得られる

Deep learning frameworks

  • 選択観点
    • 学習済みモデルが提供されているフレームワークを選ぶのがいい
    • Caffe よさそう
    • Caffe と DeCAF はほとんど違いはない
    • OverFeat は使い勝手がよくない
  • Caffe
    • 画像のリサイズなどは、まあまあフレームワークがよろしくやってくれる
      • 横長画像は正方形に変換されてしまうので、その点は注意しないといけない
    • Detection と Recognition は別
      • DNN が効くのは Recognition のほう
    • モデルの学習をする場合は、背景などが写り込んでいないものを選ぶべき
      • 分類はその限りではない

学習済みモデルの活用のアイデア

  1. 特徴抽出器として使う
    • CNNs の最後から 2 番目に得られる部分のベクトルを使う
      • pre-training 相当になっている
    • これを特徴量として、SVM などで分類する
  2. ファインチューニングをする
    • 出力層・分類数を入力画像にあわせて変更し、学習済みモデルのパラメータを最適化する
  3. 物体検出に使う

ディスカッション

  • みんなが Caffe を使い出していて、いろんな適用例が発表されはじめている
  • Caffe に Recurrent / LSTM が入るらしいということで、その手の界隈がざわついている

word2vec のご紹介 by @piroyoung さん

word2vec

Python で word2vec を使う

  • p.36
    • gensim を利用する
      • Paragraph vector も実装されている
    • コーパスをライブラリに喰わせるときに工夫が必要になる
      • ナイーブにやるとメモリが足りなくなる
      • 1 行読んでスペース区切り文を分割する… の処理をイテレーションさせる

素性のクラスタリング

  • word2vec で得られた単語のベクトルをクラスタリング
    • わりとよくクラスタリングできてる

QPR の学習

  • QPR = Quick purchase report, 消費者購買動向データ
  • p.56
    • ウィンドウサイズはものすごく大きなサイズにした
      • バスケット内の商品は、順番には意味がない
      • ただし先ほどの海野さんのツッコミにあるとおり、順番が考慮された結果となってしまった

ディスカッション

  • 文章中の助詞などを省いてみたら結果はどうなるの?
    • 単語同士の関係性を構成する要素になるので、動詞を入れても動詞が出てこなくなる
    • ものにもよるが、助詞を入れた方がいいであろう
  • 次元数 200 以外でやってみた?
    • 次元を上げて、悪くなることはなかった
    • 計算時間はその分かかる

今回の TokyoWebmining の所感

  • Deep learning、いまいちちゃんと理解できてなかったけど、雰囲気はだいぶつかめてきた
  • 画像の取り扱いにおいては Deep learning を利用するのがもはや当たり前っぽい
    • ImageNet の学習済みモデル、応用の幅が広いね!
  • 自然言語処理での Deep learning 活用、研究の進展が速いので、常にキャッチアップしていかないと置いて行かれそう…
  • とにかく Deep learning 熱の高まりっぷりがはんぱない!

2015-01-19

前回 は CTR (クリック率)、CVR (コンバージョン率) に対するいくつかの区間推定方法を、それぞれの特徴とともに列挙してみました。今回はそれらの区間推定方法による実際の信頼区間を、Java や Python, R を用いて求める方法をまとめてみます。

Java による区間推定

Java で二項比率の区間推定をするには、commons-math3org.apache.commons.math3.stat.interval パッケージ以下のクラスを使うのが手っ取り早いでしょう。

それぞれの区間推定方法に対応するクラスは以下のとおりです。

  • Wald confidence interval
    • NormalApproximationInterval クラス
  • Clopper-Pearson (or 'Exact') confidence interval
    • ClopperPearsonInterval クラス
  • Wilson (or 'Score') confidence interval
    • WilsonScoreInterval クラス
  • Adjusted Wald confidence interval / Agresti-Coull confidence interval
    • 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) メソッドで信頼区間を求めようとすると、 MathIllegalArgumentExceptionNotStrictlyPositiveException などの例外が発生してしまいます。

Wald confidence interval ではそもそも区間の幅が 0 になるケースに相当するのでどうしようもないのですが、Clopper-Pearson confidence interval では区間の片側だけでも算出することはできるはずです。

なので、このようなケースでも無理矢理に算出することはできるわけで、例えば ClopperPearsonInterval クラスの実装を以下のように修正すれば例外を生じることなく信頼区間を求めることができます。

Python による区間推定

Python の場合は、statsmodels を使います。

使い方は こちらのドキュメント を参考に、statsmodels.stats.proportion.proportion_confint(count, nobs, alpha, method) を呼び出します。CTR を算出するのであれば count にクリック数を、nobs にインプレッション数を指定し、加えて信頼水準 $100(1-\alpha)\%$ の $\alpha$ を alpha に指定します。

method には、区間推定方法を文字列で指定します。

  • Wald confidence interval
    • normal
  • Clopper-Pearson (or 'Exact') confidence interval
    • beta
  • Wilson (or 'Score') confidence interval
    • wilson
  • Adjusted Wald confidence interval / Agresti-Coull confidence interval
    • agresti_coull

サンプルコードは以下のとおり。

Java 同様に注意すべきこととして、CTR や CVR の割合が 0% or 100% の場合に、Clopper-Pearson confidence interval による区間推定の結果のうち、一方の片側が NaN になってしまうことが挙げられます。この場合、下側のエンドポイントが NaN であれば 0% と、上側のエンドポイントが NaN であれば 100% と読み替えればよいかと思います。

R による区間推定

R では、binom パッケージbinom.confintt(x, n, conf.level, methods, ...) を使って二項比率の区間推定をします。

methods"all" を指定すると、すべての区間推定方法の結果を一覧で出力してくれます。もしくは以下の文字列を指定することで、対応する区間推定方法での結果を出力してくれます。

  • Wald confidence interval
    • prop.test
  • Clopper-Pearson (or 'Exact') confidence interval
    • exact
  • Wilson (or 'Score') confidence interval
    • wilson
  • Adjusted Wald confidence interval / Agresti-Coull confidence interval
    • agresti-coull

まとめ

今回挙げた言語の各ライブラリでは、いずれもメソッド・関数を呼び出す程度の簡単なコードで二項比率の区間推定をすることができました。

ただ Java と Python については、CTR / CVR が 0% or 100% といったコーナーケースにおいてあまり好ましくない振る舞いをするため、多少の注意が必要となります。

また、各ライブラリで区間推定方法の呼称が異なることがあるため (Python / statsmodels の beta (Clopper-Pearson) や、R / binom の prop.test (Wald) など)、この点にも注意すべきかと思われます。

(次回こそはシミュレーション結果を…)

2015-01-07

わけあってクリック率・コンバージョン率の信頼区間を算出したくなったのだけど、そのやり方を調べてみたら結構ややこしかったので、調べた結果をメモに残しておきます。

はじめに

クリック率 (Click-through rate, CTR) やらコンバージョン率 (Conversion rate, CVR) を扱う仕事をしていると、少なくとも一度ぐらいはそれらの信頼区間を求めて (区間推定して) みたくなるものかと思います。

それというのも、例えば「100 回のインプレッションのうち、1 回のクリックが得られた」という標本 (サンプル) があったとして、これから CTR を点推定すると 1% になるものの、これは「サンプルサイズを増やしたときにも同様に 1% になるのか?」と言ったらそんなことは言えないわけで、ならば「どれくらいの信頼水準のときにどれくらいの範囲に真の CTR が存在しうるのか?」ということを知りたくなるわけです。この範囲を求めることがすなわち信頼区間を求める・区間推定をすることに相当します。

さて以上のようにクリック率・コンバージョン率の区間推定をしてみたいのですが、具体的にはどのようにすればいいのか? これは二項比率 (binomial proportion, この日本語訳で適切なのか、わからない…) の区間推定をすることに等しくなります。

二項比率の区間推定をする方法について、あいにく日本語で網羅的にまとまった解説が Web 上には存在しないのですが (だからこのブログエントリを書いているわけでして…) 英語 Wikipedia のページ がそこそこ充実しているので、わかる人はこちらのページを合わせて参照することをおすすめします。以降はこの Wikipedia ページと、同ページでリファレンスされているいくつかの論文 (後述) をもとに話を進めます。

二項比率の区間推定をするいくつかの方法

さて二項比率の区間推定をする方法として、今回は以下の 4 つを取り上げてみます (他にもベイズ的な確信区間などがあるわけですが、こちらは僕自身まだちゃんと理解しきれていないので割愛します)。

  • Wald confidence interval
  • Clopper-Pearson (or 'Exact') confidence interval
  • Wilson (or 'Score') confidence interval
  • Adjusted Wald confidence interval / Agresti-Coull confidence interval

以下、それぞれの区間推定方法について、数式とともにその特徴などを列挙していきます。

Wald confidence interval

$$\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)$ で近似することができるため、上記の式で信頼区間を求めることができます。

  • 数式がわりと容易である
    • つまりは実装するのも比較的楽、ということ
  • この手法により求まる信頼区間は、信頼水準から得られるそれよりも狭くなる (= 実際のカバレッジが低くなる) 傾向にある
    • サンプルサイズが小さい ($n < 100$ ぐらいの) 場合に、特にその傾向が表れる
    • また、 $p$ が 0.5 から 0 もしくは 1 に偏っているほどに顕著になる
  • 下側信頼限界が負数に、もしくは上側信頼限界が 1 を超える場合がある
    • それぞれ、 $p$ が 0 に近い場合、1 に近い場合にそのような状況になる
  • 成功回数 $x$ が、$x=0$ や $x=n$ の場合は、信頼区間を求めることができない (幅が 0 の区間になる)
    • クリック数が 0、もしくはインプレッション数に等しい場合が該当する

Clopper-Pearson (or 'Exact') confidence interval

$$\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 分布を用いることで正確な (?) 信頼区間を求めることができるそうです。

  • サンプルサイズが小さくても、求まる信頼区間は Wald confidence interval のように狭くはなく、比較してカバレッジがよい
    • むしろ逆に、ちょっと広すぎる…

Wilson (or 'Score') confidence interval

$$\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 ではその点においてバランスがとれた幅の信頼区間が求まるようです。

  • 数式が複雑である
    • 実装する際にエンバグしやすい ※参考
  • 得られる信頼区間は狭すぎず、広すぎず
  • サンプルサイズによらず、$p$ が 0 もしくは 1 に偏っている場合に信頼区間が狭くなる傾向がある

Adjusted Wald confidence interval / Agresti-Coull 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$ を加える操作をすればだいたいいい感じになってくれます。

  • 数式はまだ容易な方ではある
  • サンプルサイズが小さい場合であっても、信頼区間が狭くなるようなことはない
    • Wilson confidence interval 同様にバランスのとれた幅の信頼区間が求まる
  • 下側信頼限界が負数に、もしくは上側信頼限界が 1 を超えうる問題は健在している

クリック率・コンバージョン率の特性

ここまで区間推定方法について見てきましたが、これらを適用する先のクリック率やコンバージョン率の特性についても見ておきましょう。

  • クリック率
    • サンプルサイズ (インプレッション数) は十分な大きさとなる
      • コンテキストを考慮した CTR を算出する場合はその限りではない
    • 値は常に小さな値になりがちで、0 に近くなる
      • 1% に満たないことも十分にあり得る
  • コンバージョン率
    • サンプルサイズ (クリック数) が小さいことがある
    • 値はまちまちで、一桁 % のときもあれば二桁 % になることもある
  • 共通して言えること
    • 0% となるケースを考慮するべき

現時点での推測

長々と書いてきましたが、上記をふまえると、クリック率やコンバージョン率の区間推定には Agresti & Coull の Adjusted Wald confidence interval を利用するのがよいのではないか、と考えられます。

ただ厳密には、クリック率・コンバージョン率の特性を想定したテストケースをいくつか用意して、実際のカバレッジを測定するシミュレーションをしてみないことには胸を張って「○○ がいい!」とは言えないかな… と思います (シミュレーションは次回のブログエントリに書く予定)。

参考文献

本ブログエントリは、主に以下 2 つの論文で述べられている内容をまとめたものとなっています。より詳しく知りたい方はこれらの論文を読まれることをおすすめします。