逐次確率比検定を試す

確率 統計

はじめに

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

逐次確率比検定とは

  • ベイズ統計学の枠組みで、ベイズ更新の機能を通して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-β)/αである限り抽出を繰り返す

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

チョコボールのエンジェルとは
チョコボールのエンジェル出現確率の変動
購入判断
  • 銀のエンジェルがでる確率が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.001 1000 0 179.05 428.654
p=0.005 998 2 220.804 3607.37
p=0.01 966 34 302.276 19772.7
p=0.02 537 463 449.92 128463
p=0.03 63 937 255.586 52570.5
p=0.04 5 995 130.009 12308.3
p=0.05 1 999 87.434 5358.39
p=0.10 0 1000 32.315 436.67

参考