使用最大似然策略学习正态分布参数

最近在看计算机视觉：模型学习与推理这本书，书中讲了几十种算法，原书是使用matlab实现的，本人打算用c++重新实现一遍。书中作者也推荐自己实现一遍代码对于理解问题大有裨益。

本篇博客主要讲使用c++标准库生成正态分布的数据，然后使用最大似然的方法估计学习参数。并提供c++源码。

一维正态分布满足如下公式：

$$\mathrm{Pr}(x)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left[-0.5(x-\mu {)}^2/{\sigma }^2\right]$$

n个数据假设独立分布，则联合分布为：
$$\mathrm{Pr}\left(x_{1\dots I}|\mu ,{\sigma }^2\right)=\prod _{i=1}^I\mathrm{Pr}\left(x_i|\mu ,{\sigma }^2\right)$$

$$\begin{array}{l}={\prod }_{i=1}^I{\mathrm{Norm}}_{x_i}\left[\mu ,{\sigma }^2\right]\\ =\frac{1}{{\left(2\pi {\sigma }^2\right)}^{I/2}}\exp \left[-0.5{\sum }_{i=1}^I\frac{{\left(x_i-\mu \right)}^2}{{\sigma }^2}\right]\end{array}$$

进行极大似然估计，就是选择一组参数，使得目前观测的联合概率分布达到最大，那么这组参数就是最符合目前观测数据的。
常规方法是求对数的导数，计算极值：

$$\begin{array}{rl}\hat{\mu },{\hat{\sigma }}^2& =\underset{\mu ,{\sigma }^2}{\mathrm{argmax}}\left[\sum _{i=1}^I\log \left[{\mathrm{Norm}}_{x_i}\left[\mu ,{\sigma }^2\right]\right]\right]\\ & =\underset{\mu ,{\sigma }^2}{\mathrm{argmax}}\left[-0.5I\log [2\pi ]-0.5I\log {\sigma }^2-0.5\sum _{i=1}^I\frac{{\left(x_i-\mu \right)}^2}{{\sigma }^2}\right]\end{array}$$

求极值：

$$\begin{array}{rl}\frac{\partial L}{\partial \mu }& =\sum _{i=1}^I\frac{\left(x_i-\mu \right)}{{\sigma }^2}\\ & =\frac{{\sum }_{i=1}^I{x}_i}{{\sigma }^2}-\frac{I\mu }{{\sigma }^2}=0\end{array}$$

$$\hat{\mu }=\frac{{\sum }_{i=1}^I{x}_i}{I}$$
$${\hat{\sigma }}^2=\sum _{i=1}^I\frac{{\left(x_i-\hat{\mu }\right)}^2}{I}$$

实际上我们的最终算法就是这两个公式。算法整体流程如下:
$$\begin{array}{l}\text{ Input : Training data }{\left\{x_i\right\}}_{i=1}^I\\ \text{ Output: Maximum likelihood estimates of parameters }\theta =\left\{\mu ,{\sigma }^2\right\}\\ \text{ begin }\\ //\text{ Set mean parameter }\\ \mu ={\sum }_{i=1}^I{x}_i/I\\ //\text{ Set variance }\\ {\sigma }^2={\sum }_{i=1}^I{\left(x_i-\hat{\mu }\right)}^2/I\\ \text{ end }\end{array}$$

数据生成的代码如下：

#include <random>
#include <iostream>
#include <vector>

using std::vector;

template <typename T>
vector<T> generate_normal_distribution_data(T mu,T sigma,int number)
{
	vector<T> data;

	std::random_device rd{};
	std::mt19937 gen{ rd() };
	std::normal_distribution<> d{ mu,sigma*sigma };

	for (int n = 0; n < number; ++n) {
		data.push_back(d(gen));
	}

	return data;
}

学习代码如下:

void learning_normal_distribution_parameters()
{
	vector<double> data=generate_normal_distribution_data<double>(0, 1, 100000);
	double mu = 0.0;
	double sigma = 0.0;
	for (int i = 0; i < data.size(); i++)
	{
		mu += data[i];
	}

	mu = mu / data.size();

	for (int i = 0; i < data.size(); i++)
	{
		sigma += (data[i] - mu)*(data[i] - mu);
	}

	sigma = sigma / data.size();
	double var = sigma;
	cout << "mu: " << mu << endl;
	cout << "var: " << var << endl;
}

我们在模型的参数设计上，生成均值为0，方差为1的10000个数据点。
学习得到的参数如下：