Negative/Positive Thinking

ベータ分布の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に近い点だと誤差が大きくなっている模様)

参考