【C++】ポアソン分布を使用したスペクトルデータのシミュレーション

2026-08-25C++,スペクトル,数学

1. はじめに
2. ポアソン分布とは
- 2.1. 基本式
  - 2.1.1. パラメータ$\lambda$の意味
- 2.2. ポアソン分布の主要な統計量
3. C++におけるポアソン分布の使用方法
- 3.1. コンストラクタとパラメータ
- 3.2. パラメータを後から変更する
4. スペクトルデータのシミュレーション
- 4.1. 実装例
- 4.2. 実装上の注意点
  - 4.2.1. $\lambda$が0またはそれに近い場合
5. まとめ

はじめに

物理計測やセンサー技術の分野において、観測されるデータは理論上の数値（真の値）に統計的な揺らぎ（ノイズ）が乗ったものになります。とくに、光子計数や放射線計測のように、単位時間あたりのイベント発生をカウントするプロセスでは、そのノイズは「ポアソン分布」に従うことが知られています。

本エントリでは、C++の<random>ライブラリを用いて、理論的なモデル（背景放射線やピーク構造を持つ関数）から、いかにも「実測データらしい」ノイズを含んだスペクトルデータを生成するシミュレーション手法について説明します。

ポアソン分布とは

ある事象が一定の平均率で独立に発生すると仮定したとき、その発生回数が従う分布をポアソン分布と言います。

たとえば：

1時間あたりの平均着信件数 $\lambda = 5$のとき、ちょうど3件着信する確率はどれくらいか？
1日あたりの平均事故件数 $\lambda = 2$のとき、0件の確率はどれくらいか？
1メートルあたりの平均欠陥点数 $\lambda = 0.7$のとき、2点以上の欠陥がある確率はどれくらいか？

このような「単位時間・単位面積・単位体積あたりの平均発生率」をパラメータとして持つ分布です。

基本式

$X$がちょうど$k$回発生する確率は、

$$
P(X = k) = \frac{e^{-\lambda} \lambda^k}{k!}, \quad k = 0, 1, 2, \dots
$$

です。

ここで：

$e$ は自然対数の底（約2.71828）
$k!$は階乗
$\lambda\$は平均発生率

で、ポアソン分布は非負整数値をとります。

$$
X = 0, 1, 2, \dots
$$

パラメータ$\lambda$の意味

$\lambda$は平均発生率です。

たとえば：

1時間あたりの平均着信件数が5件なら$\lambda = 5$
1日あたりの平均事故件数が0.3件なら$\lambda = 0.3$

このとき期待値は

$$
E[X]= \lambda
$$

です。

つまり、「ポアソン分布ではパラメータそのものが平均値」になります。

ポアソン分布の主要な統計量

$X$がポアソン分布に従うとき、主な統計量は次の通りです。

項目	式
期待値（平均）	$E[X] = \lambda$
分散	$\operatorname{Var}(X) = \lambda$
標準偏差	$\sqrt{\lambda}$
変動係数	$1/\sqrt{\lambda}$
歪度（skewness）	$1/\sqrt{\lambda}$
尖度（excess kurtosis）	$1/\lambda$

とくに重要なのは、

$$
E[X]= \operatorname{Var}(X) = \lambda
$$

という点です。

ポアソン分布では平均と分散が等しくなります。

C++におけるポアソン分布の使用方法

C++においてポアソン分布を使用するには、<random>にて定義されているstd::poisson_distribution クラスを使います。

#include <random>

template<class IntType = int>
class poisson_distribution;

コンストラクタに平均値$\lambda$（mean）を渡し、乱数エンジンと組み合わせて呼び出すことで、その$\lambda$に従うポアソン分布からサンプリングされた整数値が得られます。

#include <iostream>
#include <random>

int main()
{
    // 乱数エンジンの準備
    std::random_device rd;
    std::mt19937 gen(rd());

    // 平均 λ = 5 のポアソン分布
    std::poisson_distribution<int> dist(5.0);

    // 10回サンプリングしてみる
    for (int i = 0; i < 10; ++i) {
        int k = dist(gen);
        std::cout << k << " ";
    }

    std::cout << std::endl;

    return 0;
}

実行するたびに3 4 6 5 5 7 4 5 3 6のように、平均5前後の整数がランダムに出力されます。これはたとえば「1時間あたり平均5件の電話がかかってくる窓口で、実際に観測される着信件数」をそのままシミュレーションしていることになります。

コンストラクタとパラメータ

std::poisson_distribution<IntType> dist(mean);
- meanはdouble型で、内部的にはparam_typeとして保持されます。
- IntTypeは生成される値の型で、デフォルトはintです。
dist.mean()で現在設定されている$\lambda$を取得できます。
dist.reset()は、分布が内部に持つキャッシュされた状態（実装依存）をリセットします。エンジンを取り換えたときなどに呼んでおくと安全です。

パラメータを後から変更する

スペクトルのように「$\lambda$がビンごとに変わる」データを作るとき、ビンの数だけ poisson_distributionオブジェクトを新規に作るのはムダが多くなります。そこで便利なのがparam()メソッドです。

std::poisson_distribution<int> dist(1.0); // 初期値は何でもよい

std::poisson_distribution<int>::param_type p1(5.0);
dist.param(p1);
int k1 = dist(gen); // λ=5 でサンプリング

std::poisson_distribution<int>::param_type p2(20.0);
dist.param(p2);
int k2 = dist(gen); // λ=20 でサンプリング

同じ分布オブジェクトを使い回しながら平均値だけを切り替えられるので、ループの中で$\lambda$を変化させる用途に向いています。

スペクトルデータのシミュレーション

物理計測やセンサーデータの分析では、次のようなデータがよく登場します。

横軸：エネルギーや波長、時間など（ビン番号$x$）
縦軸：そのビンで観測されたカウント数（イベント数）

理想的な（ノイズのない）スペクトルは、ある関数$\mu(x)$で表されるとします。たとえば「ピークが、なだらかなバックグラウンドの上に乗っている」というのは分光や放射線計測の典型例です。

$$
\mu(x) = \operatorname{Background} + \sum\operatorname{Peak}
$$

しかし実際の検出器では、各ビンで観測される「本当のカウント数」は決して$\mu(x)$にはなりません。放射性崩壊や光子検出のようなランダムな計数過程は、統計的にポアソン分布に従うことが知られています。つまり、

ビン$x$での期待カウント数（真の強度）が$μ(x)$であるとき
実際に観測されるカウント数$k$は$\operatorname{Poisson}(\mu(x))$からのサンプル

というモデルになります。したがって、「関数を平均値とするスペクトルデータ」を作るとは、各ビンごとに$\lambda = \mu(x)$を計算し、そのビンだけstd::poisson_distributionからサンプリングするという作業に他なりません。

実装例

具体例として、1000個のビンからなるスペクトルを作ってみます。各ビンの真の強度$\mu(x)$は「一定のバックグラウンド + ガウシアンピーク」とし、実測値はそこにポアソン揺らぎを乗せたものとします。今回は、各ビンごとに10回サンプリングし、平均を取ってみます。

#include <cmath>
#include <fstream>
#include <iostream>
#include <random>
#include <vector>

// スペクトルの「真の」平均値関数：バックグラウンド + ガウシアンピーク
double mean_function(
    double x, double x0 = 500.0, double sigma = 50.0,
    double amplitude = 200.0, double background = 20.0
)
{
    double t = x - x0;
    double gauss = amplitude * std::exp(-t * t / (2.0 * sigma * sigma));
    return background + gauss;
}

int main()
{
    const int bins = 1000;
    const int times = 10;

    std::random_device rd;
    std::mt19937 gen(rd());

    // ビンごとに mean を切り替えるので、分布オブジェクトは1つだけ用意する
    std::poisson_distribution<int> poisson;

    std::vector<double> true_mean(bins);
    std::vector<double> observed(bins);

    for (int x = 0; x < bins; x++) {
        double lambda = mean_function((double)x);
        true_mean[x] = lambda;

        // param() で平均値だけを差し替えてサンプリング
        poisson.param(std::poisson_distribution<int>::param_type(lambda));

        for(int i = 0; i < times; i++) {
            observed[x] += poisson(gen);
        }

        observed[x] /= times;
    }

    // CSVとして出力（x, 真の平均値, 観測値）
    std::ofstream ofs("spectrum.csv");

    ofs << "x, true_mean, observed\n";

    for (int x = 0; x < bins; ++x) {
        ofs << x << "," << true_mean[x] << "," << observed[x] << "\n";
    }

    return 0;
}

出力された spectrum.csvを可視化すると、true_meanは滑らかなガウシアン曲線になる一方、observedはその周りにポアソン的なギザギザのノイズが乗った、いかにも「実測データっぽい」スペクトルになります。これはまさに検出器のカウント統計をシミュレートしたい場合（放射線計測、光子計数、質量分析など）の標準的な手法です。

実装上の注意点

$\lambda$が0またはそれに近い場合

mean_function()が0を返すビンがあると、std::poisson_distributionは常に0を返します（$\lambda = 0$のポアソン分布は退化分布になるため、規格上も問題なく動作します）。ただし、$\lambda$が負になるようなパラメータ設定は未定義動作になるため、mean_function()の実装では常に非負値を返すようにしてください。

まとめ

std::poisson_distribution<IntType>は<random>ヘッダで提供される、平均$\lambda$を指定するとポアソン分布に従う整数値を返す分布クラスです。
$\lambda$をビンごとに変えたい場合は、分布オブジェクトを1つだけ用意し、param() で平均値を差し替えるのが効率的です。
「ある関数を平均値とするスペクトルデータ」は、その関数値をビンごとの$\lambda$としてポアソンサンプリングすることで作成できます。これは検出器のカウント統計をシミュレートする際の基本パターンです。

C++,スペクトル,数学

Posted by izadori

Home

【C++】いまさらだけどvectorのpush_back()とemplace_back()を比較する

コメント一覧

まだ、コメントがありません