中野智文のブログ

データ・マエショリストのメモ

「被害者が犯行可能な人を犯人の候補」とするWald法

背景

二項分布やベルヌーイ試行の確率の信頼区間を求める方法として、Wald法というものがある。 非常に気持ち悪いと思うのだが、まだ多くの人がその気持ち悪さを理解していないと思うので、自分の力不足を感じながらも書いた。

Wald 法

正規分布は n\mu > 5, n\sigma^2 > 5のとき、ベルヌーイ試行の確率によい近似となるらしい。 ここで、 \mu, \sigma^2 はそれぞれ、母平均、母分散、 n は試行回数。

正規分布 N(\mu, \sigma^{2}/n) の信頼区間

\displaystyle \mu \pm z\sqrt{\sigma^{2}/n}

である。ここで、 z は正規分布分位数。 ベルヌーイ試行の成功確率の分散は \mu(1-\mu) と平均を使って表すことができるので、 正規分布で近似したベルヌーイ試行の母平均の信頼区間は

\displaystyle \mu \pm z\sqrt{\mu(1-\mu)/n}

となる（一つ目の近似。分布の近似）。

さらに 試行回数が多くなれば 中心極限定理により、標本平均は母平均に近づく。

よって母平均を標本平均と標本分散に置き換える（二つ目の近似。今度は値の近似）。 すなわち、 m を標本平均とする。よって、母平均を標本平均と近似し、さらに分布も正規分布で近似した、 ベルヌーイ試行の成功確率の母平均の信頼区間は \displaystyle m \pm z\sqrt{m(1-m)/n}

となる。これが Wald 法。

Wilson score interval

母平均が分かっていたなら、標本平均はどのように分布するだろうか。 標本平均の分布は、正規分布 N(\mu, \sigma^{2}/n) に従うと仮定すると、確信区間*1

\displaystyle \mu \pm z\sqrt{\mu(1-\mu)/n}

に入ることになる（一度目の近似）。

ある母平均 \mu を考えたとき、観測された標本平均 m がこの確信区間の外だったとき、その母平均 \muは信頼区間の外であると考える。

これが、Wilson score interval の考えである。

続けると、この信頼区間の境界は次のようにして求める。

次の式を満たすとき、標本平均 mは確信区間に入っている。

\displaystyle |\mu - m| \le z \sqrt{\mu(1-\mu)/n}

両辺二乗して、

\displaystyle (\mu - m)^2 \le z^2 \mu(1-\mu)/n

が得られる。ところで、区間の境界は等号が成り立つところである。すなわち、

\displaystyle (\mu - m)^2 = z^2 \mu(1-\mu)/n

である。これは \mu についての二次方程式なので、それを解くと、

\displaystyle \mu = \frac{m+t/2 \pm \sqrt{ m(1-m)t + t^2/4} }{ 1+t }

ここで、 t=z^2/n (論文の式)

これが Wilson score interval である。

考察

Wilson score interval を考えると Wald 法の気持ち悪さが浮き出てくる。

母平均の信頼区間を求めるのに、母平均を標本平均で近似している。

犯人を推定するのに、被害者を一旦犯人とおいて、被害者が犯行可能な人を本当の犯人の候補とするのがWald法である。

それでうまくいことがあるとしても気持ち悪すぎる。

参考

http://www.barestatistics.nl/uploads/1/1/7/9/11797954/wilson_1927.pdf

*1:信頼区間と区別するために敢えてこのように書く

mac が重くてたまに熱か何かで落ちたのでその時のメモ

背景

mac が重くてたまに熱か何かで落ちる。vagrant で動いているものがあるので、それが原因かと思ったが、 top で確認すると、mds_stores というプロセスが動いていた。

対応

下記のページを参考にしようとした。

blog.tottokug.com

ところが、

sudo locate .git

で気がついたが、locate が動いていない。 普通は locate を使うので、エラーメッセージに従い、下記のコマンドを実行

sudo launchctl load -w /System/Library/LaunchDaemons/com.apple.locate.plist

そして、locate あれば、spotlight は使わないので、spotlight を停止することに。

qiita.com

を参考に、

 sudo mdutil -a -i off

これで止まった。（はず）

念のため、top コマンドで確認すると、mds_stores　は消えたものの、 mdworker というのが数個生きているみたい。

qiita.com

を参考に、

sudo launchctl unload -w /System/Library/LaunchDaemons/com.apple.metadata.mds.plist

を実行。確認すると消えていた。

正規分布の分位点関数（パーセント点関数）の近似

背景

正規分布の分位点（パーセント点）は数値計算ライブラリの分位点関数を使えば簡単に求められる。 例えば、python だと、scipy.stats.norm の、ppt だ。 ところが、SQL には標準でないようである。 自分で関数を作ってもいいが、ここはSQLでかける近似式をつかってもいいかもしれない。

近似式を調べる

いろいろ調べたが、Hastings 1955 は有名だが、その後多くの近似式があるようである。 全部を調べるは大変で、サーベイっぽい論文などもあるのだが、最終的には、多く場合分けをして、多く次数を増やしていけば、近似精度がよくなるという話。 これはドメイン依存で、すなわち、どのような状況でどの程度近似精度が重要なのかと、SQLなどで記述する場合に複雑になっていってもいいのかというトレードオフの問題がある。 あと、近似式だから、正解はなく、多くの大量に考え出されているので、それら全てを評価するのもコストがかかる。

というわけで、私が探した中で、一番シンプルな式の山内 1965から検証する。

正規分布の分位点関数の近似式 山内 1965

戸田 英雄, 小野令美 「正規分布のパーセント点を求めるミニマックス近似」(1993) によると式は

\begin{align} \left[ y \left( 2.0611786 - \frac{5.7262204}{y + 11.640595} \right) \right] ^{ \frac{1}{2} }, \\ y = - \log 4P(1-P) \end{align}

らしい。 y の部分がまとめられてしまっているが係数も3つで乗数もなく、入力の間違いの確認も容易である。 二乗根と対数ならSQLもサポートしている。

python での近似性能の検証

多くの場合分位点そのものに興味はない。 だから分位の近似の評価ではなく、近似された分位点から逆に計算された分位*1が元の分位と比べてどうだったかを調べる。

import numpy as np

def yamauchi1965(q):
    y = -np.log(4*q*(1-q))
    return np.sqrt(y*(2.0611786-5.7262204/(y+11.640595)))

def yamauchi1965_sql(q):
    return np.sqrt(-np.log(4*q*(1-q))*(2.0611786-5.7262204/(-np.log(4*q*(1-q))+11.640595)))
    # sqrt(-log(4*Q*(1-Q))*(2.0611786-5.7262204/(log(4*Q*(1-Q))+11.640595)))

# 検証用
from scipy.stats import norm
for i in [1, 10, 20, 50, 100, 1000, 10000, 100000]:
    q = 1.0-1.0/2/i
    z = norm.ppf(q)
    z1 = yamauchi1965(q) # これが近似の分位点z1
    q1 = norm.cdf(z1)    # 近似の分位点z1から、元の分位を逆計算q1
    r = 1-(1-q1)/(1-q)   # 分位(1-q)と近似から逆計算した分位(1-q1)の割合の違い（割合）からさらに１を引いて相対的な差を求める
    print("q: %f, z: %f, z1: %f, q1: %f, r: %f" % (q, z, z1, q1, r))

実行結果

q: 0.500000, z: 0.000000, z1: -0.000000, q1: 0.500000, r: 0.000000
q: 0.950000, z: 1.644854, z1: 1.645636, q1: 0.950081, r: 0.001613
q: 0.975000, z: 1.959964, z1: 1.960594, q1: 0.975037, r: 0.001471
q: 0.990000, z: 2.326348, z1: 2.326358, q1: 0.990000, r: 0.000027
q: 0.995000, z: 2.575829, z1: 2.575264, q1: 0.994992, r: -0.001635
q: 0.999500, z: 3.290527, z1: 3.288966, q1: 0.999497, r: -0.005561
q: 0.999950, z: 3.890592, z1: 3.890511, q1: 0.999950, r: -0.000335
q: 0.999995, z: 4.417173, z1: 4.420750, q1: 0.999995, r: 0.016409

0.5 がエラーになると思っていたが、特に問題はなさそう。 元の分位に戻したときの相対的な差は多くの場合1％以下に抑えられている。 個人的はこれで満足である。

ただし、標準正規分布表と比べたときに、一番下の段の方では、違いが出てくる。

まとめ

山内 1965

sqrt(-log(4*Q*(1-Q))*(2.0611786-5.7262204/(log(4*Q*(1-Q))+11.640595)))

ただし、Q は分位（割合・確率）。

ちなみに、mysql では、これらの整数を4 の代わりに 4.000000 としてやらないと、計算してくれないことがあるので注意

*1:分位という言葉は正しいかわからないが、ここでは分位点を調べるための割合、確率をいう

scikit learn の Kfold, StratifiedKFold, ShuffleSplit の違い

背景

scikit learn の cross validation にて、テスト事例の分割方法に Kfold, StratifiedKFold, ShuffleSplit, StratifiedShuffleSplit というのがある（他にもある）。 その違いがわかりにくい。

Kfold と StratifiedKFold の違い

Kfold は知っているという前提で、StratifiedKFold は、その説明にある通り、

The folds are made by preserving the percentage of samples for each class.

その分割は、サンプルのそれぞれのクラスの割合を保ったまま作られる（超意訳）。

ということらしい。100個のデータがあって、クラス１とクラス０に、それぞれ5個と95個属するとき、それぞれのクラスの割合は、5%と95％である。 よって、その5%と10％を保ったまま分割されるということだ。

でもこのデータを10-Fold にする場合はどうなるのだろうか？ 10個づつになるわけだから、0.5個 と 9.5個 に分割しないとクラスの割を保ったままに出来ないはずだが、そのような分割はできない。 だから、そもそもこのようなマイナーなクラスが小さなデータの場合、理想的にクラスの割合が等しく分割することは出来ない。 マイナーなクラスのための機能のような気がするのに、だ。 結果として、

The least populated class in y has only 5 members, which is too few. The minimum number of labels for any class cannot be less than n_folds=10.

といったようなメッセージを見るかもしれない。 これはコード内のコメントを読むと、例外を出す代わりに trim すると言っている。 コードを読み切れていないが、すでに k-fold とはよべないような状況が起きている気がする。 （例えばマイナーなクラスの事例が別の分割(fold)にも重複して使われたり…）

KFold と ShuffleSplit の違い

こちらも、 KFold は知っているという前提で話をすすめる。 KFoldのとの大きな違いは ShuffleSplit は繰り返し再サンプリングする。 よって、k 回行うと重複して同じ事例がでてくるかもしれない。 この ShuffleSplit は事例の数を減らして（ダウンサンプリングして）学習させる際に、別のダウンサンプリングでも結果がそれほど変わらないかを検証する際に使うらしい。

だから、KFold の オプション shuffle=True と同じ結果になるわけではないことに注意。

ちなみに shuffle オプションはデフォルトで、False なので、 shuffle=True をつけないと、事例をランダムにソートしてくれない。

まとめ

  • クラスの割合を保ったまま分割したときは、StratifiedKFold を使う。
  • ただし、マイナーなクラスの事例数が少ない（特に分割数より少ない）と警告が出て、その場合は想定外の結果となるだろう。
  • ShuffleSplit は KFold の shuffle=True ではなく、特にデータが大きい場合のダウンサンプリングのばらつきなどを検証するために使う（複数の分割において重複するかもしれない）

なお、GridSearchCVの説明 によると、推定器か教師データが二値か多クラスの場合、デフォルトで、StratifiedKFold になるそうである。 もちろん、shuffle は Falseになっている。 scikit learn は初心者向けのようで実は…

参考

stackoverflow.com

https://www.kaggle.com/ogrellier/kfold-or-stratifiedkfold

scikit learn の GridSearchCV で検証事例に class_weight（sample_weight) をつける

背景

GridSearchCV で検証事例に sample_weight をつけるような引数はまだ存在しない。

github.com

でも使いたい。（使わないとうまく行かねー）

metrics 関数自体は sample_weight に対応しているんだよね〜。

対応

対応した metrics 関数をつくり、make_scorer で呼び出そう。

以前の記事

nakano-tomofumi.hatenablog.com

にちなんだ例で作る。

from sklearn.model_selection import GridSearchCV
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import KFold
from sklearn.metrics import accuracy_score, log_loss, make_scorer
from sklearn.utils import class_weight 


def balanced_log_loss(y_true, y_pred, cw, eps=1e-15, normalize=True, labels=None):
    sample_weight = [cw[y] for y in y_true]
    return log_loss(y_true=y_true,
                              y_pred=y_pred,
                              eps=eps,
                              sample_weight=sample_weight,
                              labels=labels)



...
    labels = [0, 1] # np.unique(y)
    cw_array = class_weight.compute_class_weight(y=y,
                                                                                    class_weight='balanced',
                                                                                    classes=labels)
    cw = dict(zip(labels, cw_array))
    clf = GridSearchCV(
        LogisticRegression(penalty='l1',
                           class_weight='balanced',
                           fit_intercept=False),
        {'C': [0.01, 0.1, 1, 10, 100]},
        cv=KFold(n_splits=20, shuffle=True),
        scoring={'neg_log_loss': make_scorer(balanced_log_loss, 
                                                                        labels=labels,
                                                                        cw=cw, 
                                                                        greater_is_better=False), 
                 'accuracy': 'accuracy'},
        n_jobs=-1,
        refit='neg_log_loss'
    )
    clf.fit(X, y)
...

ちょっとだけ説明すると、balanced_log_loss という独自の metrics 関数を作っている。 その引数に、dict 型の cw というものを用意した。class_weight の略だ。 そして、評価される真の値(y_true)に基づいて、重みを与えている。 もし、予想の値に対して重みを与えたければ、y_pred を[cw[y] for y in y_true]の中のy_trueと書き換えれば良い。

さて、この cw だが、学習オプションのclass_weight='balanced'を設定すると内部で呼ばれる compute_class_weight を使っている。 今回は訓練事例全体で cw を求めたが、cross validation 中の訓練事例のみで求めることは難しいようである。 検証事例のみで求めたければ、先の自作関数の中で y_true や y_pred から作れば良いのだが、それはナンセンスだろう。

このようなclass weight と、先に作った自作関数を、make_scorer で呼び出す。 今回も neg_log_loss に変換するために、greater_is_better=False オプションを付与して呼び出した。

まとめ

  • 学習時だけでなく、cross validation の評価時にも、sample_weight を使いたい。
  • GridSearchCV にはそのようなオプションはない。
  • そこで、make_scorer から呼び出すことができる自作の metrics 関数にそのようなオプションをつけて対応する。

正規分布の分位数（標準正規分布）を python で求める

背景

あの Z=1.96 などの数値は、正規分布の分位数とか、分位点とか、パーセントを使っている場合には、パーセント点とか、百分位点とか、パーセンタイル とかよばれる。 さてこの正規分布の分位数であるが、計算環境が十分でない状況では標準正規分布表というものを使って求めるのが普通だったかもしれない（少なくとも昔の統計の本にはそのようなヒョウがついてくるのが普通だった）が、今なら直接計算できる。

正規分布の分位数の計算

片側95%点（両側だと90%）の求め方

from scipy.stats import norm
norm.ppf(0.95)

片側97.5%点（両側だと95%）の求め方

from scipy.stats import norm
norm.ppf(0.975)

参考

stackoverflow.com

ValueError: y_true contains only one label (0). Please provide the true labels explicitly through the labels argument. というエラー

背景

GridSearchCV をすると、

ValueError: y_true contains only one label (0). Please provide the true labels explicitly through the labels argument.

というようなエラーに遭遇することがある。 これは「評価用のデータに１つの正解ラベルしか含まれていないから、複数の正解ラベルを labels 引数に入れてね」というメッセージなのだが、どこに入れるか分からない。

labels は score 関数のラベルに入れる

結論は表題の通りなのだが、scoring のパラメータに入れる関数のパラメータとなる。scoring パラメータは、string, list, dict, None の形式に対応しているが、 string や list は、score関数の文字列で、引数を指定できない。よってdict形式による関数の指定ということなる。 dict 関数の指定は、文字列と関数の二種類があるが、この例のprecision のようにstringで指定したら引数を指定できないので、 make_scorer(accuracy_score) のように指定する。この make_scorer は、 sklearn.metrics.make_scorer — scikit-learn 0.19.1 documentation のとおりだが、第一引数に、メトリックス関数を指定できるが、その他にそのメトリックス関数の引数（**kwargsと書かれている）を指定できる。

例

neg_log_loss に labels=[0, 1] に引数で与える場合。accuracy は一致をみるので、他のクラスの情報をしらなくても問題ないが、neg_log_loss は他に何クラスあるのか知っていないと計算できない模様。

from sklearn.model_selection import GridSearchCV
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import KFold
from sklearn.metrics import accuracy_score, log_loss, make_scorer

...

    clf = GridSearchCV(
        LogisticRegression(penalty='l1')
        {'C': [0.01, 0.1, 1, 10, 100]},
        cv=KFold(n_splits=20),
        scoring={'neg_log_loss': make_scorer(log_loss, labels=[0, 1], greater_is_better=False), 
                 'accuracy': 'accuracy'},
        n_jobs=-1,
        refit='neg_log_loss'
    )
...

greater_is_better=False は、log_loss から、neg_log_loss に変換するときに必要。

まとめ

make_scorer 関数を使って、labels 引数を与えたものを scoring 引数のdict型で与える。

ちなみに、cross validation に、KFold を使わなかったり、cv数(n_splits)を小さくすれば、このエラーを回避することもできる場合もある（検証事例が単独のクラスになりくくなる）が、それは誤った解決法である。

参考

github.com