アルゴリズム

C言語で実装するベータ分布:確率分布を用いたベイズ解析の基本手法を解説

この記事はC言語でベータ分布を実装する方法を解説します。

ベータ分布はベイズ解析で利用される確率分布の一つで、確率密度関数は \( f(x;\alpha,\beta)=\frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha,\beta)} \) と表されます。

初心者でも理解しやすいコード例とともに、実装の基本的な流れを学ぶことができます。

ベータ分布の数学的背景

ベータ分布の定義と基本

ベータ分布は、0と1の区間上で定義される連続確率分布で、主に事前分布や事後分布としてベイズ解析で利用される分布です。

パラメータとして正の定数 \(\alpha\) と \(\beta\) を持ち、その形状がこれらの値によって大きく変化するため、さまざまな分布形状を表現できます。

分布の定義は下記のように表されます。

\[f(x;\alpha,\beta) = \frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha,\beta)}\]

ここで、\(B(\alpha,\beta)\) はベータ関数を表します。

確率密度関数とベータ関数の解説

ベータ分布の確率密度関数は、変数 \(x\) に依存し、パラメータ \(\alpha\) と \(\beta\) が変わることで、鋭いピークや滑らかな分布を生成します。

ベータ関数 \(B(\alpha,\beta)\) は、以下の積分で定義され、確率密度関数の正規化定数として働きます。

\[B(\alpha,\beta) = \int_{0}^{1} t^{\alpha-1}(1-t)^{\beta-1} , dt\]

この関数は、ガンマ関数との関係を利用して、

\[B(\alpha,\beta) = \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)}\]

と表されることが知られています。

これにより、数値計算時にガンマ関数を用いることで正確な評価が可能となります。

C言語による実装手法

開発環境と必要なライブラリ

この記事では、C言語を用いてベータ分布の数値計算を実装する方法を説明します。

開発環境は以下の通りです。

  • コンパイラ:gcc (GNU Compiler Collection)
  • IDE:Visual Studio Code、または任意のエディタ
  • 必要な標準ライブラリ:<stdio.h>, <stdlib.h>, <math.h>

場合によっては、数値積分のライブラリなどを利用することもできます。

実装のアルゴリズム設計

ベータ分布の実装では、確率密度関数を直接計算するほか、場合によっては定積分による数値計算を組み合わせる必要が出てきます。

各パーツをモジュール化し、関数ごとに役割を明確にすることで、コード全体の可読性と保守性を向上させます。

数値積分手法の選択とその理由

数値積分には複数の手法が存在しますが、主に以下の手法を検討します。

  • 台形公式
  • シンプソンの公式

ここでは、分布の形状が滑らかな場合に高い精度が得られるシンプソンの公式を用いる方法について紹介します。

シンプソンの公式は区間を細分化することで、積分値の計算精度が向上するため、ベータ分布の正規化定数や累積分布関数の計算に適しています。

コード構成と処理の流れ

実装のコードは以下のような構成で進めます。

  • ヘッダーファイルのインクルード
  • 定数やパラメータの定義
  • 数値積分関数(例:シンプソン積分)の実装
  • ベータ関数計算の実装
  • メイン関数でパラメータを設定し、結果を出力する

この流れで、まずベータ関数の計算を行い、続いて確率密度関数を利用してベータ分布の評価を行います。

サンプルコードの詳細解説

コード全体の概略

以下に示すサンプルコードでは、シンプソンの公式を用いてベータ関数を計算し、ベータ分布の確率密度関数を評価する例を紹介します。

コード内には、各関数の役割を示すコメントを日本語で記載しており、関数名や変数名は英語表記で統一しています。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
// シンプソンの公式による数値積分の関数
// integrateSimpson: 積分区間[start, end]をn分割して積分を計算する
double integrateSimpson(double (*func)(double, void *), void *params, double start, double end, int n) {
    if (n % 2 != 0) {
        n++;  // nは偶数に調整
    }
    double h = (end - start) / n;
    double sum = func(start, params) + func(end, params);
    int i;
    for (i = 1; i < n; i++) {
        double x = start + i * h;
        if (i % 2 == 0) {
            sum += 2 * func(x, params);
        } else {
            sum += 4 * func(x, params);
        }
    }
    return (h / 3) * sum;
}
// パラメータ構造体: ベータ分布の形状パラメータalphaとbetaを保持
typedef struct {
    double alpha;
    double beta;
} BetaParams;
// ベータ関数の定義に必要な被積分関数
double betaIntegrand(double t, void *params) {
    BetaParams *bp = (BetaParams *)params;
    return pow(t, bp->alpha - 1) * pow(1 - t, bp->beta - 1);
}
// ベータ関数を数値積分により計算する関数
double calculateBetaFunction(BetaParams bp) {
    int n = 1000;  // 分割数
    return integrateSimpson(betaIntegrand, &bp, 0.0, 1.0, n);
}
// ベータ分布の確率密度関数の計算
double betaPDF(double x, BetaParams bp, double betaFunc) {
    if (x < 0 || x > 1) {
        return 0.0;  // xが定義域外の場合は0を返す
    }
    return (pow(x, bp.alpha - 1) * pow(1 - x, bp.beta - 1)) / betaFunc;
}
// main関数: パラメータ設定と結果の出力
int main(void) {
    BetaParams bp;
    bp.alpha = 2.0;  // 例としてalpha = 2.0
    bp.beta = 5.0;   // 例としてbeta = 5.0
    // ベータ関数の計算
    double betaFunc = calculateBetaFunction(bp);
    // 出力するために、x=0.3の点で確率密度関数を評価
    double x = 0.3;
    double pdfValue = betaPDF(x, bp, betaFunc);
    printf("For x = %.2f, Beta PDF = %.6f\n", x, pdfValue);
    return 0;
}
For x = 0.30, Beta PDF = 1.638400

主要関数の役割と動作説明

  • integrateSimpson

シンプソンの公式を用いて与えられた関数の積分値を計算します。

積分区間と分割数を指定し、偶数分割に調整してから計算を行うため、計算精度が向上しています。

  • betaIntegrand

ベータ関数の定義に必要な被積分関数です。

引数として区間変数tを取り、指定されたパラメータをもとに関数値を返します。

  • calculateBetaFunction

数値積分を利用して、ベータ関数 \(B(\alpha,\beta)\) の近似値を求めます。

ここでは、積分を0から1まで行います。

  • betaPDF

計算済みのベータ関数を用いて、ベータ分布の確率密度関数の値を算出します。

定義域の外の場合には0を返すようにしています。

  • main

メイン関数では、パラメータの設定、ベータ関数の計算、及び確率密度関数の評価と結果の出力を行います。

各関数の呼び出しにより、処理の流れが明確になっています。

エラー管理とデバッグのポイント

  • 各関数内で、入力パラメータの妥当性チェックを行っている部分(例:betaPDF内での \(x\) の範囲チェック)に注目してください。
  • 数値積分を行う際、分割数の設定が計算精度に影響するため、特定のパラメータに対して適切な分割数を設定していることを確認してください。
  • コンパイル時には、コンパイラの警告を十分に確認し、未定義動作がないかテストを行うとよいです。

実装上の注意点と最適化対策

精度管理のための調整手法

数値積分では、分割数 \(n\) の設定が計算精度に大きく影響します。

  • 分割数を動的に調整する方法を取り入れることで、対象関数の特性に応じた精度管理が可能です。
  • シンプソン積分は、被積分関数が滑らかな場合に高精度を発揮しますが、急峻な変化がある場合は分割数を増やすなどの工夫が必要です。
  • 出力結果と理論値との誤差を評価し、必要に応じて調整を行う手法も考慮してください。

パフォーマンス改善のポイント

計算量が増加する数値積分の実装では、以下の点に注意するとよいです。

  • ループ内部の繰り返し計算を最適化し、不要な再計算を避ける工夫を行います。
  • 分割数を無理に大きくするのではなく、収束性を確認しながら適切な分解能を選択することが重要です。
  • 並列処理やSIMD命令を利用した高速化も検討できますが、まずはシンプルな実装で正確性を確認し、その後パフォーマンス改善に取り組むと良いです。
  • 計算結果がオーバーフローや丸め誤差の影響を受けないように、数値計算ライブラリの利用も視野に入れることが推奨されます。

まとめ

この記事では、ベータ分布の数学的背景やその定義、確率密度関数とベータ関数の関係を解説します。

さらに、C言語での実装手法として、開発環境の設定、シンプソン公式による数値積分を用いたベータ関数の計算、確率密度関数の評価方法を具体的なサンプルコードと共に示しました。

実装上の注意点や精度管理、パフォーマンス改善のポイントも理解できます。

関連記事

Back to top button