ギブスサンプリングを試す(2次元正規分布の場合)

プログラミング 確率

はじめに

ギブスサンプリングで2次元正規分布からサンプリングを試してみる。

サンプリングの流れ

2次元正規分布なので(x,y)がサンプリングされる。

  • 初期値(x0,y0)を決める
  • y0を固定した場合の条件付き正規分布からx1をサンプリングする→(x1,y0)
  • 次に、x1を固定した場合の条件付き正規分布からy1をサンプリングする→(x1,y1)
  • 繰り返し
  • 1次元の(条件付き)正規分布からのサンプリングはBox-Muller法で行ってみる

条件付き正規分布

  • 多次元正規分布は、いくつかの変数集合が与えられた場合の条件付き分布も正規分布になる性質がある(PRML第2章)
  • 2次元正規分布のとき、(x,y)のうちyが固定された場合、xは以下の1次元正規分布に従う
    • 平均：\mu_{x|y} = \mu_{x} + \Sigma_{xy} \Sigma_{yy}^{-1} (y-\mu_y)
    • 分散：\Sigma_{x|y} = \Sigma_{xx} - \Sigma_{xy} \Sigma_{yy}^{-1} \Sigma_{yx}

コード

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

//2次元正規乱数(Gibbsサンプリング)
class normal_rand_gibbs {
  static const double PI = 3.14159265358979323846264338;
  //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/10000000000.0;
    //上の返り値は偏りがありそうなので、以下を用いた方がよさそう
    //return xor128()%ULONG_MAX/static_cast<double>(ULONG_MAX); 
  }
  //1次元正規乱数(Box-Muller法)
  double normal_rand(double mu, double sigma2){
    double sigma = sqrt(sigma2);
    double u1 = frand(), u2 = frand();
    double z1 = sqrt(-2*log(u1)) * cos(2*PI*u2);
    //double z2 = sqrt(-2*log(u1)) * sin(2*PI*u2);
    return mu + sigma*z1;
  }
public:
  normal_rand_gibbs(double x0, 
		    double y0, 
		    double mu_x, 
		    double mu_y,
		    double sig_xx,
		    double sig_xy,
		    double sig_yy
		    ):
    tx(x0),
    ty(y0),
    mu_x(mu_x),
    mu_y(mu_y),
    sig_xx(sig_xx),
    sig_xy(sig_xy),
    sig_yy(sig_yy),
    turn(false){}

  //2次元正規乱数生成
  std::pair<double, double> rand(){
    if(!turn){ //xをサンプリング
      tx = normal_rand(mu_x + (sig_xy/sig_yy)*(ty-mu_y),
		       sig_xx + sig_xy*sig_xy/sig_yy
		       );
    }else{ //yをサンプリング
      ty = normal_rand(mu_y + (sig_xy/sig_xx)*(tx-mu_x),
		       sig_yy + sig_xy*sig_xy/sig_xx
		       );
    }
    
    if(turn) turn = false;
    else turn = true;

    return std::make_pair(tx,ty);
  }
  
private:
  double tx, ty; //現在の位置
  bool turn; //false->x, true->y
  double mu_x, mu_y; //平均
  double sig_xx, sig_xy, sig_yy; //分散
};

int main(){
  // 平均 = ( mu_x, mu_y )
  const double mu_x = 2.0;
  const double mu_y = 3.0;
  // 共分散行列 = ( sig_xx  sig_xy )
  //              ( sig_xy  sig_yy )
  const double sig_xx = 5.0;
  const double sig_xy = -1.5;
  const double sig_yy = 3.0;
  //スタート位置
  const double start_x = 0.0;
  const double start_y = 0.0;


  normal_rand_gibbs nrg(start_x, start_y, mu_x, mu_y, sig_xx, sig_xy, sig_yy);

  std::cout << start_x << "\t" << start_y << std::endl;
  for(int i=2; i<=3000; i++){
    std::pair<double,double> pos = nrg.rand();
    std::cout << pos.first << "\t" << pos.second << std::endl;
    //if(i%100==0) std::cout << "\n\n" << std::endl;
  }

  return 0;
}

結果

  • ケース1
    • 平均：(0, 0)
    • 分散：( (1.0, 0.5), (0.5, 1.0) )
    • 初期位置：(10.0, 10.0)
    • 100ターンごとに色を変えて


  • ケース2
    • 平均：(2.0, 3.0)
    • 分散：( (5.0, -1.5), (-1.5, 3.0) )
    • 初期位置：(0.0, 0.0)
    • 3000回分

一応、分布に従ってるよう。
例えば、サンプリングした全ての点を使って平均や分散を計算すれば指定の分布のパラメータの近似になる(はず)
ケース1のように、初期値によっては適切なサンプリングに収束するまである程度かかる(burn-in)。
なので、最初の何ターンかは捨てる必要がある。

参考文献