【C#】MathNet.Numericsの最適化を使って最尤推定してみる

これまで、確率分布のパラメータ推定などの最尤推定をするときは、RのOptim関数などを使っていました。

たとえばこちらの記事:[R]尤度を可視化してみた&最尤推定について

だいたいの場合はRで十分なのですけど、多数の標本セットを対象とするような場合だと手続き型の言語の方が都合が良かったりします。私個人としてはC#が書き慣れているという事情もあって、C#のライブラリ「MathNet.Numerics」を使って、尤度関数の最適化によって正規分布のパラメータを推定してみました。

using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.Optimization;
using MathNet.Numerics.Distributions;

const int Length = 1000; //標本の長さ
const double Mu = 5d; //答え(μ)
const double Sd = 2d; //答え(σ)

static void Main(string[] args) {
	Random rnd = new Random(0);

	//μ=Mu、σ=Sdの正規分布に従う乱数
	double[] dX = new double[Length];
	for (int i = 0; i < Length; i++) {
		double dP = rnd.NextDouble();
		dX[i] = Normal.InvCDF(Mu, Sd, dP);
	}

	//標本と尤度関数を持つオブジェクトのインスタンス
	OptimNormDist nd = new OptimNormDist(dX);

	//最適化の準備
	var fLH = new Func<Vector<double>, double>(x => nd.LogLikelifood(x));
	var obj = ObjectiveFunction.Value(fLH);
	var init = Vector<double>.Build;

	//最適化(Nelder-Mead法…微分に頼らずに極値を求める)
	var res = NelderMeadSimplex.Minimum(obj, init.Dense(new double[] { 0d, 1d }));

	//パラメータの取得
	double dMu_est = res.FunctionInfoAtMinimum.Point[0];
	double dSd_est = res.FunctionInfoAtMinimum.Point[1];

	Console.WriteLine("μ:{0:f2}, μ_推定:{1:f4}", Mu, dMu_est);
	Console.WriteLine("σ:{0:f2}, σ_推定:{1:f4}", Sd, dSd_est);
}

/// <summary>
/// 標本と尤度関数を持つオブジェクト
/// </summary>
class OptimNormDist {
	public double[] Samples { get; }
	public OptimNormDist(double[] Samples) {
		this.Samples = Samples;
	}

	//対数尤度
	public double LogLikelifood(Vector<double> Params) {
		if (Params[1] <= 0d) { return (double.MinValue); }
		else {
			double dLogLH = 0d;
			for (int i = 0; i < Samples.Length; i++) {
				dLogLH += Math.Log(Normal.PDF(Params[0], Params[1], Samples[i]));
			}
			return (-dLogLH); //最大化したいので負の値を返す
		}
	}
}

謝辞:記事の作成にあたっては下記を参考にさせて頂きました。
Math.NET Numericsに機能が増えている