Negative/Positive Thinking

2014-11-19

ベータ分布のquantile

はじめに

boostには、確率分布のquantile(分位数、分布をp:1-pに分割する点)を計算するものが用意されている。
ベータ分布の場合について、自分でも書いてみる。

コード

#include <iostream>
#include <cmath>
//#include <boost/math/distributions/beta.hpp>

class BetaDistribution {
  double eps;
  double val_a, val_b;

  //Ix : regularized incomplete beta function
  double Ix(double x, double a, double b){
    if(a <= 0) return -1;
    if(b <= 0){
      if(x < 1) return 0;
      if(fabs(x-1) < eps) return 1.0;
      return -1;
    }

    if(x > (a+1) / (a+b+2)) return 1 - Ix(1-x, b, a);
    if(x <= 0) return 0;
    double p1 = 0, q1 = 1;
    double p2 = exp(a*log(x) + b*log(1-x) + lgamma(a+b) - lgamma(a) - lgamma(b)) / a, q2 = 1;
    double prev, d;
    for(int k=0; k<500;){
      prev = p2;
      d = -((a+k) * (a+b+k) * x) / ((a+2*k) * (a+2*k+1));
      p1 = p1*d + p2;
      q1 = q1*d + q2;
      k++;
      d = (k * (b-k) * x) / ((a+2*k-1) * (a+2*k));
      p2 = p2*d + p1;
      q2 = q2*d + q1;
      if(fabs(q2) < eps){
        q2 = 1e+9;
        continue;
      }
      p1 /= q2;
      q1 /= q2;
      p2 /= q2;
      q2 = 1;
      if(fabs(p2-prev) < eps) return p2;
    }
    return -1;
  }
public:
  BetaDistribution(double val_a, double val_b):
    val_a(val_a),
    val_b(val_b),
    eps(1e-5) 
  {}
  
  //Int_0^Q { dPhi(x) } >= p
  double quantile(double p){
    double lb = 0.0, ub = 1.0;
    for(int i=0; i<200; i++){
      double m = (lb+ub) / 2.0;
      if(Ix(m, val_a, val_b) < p) lb = m;
      else ub = m;
    }
    return lb;
  }
};

int main(){
  double a, b, p;

  while(std::cin >> a >> b >> p){
    //boost
    //boost::math::beta_distribution<> dist(a, b);
    //std::cout << quantile(dist, p) << std::endl;
    
    //mybeta
    BetaDistribution mybeta(a, b);
    std::cout << mybeta.quantile(p) << std::endl;
  }

  return 0;
}

regularized不完全ベータ関数Ix(a,b)は「C言語による最新アルゴリズム辞典」で紹介されている連分数展開による方法。
quantileの計算は、2分探索で探している・・・
(一応、boostの結果と値はだいたい同じみたいだが、pが0や1に近い点だと誤差が大きくなっている模様)

参考

2014-10-15

Friedman testとNemenyi testメモ

はじめに

複数のアルゴリズムの結果の有意差検定に使用されていたので、メモ。

より詳細に紹介されているのは以下の論文
Demsar, Statistical Comparisons of Classifiers over Multiple Data Sets, 2006
http://jmlr.csail.mit.edu/papers/volume7/demsar06a/demsar06a.pdf
1999年から2003年のICMLの論文で使われた検定について調査して、アルゴリズムやデータセットの種類によって使うべき検定方法を紹介している。

Friedman test

  • 対応のある順序尺度の3群以上の平均順位に有意差があるかどうかを検定(ノンパラメトリック検定)
    • repeated-measures ANOVAのノンパラメトリックに相当
  • 帰無仮説(H0) : μ_0=μ_1=…=μ_k (すべての平均順位が等しい)
  • 検定方法
    • χ^2_F = 12n/{k(k+1)} Σ_{i=1}^k { (R_i)^2 - k(k+1)^2 / 4 }を計算
    • 自由度(k-1)のχ^2検定をする。nやkが小さい場合はχ^2分布に近似しないので、検定表などを利用
    • k : 群数
    • n : データ数
    • Ri : 群iの順位平均(タイがある場合はそれらの平均順位を利用)
      • r_i^j : i番目のデータセットで、j番目のアルゴリズムの順位(1位、2位・・・)
      • Ri = 1/n * Σ_i { r_i^j }

χ^2_Fは、nとkが十分大きいとき(n>10, k>5)、k-1自由度のχ^2分布になるとされるが、
F_F = { (n-1)χ^2_F } / { n(k-1) - χ^2_F }がk-1と(k-1)(n-1)自由度のF分布になることも示されており、こちらの方がよりよい統計量であることが示されているっぽい。

平均順位とは、「1.0、0.5、0.5、0.1」のようなスコアが出たとして、タイの順位は2位と3位に該当するので、その平均を取って、「1位、2.5位、2.5位、4位」のようにする。

http://www.snap-tck.com/room04/c01/stat/stat04/stat0403.html
http://www2.hak.hokkyodai.ac.jp/fukuda/lecture/SocialLinguistics/12nonpara.html

post-hoc test (事後検定)

Nemenyi test

2013-10-05

標本抽出メモ

はじめに

大量(または無限)のデータがあっても、人が確認するだとか、1つのデータあたりのなんらかのコストが高い場合、少量のデータを選んで利用する事が多い。(大量に収集されたログデータの分析をするとか、あるプログラムのパフォーマンスを見るために速度を測るとか)
しかし、「少量のデータ」の選び方やその妥当性の判断はかなり難しい。用語などメモしておく。

統計的推測

  • 母集団から標本を抽出して、その標本に対して分析し、母集団について推測すること
  • 統計学や統計的手法を適用する事で、妥当性の判断などが行えるが、適用を間違えると間違った判断を下しかねない
    • 測定量、適用方法、適用手順が適切でないと、母集団と標本が大きく異なってしまう

母集団
  • 対象となるデータすべての集まり

  • 有限母集団 : データ数が有限の場合
    • ログデータ、など
  • 無限母集団 : データ数が無限の場合
    • プログラムの速度測定結果、など

  • パラメトリック : 母集団について、分布の形とそのパラメータがわかるような場合
    • 分布に関連するパラメータや指標を用いて分析などが行える
  • ノンパラメトリック : 母集団について、分布の形やパラメータで分布を決定する事ができないような場合
    • 分布に依らないパラメータを使って分析することになる

標本
  • 母集団から部分的に抽出されたデータ
  • 母集団の縮図的なものと見なされるが、標本における統計量は母集団の統計量とは必ずしも一致しない
    • 一般に、誤差が含まれる
    • そのため、母集団と標本での統計量は区別して扱われる

  • 不完全標本 : ある事象が観測されずに、標本で拾えていない標本(切断、打ち切り)
  • 混合標本 : 2つ以上の事象が混同されて抽出された標本
  • 観測過程で時間/ランダムに変化する標本 : 時系列で変化があるものやランダムに変化する事象から抽出された標本

標本抽出
  • 母集団から標本を選び出すこと
    • 標本の分析結果は「どのような標本抽出をするか」に依存

  • 復元抽出 : 抽出した標本を母集団に戻して、次の抽出対象としても扱う場合
  • 非復元抽出 : 抽出した標本は母集団に戻さない場合

標本サイズの決め方

  • 母集団から何個抽出したらよいか?
  • 基本的に、用途によって異なる
  • 標本サイズが大きくなれば検出力が大きくなる
  • 仮説検定などを行う場合、定義から逆算して求める事ができる

標本抽出方法

単純無作為抽出法
  • 母集団の各要素が標本に含まれる確率を等しいとして、等確率無作為に抽出する方法

確率比例抽出
  • 母集団の各要素に均一でない抽出確率が割り当てられる場合、その確率に基づいて必要数抽出する方法

系統抽出法・等間隔抽出法
  • 無作為に一つ選び、そこからシステマチックに抽出する要素を選んでいく方法
  • 母集団に偏りがある場合に、抽出時にその部分のみを取り出してしまい、偏った標本ができてしまう場合がある

多段抽出法
  • 母集団からの抽出をn段階に分けて無作為抽出を繰り返す方法

層化抽出法
  • 母集団がクラスなどに分割できる場合、それらクラスに分類して、そのクラスごとに無作為抽出する方法
  • 比例割当法は、各クラスから抽出する標本数の大きさを、母集団でのクラスの要素数の割合を使って変える
  • 最適割当法は、母集団でのクラスを標準偏差でわけ、分散の大きいクラスから多く抽出する

標本抽出の問題

信頼性・再現性
  • 同じ手順、条件で抽出した際にいつも同様の結果になるかどうか

妥当性
  • 目的にあった抽出であるかどうか

  • 内的妥当性
    • 抽出した標本が適切な標本だとみなせるかどうか
  • 表面的妥当性
    • 抽出が外見上適切に行われているように見えるか
  • 基準関連妥当性
    • いくつかの抽出での結果で得られる指標とどれだけ一致するか
  • 構成概念妥当性
    • 構成した理論での推定値と標本での結果が一致するかどうか


参考資料

逐次確率比検定を試す

はじめに

あらかじめ標本サイズを決めるのではなく、十分と判断されるまでダイナミックに判断を繰り返す逐次確率比検定を参考に、
チョコボールの銀のエンジェルの出現確率について判断するとどうなるか試してみる。

逐次確率比検定とは

  • ベイズ統計学の枠組みで、ベイズ更新の機能を通して1つずつ標本抽出していきながら同時に検定にも用いる事ができる
  • 逐次決定過程 : 標本抽出をするたびに判断を行い、結論がでたと認められるタイミングで停止する過程
  • 行動
    • action0 : 結論を保留し、標本抽出を再度行う
    • action1 : 帰無仮説H1を採択
    • action2 : 対立仮説H2を採択

尤度比検定(Likelihood Ratio Test)
  • 「尤度比」を検定統計量として行う統計学的検定の総称
    • 尤度比λ=(Π^n_{i=1}{f(Xi|θ1}) / (Π^n_{i=1}{f(Xi|θ2})
    • 帰無仮説H1 : θ=θ1
    • 対率仮説H2 : θ=θ2
  • 尤度比検定 : 任意のC>0に対して、
    • λ<=Cならば、H1を採択
    • λ>Cならば、H2を採択

逐次確率比検定(Sequential Probability Ratio Test)
  • あるA,B>0に対して、n回目まで抽出した標本での尤度比λnに対し、以下のようにする
    • B<λn<A : action0
    • λn <= B : action1
    • A <= λn : action2
    • B = β / (1 - α)
    • A = (1 - β) / α
      • α:第1種の誤りの確率(有意水準)
      • β:第2種の誤りの確率(危険率,1-β=検出力)

二項分布の場合
  • n回の独立試行のうち、k回が事象Tで、(n-k)回が事象Sで、各試行における事象Tである確率pは一定であるような場合、事象Tは二項分布となる
    • P[X=k;p] = nCk * p^k * (1-p)^{n-k}
  • H1 : p = θ1 (出現確率pがθ1より小さい)
  • H2 : p = θ2 (出現確率pがθ2より大きい)
  • 尤度比λn = {(1-θ2)/(1-θ1)}^n * {(θ2/(1-θ2))/(θ1/(1-θ1))}^k
    • β/(1-α) < λn < (1-β)/αである限り抽出を繰り返す


チョコボールのエンジェル

チョコボールのエンジェルとは
  • http://ja.wikipedia.org/wiki/%E3%83%81%E3%83%A7%E3%82%B3%E3%83%9C%E3%83%BC%E3%83%AB
  • チョコボールは森永製菓から発売されているチョコレート菓子
    • 菓子はすべて栃木県の小山市の森永製菓小山工場で生産されている
  • パッケージに「くちばし」と呼ばれる取り出し口がついており、そこに低い確率で「エンジェル」のイラストがついている
    • 金のエンジェル、銀のエンジェルがあり、かつ、金のエンジェルの出現確率<銀のエンジェルの出現確率であることがわかっている
    • 金のエンジェルは1枚で「おもちゃの缶詰」と呼ばれるプレゼントと交換可能
    • 銀のエンジェルは5枚で「おもちゃの缶詰」と交換可能
  • チョコボールの種類として、以下がある
    • ピーナッツ
    • キャラメル
    • いちご
    • クッキー&クリーム
    • 北海道限定夕張メロンチョコボール
    • コンビニ限定商品
  • これまで、このエンジェルの出現確率を求めるため、数々の検証が行われている
チョコボールのエンジェル出現確率の変動
購入判断
  • 銀のエンジェルがでる確率が5%以上ならばチョコボールを買い続けるけど、1%以下ならば買い続けるのはバカらしいのでやめたい、とする
  • 判断方法として、逐次確率比検定の二項分布の場合に照らし合わせて考えてみる
    • 常にエンジェルが入っている確率は、抽出期間内において、時期に依らず一定であると仮定

チョコボールの標本抽出方法
  • チョコボールの標本抽出には以下があると考えられる
    • 母集団を販売されているチョコボールとして、定期的に1つずつ無作為に選んで購入する(無作為抽出)
    • ある店舗において在庫としておいてある箱を複数個無作為に選んで購入し、一つずつ抽出する
    • 店舗、箱、と段階に分けて無作為抽出を繰り返す
    • など

定義
  • 以下と定義して計算する
    • H1 : p = 0.1
    • H2 : p = 0.5
  • チョコボールの値段は安いので、標本サイズが増えたとしても、厳しめに判定したいとする
    • α = 0.05
    • β = 0.005

コード

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

#include <iostream>
#include <vector>
#include <cmath>

static const double PI = 3.141592653589793238;
static const double EPS = 1e-9;

//xorshift
// 注意: longではなくint(32bit)にすべき
unsigned long xor128(){
  static unsigned long x=123456789, y=362436069, z=521288629, w=88675123;
  unsigned long t;
  t=(x^(x<<11));
  x=y; y=z; z=w;
  return w=(w^(w>>19))^(t^(t>>8));
}
//[0,1)の一様乱数
// 注意: int_maxぐらいで割った方がよい
double frand(){
  return xor128()%10000000000/static_cast<double>(10000000000); 
}


bool getChocoBall(double p){
  double r = frand();
  if(r < p) return true;
  return false;
}


int main(){
  int n = 0, k = 0;
  double p = 0.005; //実際の出現確率
  
  double theta1 = 0.01; //H1
  double theta2 = 0.05; //H2
  double alpha = 0.05; //有意水準
  double beta = 0.005; //危険度

  while(true){
    n++;

    //購入
    bool res = getChocoBall(p);
    if(res) k++;

    //対数尤度比の計算
    double cal = 0.0;
    cal += n * log( (1.0-theta2)/(1.0-theta1) );
    cal += k * log( (theta2/(1.0-theta2)) / (theta1/(1.0-theta1)) );

    //判定
    double B = log( beta / (1.0-alpha) );
    double A = log( (1.0-beta) / alpha );
    //std::cout << B << "\t" << cal << "\t" << A << std::endl;
    if(cal+EPS < B){
      std::cout << "STOP" << std::endl; //買うのをやめたい
      break;
    }
    else if(A+EPS < cal){
      std::cout << "CONTINUE" << std::endl; //これからも買い続けよう
      break;
    }

    std::cout << n << "(" << k << ") : MEASURING..." << std::endl;
  }
  return 0;

}

結果

1000回ずつ、買い続けるか(CONTINUE)、買うのをやめるか(STOP)を判定したものを比較してみる。



実際の確率STOPになった個数CONTINUEになった個数抽出回数についての標本平均抽出回数についての標本分散
p=0.00110000179.05428.654
p=0.0059982220.8043607.37
p=0.0196634302.27619772.7
p=0.02537463449.92128463
p=0.0363937255.58652570.5
p=0.045995130.00912308.3
p=0.05199987.4345358.39
p=0.100100032.315436.67

参考

2013-09-28

Mann-WhitneyのU検定メモ

はじめに

ROC曲線のAUCとも関係のあるらしいマン・ホイットニーのU検定についてメモ。

マン・ホイットニーのU検定とは

  • Mann-Whitney's U test
    • ウィルコクソンの順序和検定も同等のため、マン・ホイットニー・ウィルコクソン(MWW)検定とも
  • ノンパラメトリックな統計的検定の一種
    • 母集団のパラメータを仮定しないdistribution-freeな方法
    • 一般に、標本数が少ない場合は帰無仮説が偽で、帰無仮説が棄却される可能性が低くなる
  • ある2つの(対応関係のない)標本(群)の母代表値に差があるかどうか」を検定
    • 帰無仮説:両群に差はない(同じ分布をしている)
    • 対立仮説:両群に差がある

統計量U

Uの定義はいくつか種類があるよう。

標本数が少ないとき
  • 標本数が小さい方の群をA、他方の群をBとすると、Aの標本iについて、Bの中で標本iより小さい値の個数を求める
  • すべてのiについて個数の総和をUとする

標本数が少し多めのとき
  • 標本を順位に直し、それぞれの標本群についてその順位の和を求め、その小さい方をUとする
  • 統計量U = min(U1, U2)
    • U1 = n1*n2 + ( n1*(n1+1) )/2 - R1
    • U2 = n1*n2 + ( n2*(n2+1) )/2 - R2
    • n1,n2 : 群1、群2の標本数
    • R1,R2 : 群1、群2の標本の順位の和
  • 同順位がある場合はこの方法は使えないので、t検定(ウェルチのt検定)/z検定を用いる
  • データ数は20程度まで。それ以上あるならば、正規分布で近似できるので、z検定する

検定
  • Uが「棄却限界値以下の時」に帰無仮説が棄却される
    • U<=棄却限界値なら帰無仮説が棄却される = 両群(の平均順位)に差がある
    • U>棄却限界値なら帰無仮説は棄却されない = 両群(の平均順位)に差があるとはいえない
  • 棄却限界値の数表

z検定する場合

  • 順位に対してz検定する
  • 統計量z = | U-E(U) | / sqrt(V(U))
    • 平均 : E(U) = n1*n2/2
    • 分散 : V(U) = n1*n2*(n1+n2+1)/12
      • ただし、同順位がある場合は、V(U)=n1*n2/(12*(n^2-n)) * {n^3-n-Σ_{i=1}^{m}{t^3_i-t_i}}を用いる
      • m : 同順位の種類数
      • t : 同順位の個数
  • z値からp値を求め、有意水準以下の場合、帰無仮説を棄却 = 両群に差がある

ROC曲線AUCと統計量U

AUCと統計量Uが関係するとのことで、(Continuous forecast probabilities, no tiesだけ)ちょっと以下の論文を追ってみた。

S.J.Manson and N.E.Graham, Areas beneath the relative operating characteristics (ROC) ans relative operating levels (ROL) curves: Statistical significance and interpretation
http://www.inmet.gov.br/documentos/cursoI_INMET_IRI/Climate_Information_Course/References/Mason+Graham_2002.pdf

  • ROC空間に点が打たれて、それらをつないだものがROC曲線
    • その曲線の下側の面積がAUC
  • AUCを計算するために、ROC曲線の隣り合う点間の面積を計算すると、
    • 真陽率(y軸)方向では、TPの数が1つ分あがることになるので、1/(TP+FN)
      • TP+FN=eとしておく。1/e。
    • 偽陽率(x軸)方向では、もし偽陽率があがる場合、(TN+FP)中f上がるとすると、f/(TN+FP)
      • TN+FP=e'としておく。f/e'。また、e+e'=TP+FN+TN+FP=n。
    • 下側面積なので、点間あたり、(1/e) * (1-f_i/e') = (e'-f_i)/(e*e')だけ面積が増加する
  • なので、AUC=Σ^{e}_{i=1}{(e'-f_i)/(e*e')}となる
    • すこし変形し、AUC=1-1/(e*e')*Σ^{e}_{i=1}{f_i}
    • ここで、Σf_iについて、順位差を導入してΣ(r_i-i)と書くと、Σ^{e}_{i=1}{f_i}=Σ^{e}_{i=1}{r_i} - e*(e+1)/2
  • U統計量=Σ^{e}_{i=1}{r_i}-e*(e+1)/2と書けることから、
  • AUC=1-1/(e*e')*U=(e*e'-U)/(e*e')、または、U=(e*e')*(1-AUC)の関係がある

参考資料