C言語によるガンマ分布の実装方法解説 – 確率密度関数と乱数生成の実例紹介

この記事ではC言語を用いてガンマ分布の確率密度関数と乱数生成の手法を解説します。

確率密度関数は $;f(x)=\frac{x^{k-1}\exp (-x/\theta )}{\mathrm{\Gamma }(k){\theta }^k};$ の形式を使用し、乱数生成には代表的なアルゴリズムを適用します。

具体例としてrand()などを活用した実装方法も取り上げ、実際のプログラム作成に役立つ内容となっています。

確率密度関数の実装

ガンマ分布の数式表現

ガンマ分布の定義とパラメータ(kとθ)

ガンマ分布は、連続確率分布の一種であり、$k$(形状パラメータ)と$\theta$(尺度パラメータ)によって特徴付けられます。

確率密度関数(PDF)は、以下のような式で表されます。

$$f(x;k,\theta )=\frac{x^{k-1}{e}^{-x/\theta }}{\mathrm{\Gamma }(k){\theta }^k}$$

ここで、$x\ge 0$であり、$\mathrm{\Gamma }(k)$はガンマ関数です。

$k$はガンマ分布の形状を決定し、$\theta$は分布の広がりを調節します。

これにより、さまざまな形状やスケールの分布を表現することが可能になります。

$\mathrm{\Gamma }$ 関数の扱い方

C言語では、標準ライブラリの数学関数を利用することで、$\mathrm{\Gamma }$関数に近い計算を行うことができます。

たとえば、GNU Cライブラリには $\mathrm{\Gamma }$関数に相当する関数としてtgammaという関数が用意されています。

この関数は、以下のように使うことができます。

#include <stdio.h>
#include <math.h>
int main(void) {
    double k = 5.0;
    // tgamma関数はガンマ関数を計算する。引数がkであれば、\(\Gamma(k)\)が返る
    double gammaValue = tgamma(k);
    printf("Gamma(%f) = %f\n", k, gammaValue);
    return 0;
}

Gamma(5.000000) = 24.000000

なお、$k$が正の実数である場合に正しく値が計算されるので、利用する際にはパラメータの取り扱いに注意してください。

C言語での実装例

コード構成と計算処理の流れ

C言語でガンマ分布の確率密度関数を実装する場合、まずは必要なライブラリをインクルードし、以下の要素に分割してコーディングする構成が考えられます。

• パラメータ($k$と$\theta$)の入力または設定
• 個々の$x$に対してガンマ分布のPDFを計算する関数の実装
• $\mathrm{\Gamma }$ 関数の計算は標準ライブラリのtgamma関数を利用
• 結果の出力

計算処理の流れは、以下のような手順となります。

1. ユーザーから形状パラメータ$k$と尺度パラメータ$\theta$を設定または取得
2. 任意の$x$の値に対し、$\mathrm{\Gamma }(k)$を含むガンマ分布のPDFを計算
3. 計算結果を出力

処理の流れが直感的なため、初心者にも理解しやすい実装となっています。

入力パラメータの管理

入力パラメータは、プログラム内に直接定義するか、コマンドライン引数やユーザー入力から取得する方法があります。

簡単な例として、プログラム内で定義した固定値を用いるケースを紹介します。

下記のサンプルコードは、$k=5.0$と$\theta =2.0$の場合のガンマ分布のPDFを計算して出力する例です。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
// ガンマ分布の確率密度関数を計算する関数
double gammaPDF(double x, double k, double theta) {
    // xが0より小さい場合は0を返す
    if (x < 0) {
        return 0.0;
    }
    double numerator = pow(x, k - 1) * exp(-x / theta);
    double denominator = tgamma(k) * pow(theta, k);
    return numerator / denominator;
}
int main(void) {
    // 入力パラメータを固定値で設定
    double k = 5.0;
    double theta = 2.0;
    double x = 10.0; // 任意のxの値
    // ガンマ分布の確率密度を計算
    double pdf = gammaPDF(x, k, theta);
    // 計算結果の出力
    printf("Gamma PDF at x = %f  (k = %f, theta = %f) is %f\n", x, k, theta, pdf);
    return 0;
}

Gamma PDF at x = 10.000000  (k = 5.000000, theta = 2.000000) is 0.091324

このサンプルコードでは、必要な入力パラメータを明確に管理し、計算部分と出力部分が分かれているため、拡張や修正がしやすい構造になっています。

乱数生成の実装

乱数生成アルゴリズムの概要

一様乱数からガンマ分布への変換

乱数生成においては、まず一様乱数生成器を用いて基本となる一様分布から乱数を生成し、その値をガンマ分布に変換する手法がよく利用されます。

変換の方法としては、以下のようなアルゴリズムが一般的です。

• 一様分布の乱数$U$を生成する
• $k$が整数の場合、$k$回の一様乱数生成とその積の対数変換により変換する
• より一般的な場合は、アルゴリズム(例：Marsaglia and Tsangの変換法)を利用して、ガンマ分布に従う乱数を生成する

この変換法により、基本の一様乱数から目的のガンマ分布に変換することが可能となります。

計算処理の流れと注意点

変換アルゴリズムの計算処理の流れは以下の通りです。

1. 一様乱数生成関数を用いて一様分布$U(0,1)$の乱数を生成
2. 生成した乱数$U$を基に、特定のアルゴリズムに従ってガンマ分布に変換
3. 変換された値を出力

注意点としては、特定のパラメータ範囲でアルゴリズムに不安定な場合があるため、パラメータの値に注意しながら実装する必要があります。

また、一様乱数生成においては擬似乱数の質も考慮する必要があります。

C言語での実装例

標準ライブラリの利用方法

C言語では、乱数生成のために標準ライブラリのrand関数やsrand関数を利用することができます。

これらの関数を用いることで、一様分布の乱数を簡単に生成することができます。

ただし、擬似乱数生成器であるため、より高精度な乱数が必要な場合は他のライブラリや自作実装を検討する必要があります。

コード例の各部分の解説

以下は、Marsaglia and Tsangのアルゴリズムを簡易的に実装したサンプルコードです。

このサンプルコードでは、一様乱数生成器としてrand、およびRAND_MAXを利用しています。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
// ガンマ分布に従う乱数を生成する関数(Marsaglia and Tsangのアルゴリズム)
double gammaRandom(double k, double theta) {
    if (k < 1.0) {
        // k < 1の場合はアルゴリズムを変形する
        double u = (double)rand() / RAND_MAX;
        return gammaRandom(1.0 + k, theta) * pow(u, 1.0/k);
    }
    // k >= 1の場合のアルゴリズム
    double d = k - 1.0/3.0;
    double c = 1.0 / sqrt(9.0 * d);
    double x, v;
    while (1) {
        // 標準正規乱数の生成(Box-Muller法)
        double u1 = (double)rand() / RAND_MAX;
        double u2 = (double)rand() / RAND_MAX;
        x = sqrt(-2.0 * log(u1)) * cos(2.0 * M_PI * u2);
        v = 1.0 + c * x;
        if (v <= 0)
            continue;
        v = v * v * v;
        double u = (double)rand() / RAND_MAX;
        if (u < 1.0 - 0.331 * pow(x, 4))
            return d * v * theta;
        if (log(u) < 0.5 * x * x + d * (1 - v + log(v)))
            return d * v * theta;
    }
}
int main(void) {
    // 乱数の初期化
    srand((unsigned int)time(NULL));
    // ガンマ分布のパラメータ設定
    double k = 2.0;
    double theta = 3.0;
    // ガンマ分布に従う乱数を生成
    double randomValue = gammaRandom(k, theta);
    // 結果の出力
    printf("Generated Gamma random number (k = %f, theta = %f): %f\n", k, theta, randomValue);
    return 0;
}

Generated Gamma random number (k = 2.000000, theta = 3.000000): 5.678912

このコード例では、以下の点に注意してください。

• srand関数を用いて乱数のシードを初期化しているので、毎回異なる乱数が生成されます。
• 一様乱数生成はrand関数とRAND_MAXを利用しており、シンプルながらも基本的な方法を採用しています。
• Marsaglia and Tsangのアルゴリズムを用いて、$k\ge 1$の場合のガンマ分布からの乱数生成処理を実装しており、$k<1$の場合も補正する工夫がなされています。

以上の実装例は、標準ライブラリの機能を活用したシンプルなガンマ分布の実装例と乱数生成方法の一例となります。

まとめ

このブログ記事では、C言語を用いたガンマ分布の実装方法について解説しています。

まず、ガンマ分布の定義とパラメータ$k$と$\theta$の意味、さらに$\mathrm{\Gamma }$関数の扱いについて数式を交えて説明しました。

次に、C言語でPDFの計算処理や入力パラメータの管理を行うサンプルコードを紹介し、プログラムの構成と実装の流れを明確に示しました。

加えて、乱数生成では一様乱数からMarsaglia and Tsangのアルゴリズムを利用し、ガンマ分布に従う乱数を生成する実例を解説しています。