C++でガウス分布

はじめに

C++初心者が自分でガウス分布を実装します。別記事で実装した行列演算クラス(C++で行列計算実装)を用いて多次元版のガウス分布も同時に実装します。

なお、行列計算を自分で実装したものを用いていますが、基本的には行列計算ライブラリのEigenなどで置き換えても同様の実装になるかと思います。

使い方のイメージ

例えば平均０、分散１のガウス分布において0.3における値を生成したい場合、

double x = gauss(0.3,0,1) //xに0.381が入る

のように計算します。

実装

1次元ガウス分布

1次元の場合ガウス分布の式は以下で与えられます。

\mathcal{N}(x|\mu,\sigma^2)=\frac{1}{(2\pi\sigma^2)^{1/2}}\exp \biggl( -\frac{1}{2\sigma^2}(x-\mu)^2\biggl)

定義通りに実装して、下記になります。

double Distributions::gauss(double x, double mu, double sigma){
  return 1 / pow(2 * M_PI * pow(sigma,2), 0.5) * exp(-1/(2 * pow(sigma,2)) * pow(x - mu,2));
}

ただし、Distributionsという名前空間内で実装しました。

多次元ガウス分布

多次元の場合ガウス分布の式は下記です。

\mathcal{N}(\bf{x}|\boldsymbol{\mu},\boldsymbol{\Sigma})=\frac{1}{(2 \pi)^{D/2}}\frac{1}{|\boldsymbol{\Sigma}|^{1/2}}\exp \biggl( -\frac{1}{2}(\bf{x}-\boldsymbol{\mu})^{T}\Sigma^{-1}(\bf{x}-\boldsymbol{\mu})\biggl)

定義通りに実装すると、

double Distributions::gauss(Ma x, Ma mu, Ma sigma){
  int D = sigma.shape()[0];
  return 1 / pow(pow(2 * M_PI, D),0.5) * 1 / pow(sigma.det(), 0.5) * exp(-0.5 * ((x-mu).T() * sigma.inv() * (x-mu))(0,0));
}

ここで以前投稿した記事で実装した行列計算クラスMa(C++で行列計算実装)を入力として受け取っています。

行列計算クラスを使うと多次元の場合であっても簡単に実装できました。

まとめ

ガウス分布をC++で実装しました。ガウス分布以外のベータ分布やディリクレ分布なども簡単に実装できそうです。

コード

dist.h

#include "mat.h"
namespace Distributions{
  double gauss(Ma x, Ma mu, Ma sigma);
  double gauss(double x, double mu, double sigma);
}

dist.cpp

#include "mat.h"
#include <math.h>
#include "dist.h"


double Distributions::gauss(double x, double mu, double sigma){
  return 1 / pow(2 * M_PI * pow(sigma,2), 0.5) * exp(-1/(2 * pow(sigma,2)) * pow(x - mu,2));
}
double Distributions::gauss(Ma x, Ma mu, Ma sigma){
  int D = sigma.shape()[0];
  return 1 / pow(pow(2 * M_PI, D),0.5) * 1 / pow(sigma.det(), 0.5) * exp(-0.5 * ((x-mu).T() * sigma.inv() * (x-mu))(0,0));
}