歪んだサイコロでベイズ - Negative/Positive Thinking

はじめに

歪んだサイコロを使ってベイズ統計学の初歩を勉強してみる。

歪んだサイコロ

6面サイコロがある。しかし、よく見ると各面の大きさが違う。
以下、疑似的に出目を生成してみる。

歪んだサイコロの出目

コード

#include <iostream>
#include <vector>
#include <climits>

static const int NUM_SIDE = 6; //サイコロの面の数
static const double prob[NUM_SIDE] = {6.0/21, 5.0/21, 4.0/21, 3.0/21, 2.0/21, 1.0/21}; //各サイコロ面の出現確率
//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()%ULONG_MAX/static_cast<double>(ULONG_MAX); 
}
//サイコロを振る
int roll_dice(){
  int ret = 0;
  double r = frand(), sum = 0;
  for(int i=0; i<NUM_SIDE; i++){
    sum += prob[i];
    if(r < sum) return i+1;
  }
  return NUM_SIDE;
}

int main(){
  
  int num[NUM_SIDE+1];
  for(int i=0; i<NUM_SIDE+1; i++) num[i] = 0;
  
  for(int t=0; t<100000; t++){
    int a = roll_dice();
    if(t<10) std::cout << a << " ";
    num[a]++;
  }
  std::cout << std::endl;

  for(int i=1; i<NUM_SIDE+1; i++) std::cout << num[i] << std::endl;
  
  return 0;
}

サイコロを振ったときの出目

5 1 3 4 1 3 3 1 1 1 ...

出目の回数

#| 1 | 2 | 3 | 4 | 5 | 6 |
# 100回
| 34 | 20 | 19 | 16 | 7 | 4 |
# 100000回
| 28669 | 23738 | 19074 | 14164 | 9510 | 4845 |

ベイズ推定

サイコロの目の出る確率はわからないものとして、この確率をベイズ推定してみる。

尤度関数

N回の独立したベルヌーイ試行の結果、K個の値のどれかをとり、それぞれの値が出る確率はμ_1,...,μ_K (Σμ_i = 1)のとき、N回試行後、値iが出る回数X_iは「多項分布」に従う。
よって、N回サイコロを振って、各出目の回数は「多項分布」に従う。

多項分布
$p ( x | \mu ) = \prod_{k=1}^{M}{\mu_{k}^{x_k}}$

自然共役分布

「多項分布」の自然共役分布は「ディリクレ分布」なので、これを事前分布にする。すると、事後分布もディリクレ分布となる。

ディリクレ分布
$Dir ( \mu | \alpha ) = \frac{\Gamma ( \sum_{k=1}^{K}{\alpha_k} )}{\Gamma(\alpha_1)...\Gamma(\alpha_K)}{\prod_{k=1}^{K}{\mu_{k}^{\alpha_k -1}}}$

ディリクレ分布グラフ(3変数上)
- コード

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

//ディリクレ分布Dir(alpha, p)
double dirichlet_P(const std::vector<double>& alpha, const std::vector<double>& p){
  double ret = 0.0;
  double sum = 0.0;
  for(int i=0; i<alpha.size(); i++) sum += alpha[i];
  
  ret += lgamma(sum);
  for(int i=0; i<alpha.size(); i++) ret -= lgamma(alpha[i]);
  for(int i=0; i<alpha.size(); i++) ret += (alpha[i]-1) * log(p[i]);
  
  return exp(ret);
}

int main(){
  std::vector<double> alpha, p;
  for(int i=0; i<3; i++){
    alpha.push_back(0.0);
    p.push_back(0.0);
  }
  //alpha
  alpha[0] = 0.1; //alpha_1
  alpha[1] = 0.1; //alpha_2
  alpha[2] = 0.1; //alpha_3
  
  for(double p1=0.01; p1<1.0; p1+=0.01){
    for(double p2=0.01; p2<1.0; p2+=0.01){
      if(p1+p2>1.0) continue;
      p[0] = p1;
      p[1] = p2;
      p[2] = 1.0-p1-p2;
      std::cout << p1 << "\t" << p2 << "\t" << dirichlet_P(alpha, p) << std::endl;
    }
  }
  
  return 0;
}

- α=(0.1, 0.1, 0.1)

- α=(1.0, 1.0, 1.0)

- α=(10, 10, 10)

事後分布

「事後分布∝事前分布*尤度関数」なので、ディリクレ分布と多項式を掛け合わせると、事後分布は、
$p( \mu | D, \alpha ) = \frac{\Gamma ( M + \sum_{k=1}^{K}{\alpha_k} )}{\Gamma(\alpha_1 + m_1)...\Gamma(\alpha_K + m_K)}{\prod_{k=1}^{K}{\mu_{k}^{\alpha_k +m_k -1}}}$
したがって、
$p( \mu | D, \alpha ) = Dir(\mu |\alpha + m)$
となる。

ということで、データとしてサイコロの出目が与えられると、対応するディリクレ分布のパラメータが「α_i→α_i + 1」とベイズ更新される。

分布の代表値

パラメータの分布がわかっても、具体的に「目iの出現確率はp_i」というところまでは言っていない。
推定誤り損失の定義によって、分布の代表値として「中央値」「平均値」「最頻値」「パーセンタイル」などがあるが、ここではディリクレ分布の「平均値(期待値)」「最頻値(モード)」を見てみる。

ディリクレ分布の平均値
$E[\mu_i] = \frac{\alpha_i}{\sum {\alpha_k}}$

ディリクレ分布の最頻値
$mode[\mu_i] = \frac{\alpha_i - 1}{\{\sum {\alpha_k}\} - K}$

コード

#include <iostream>
#include <vector>
#include <climits>

static const int NUM_SIDE = 6; //サイコロの面の数
static const double prob[NUM_SIDE] = {6.0/21, 5.0/21, 4.0/21, 3.0/21, 2.0/21, 1.0/21}; //各サイコロ面の出現確率
//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()%ULONG_MAX/static_cast<double>(ULONG_MAX); 
}
//サイコロを振る
int roll_dice(){
  int ret = 0;
  double r = frand(), sum = 0;
  for(int i=0; i<NUM_SIDE; i++){
    sum += prob[i];
    if(r < sum) return i+1;
  }
  return NUM_SIDE;
}

int main(){
  
  int alpha[NUM_SIDE+1];
  for(int i=0; i<NUM_SIDE+1; i++) alpha[i] = 1;

  int N = 0;
  for(int t=0; t<5000; t++){
    int a = roll_dice();
    alpha[a]++;
    N++;
    //平均値
    for(int i=1; i<NUM_SIDE+1; i++) std::cout << alpha[i]/static_cast<double>(NUM_SIDE+N) << "\t";
    //最頻値
    //for(int i=1; i<NUM_SIDE+1; i++) std::cout << (alpha[i]-1)/static_cast<double>((NUM_SIDE+N)-NUM_SIDE) << "\t";
    std::cout << std::endl;
  }
  
  return 0;
}