在数据分析与清洗中经常遇到这样的数据：数据不仅仅向单个中心靠拢，而是类似分段的向两个甚至多个中心靠拢。数据向单个中心靠拢时，可以简单地用高斯函数去拟合，得到相关参数，那么可以通过3sigma准则对数据进行去噪以为后续处理做准备。当数据向两个中心靠拢时，单峰高斯函数就不再适用了，这时候就需要双峰高斯函数出马了。

高斯函数（Gaussian Function）是一种常见的数学函数，广泛应用于概率论、统计学、信号处理、图像处理等领域。

---

一维高斯函数

一维高斯函数的公式为：
$$f(x)=A\cdot \exp \left(-\frac{(x-\mu {)}^2}{2{\sigma }^2}\right)$$

参数说明：

• $A$：振幅（Amplitude），控制函数的最大值。
• $\mu$：均值（Mean），表示函数的中心位置（即峰值所在的位置）。
• $\sigma$：标准差（Standard Deviation），控制函数的宽度（越小越尖锐，越大越平缓）。
• $x$：自变量。

归一化形式：
如果需要将高斯函数归一化（使其在 ($-\infty$) 到 ($+\infty$) 的积分等于 1），则公式变为：
$$f(x)=\frac{1}{\sqrt{2\pi }\sigma }\cdot \exp \left(-\frac{(x-\mu {)}^2}{2{\sigma }^2}\right)$$
此时，( A = \frac{1}{\sqrt{2\pi}\sigma} )。

---

二维高斯函数

二维高斯函数通常用于描述二维空间中的分布，例如图像处理中的模糊滤波器。其公式为：
$$f(x,y)=A\cdot \exp \left(-\frac{(x-x_0{)}^2}{2{\sigma }_x^2}-\frac{(y-y_0{)}^2}{2{\sigma }_y^2}\right)$$

各向同性高斯函数
当 ( \sigma_x = \sigma_y = \sigma ) 时，二维高斯函数变为：
$$f(x,y)=A\cdot \exp \left(-\frac{(x-x_0{)}^2+(y-y_0{)}^2}{2{\sigma }^2}\right)$$

---

多维高斯函数

对于 n 维高斯函数，其通用形式为：

$$f(\mathbf{x})=A\cdot \exp \left(-\frac{1}{2}(\mathbf{x}-\mathbf{\mu }{)}^{\top }{\Sigma }^{-1}(\mathbf{x}-\mathbf{\mu })\right)$$

参数说明：

• $\mathbf{x}$：( n )-维向量，表示自变量。
• $\mathbf{\mu }$：( n )-维向量，表示均值。
• $\Sigma$：协方差矩阵，控制分布的形状和方向。
• ${\Sigma }^{-1}$：协方差矩阵的逆矩阵。

各向同性多维高斯函数
当协方差矩阵为对角矩阵且所有对角元素相等时（即$\Sigma ={\sigma }^{2}I$，其中 $I$ 是单位矩阵），公式简化为：
$$f(\mathbf{x})=A\cdot \exp \left(-\frac{\Vert \mathbf{x}-\mathbf{\mu }{\Vert }^2}{2{\sigma }^2}\right)$$
其中 $\Vert \mathbf{x}-\mathbf{\mu }\Vert$ 表示欧几里得距离。

---

一维双峰高斯函数

对于双峰高斯函数，我们可以将其视为两个这样的高斯函数的叠加，每个都有自己的振幅，均值及标准差。因此，双峰高斯函数的形式可以表示为：

$$f(x)=A_1\exp \left(-\frac{(x-{\mu }_1{)}^2}{2{\sigma }_1^2}\right)+A_2\exp \left(-\frac{(x-{\mu }_2{)}^2}{2{\sigma }_2^2}\right)$$

可以用于拟合或描述具有双峰特性的数据集。

---

代码实现

这里基于ceres，在C++中对双峰高斯拟合进行实现。
思路很简单，输入一维的数据和直方图的宽度，构建直方图，然后通过直方图的顶端坐标来拟合双峰高斯函数，稍微麻烦的是找到合适的初始值：

1. 建直方图
2. 确定初值
3. 优化求解

struct BimodalGaussianResidual
{BimodalGaussianResidual(double x, double y) : x_(x), y_(y){}template <typename T>bool operator()(const T* const params, T* residual) const{const T A1 = params[0];const T mu1 = params[1];const T sigma1 = params[2];const T A2 = params[3];const T mu2 = params[4];const T sigma2 = params[5];T f1 = A1 * ceres::exp(-(x_ - mu1) * (x_ - mu1) / (T(2) * sigma1 * sigma1));T f2 = A2 * ceres::exp(-(x_ - mu2) * (x_ - mu2) / (T(2) * sigma2 * sigma2));residual[0] = y_ - (f1 + f2);return true;}private:const double x_;const double y_;
};bool FitBimodalGaussian(std::vector<float>& data, double* params, const double  bin_width)
{//==========================建直方图==========================double min_val = *std::min_element(data.begin(), data.end());double max_val = *std::max_element(data.begin(), data.end());std::map<double, int> histogram;for (double value : data){double bin = std::floor((value - min_val) / bin_width) * bin_width; // 找到所属的binhistogram[bin]++;}std::vector<double> x_data;std::vector<double> heights;for (const auto& entry : histogram){x_data.push_back(entry.first + min_val + 0.5 * bin_width);heights.push_back(entry.second);}//==========================确定初值==========================// 找到直方图中的两个峰值double max_height1 = 0, max_height2 = 0;double mu1 = 0, mu2 = 0;size_t peak1_index = 0; // 第一个峰值的索引// 找第一个峰值for (size_t i = 0; i < heights.size(); ++i){if (heights[i] > max_height1){max_height1 = heights[i];mu1 = x_data[i];peak1_index = i;}}// 找第二个峰值（避开第一个峰值附近的区域）const double separation_threshold = 5; // 峰值之间的最小距离      // 根据数据确定for (size_t i = 0; i < heights.size(); ++i){if (std::abs(x_data[i] - mu1) > separation_threshold && heights[i] > max_height2){max_height2 = heights[i];mu2 = x_data[i];}}// 估计标准差 sigma1 和 sigma2double sigma1 = 0.5, sigma2 = 0.5; // 初始猜测值const double half_max_height1 = max_height1 / 2.0;const double half_max_height2 = max_height2 / 2.0;// 找第一个峰值的半高宽double left_edge1 = mu1, right_edge1 = mu1;for (size_t i = 0; i < x_data.size(); ++i){if (x_data[i] < mu1 && heights[i] >= half_max_height1){left_edge1 = x_data[i];}if (x_data[i] > mu1 && heights[i] >= half_max_height1){right_edge1 = x_data[i];break;}}sigma1 = (right_edge1 - left_edge1) / (2 * std::sqrt(2 * std::log(2)));// 找第二个峰值的半高宽double left_edge2 = mu2, right_edge2 = mu2;for (size_t i = 0; i < x_data.size(); ++i){if (x_data[i] < mu2 && heights[i] >= half_max_height2){left_edge2 = x_data[i];}if (x_data[i] > mu2 && heights[i] >= half_max_height2){right_edge2 = x_data[i];break;}}sigma2 = (right_edge2 - left_edge2) / (2 * std::sqrt(2 * std::log(2)));// 设置初始参数 [A1, mu1, sigma1, A2, mu2, sigma2]params[0] = max_height1;params[1] = mu1;params[2] = sigma1;params[3] = max_height2;params[4] = mu2;params[5] = sigma2;#if MY_LOGstd::cout << std::fixed << std::setprecision(6);std::cout << "\nInitial Parameters:" << std::endl;std::cout << "A1 = " << params[0] << ", mu1 = " << params[1] << ", sigma1 = " << params[2] << std::endl;std::cout << "A2 = " << params[3] << ", mu2 = " << params[4] << ", sigma2 = " << params[5] << std::endl;
#endif//==========================优化求解==========================ceres::Problem problem;for (size_t i = 0; i < x_data.size(); ++i){problem.AddResidualBlock(new ceres::AutoDiffCostFunction<BimodalGaussianResidual, 1, 6>(new BimodalGaussianResidual(x_data[i], heights[i])), nullptr, params);}ceres::Solver::Options options;options.linear_solver_type = ceres::DENSE_QR;options.minimizer_progress_to_stdout = false;ceres::Solver::Summary summary;ceres::Solve(options, &problem, &summary);#if MY_LOGstd::cout << "\nFinal Parameters:\n";std::cout << "A1: " << params[0] << ", mu1: " << params[1] << ", sigma1: " << params[2] << "\n";std::cout << "A2: " << params[3] << ", mu2: " << params[4] << ", sigma2: " << params[5] << "\n";
#endifreturn true;
}

用一组真实数据测试，可得：

Initial Parameters:
A1 = 16223.000000, mu1 = -1.018652, sigma1 = 0.084932
A2 = 2868.000000, mu2 = 7.181348, sigma2 = 0.084932Final Parameters:
A1: 15804.267874, mu1: -1.138264, sigma1: 0.744859
A2: 2530.132897, mu2: 6.820187, sigma2: 0.759935

二维双峰高斯函数

类似的，二维双峰高斯函数可表达式如下：
$$f(x,y)=A_1\exp \left(-\frac{(x-{\mu }_{1x}{)}^2}{2{\sigma }_{1x}^2}-\frac{(y-{\mu }_{1y}{)}^2}{2{\sigma }_{1y}^2}\right)+A_2\exp \left(-\frac{(x-{\mu }_{2x}{)}^2}{2{\sigma }_{2x}^2}-\frac{(y-{\mu }_{2y}{)}^2}{2{\sigma }_{2y}^2}\right)$$

代码实现

#include <iostream>
#include <vector>
#include <cmath>
#include <random>
#include <map>
#include <iomanip>  
#include "ceres/ceres.h"
#include "gnuplot-iostream.h"// 二维双峰高斯残差函数
struct BimodalGaussianResidual2D
{BimodalGaussianResidual2D(double x, double y, double z) : x_(x), y_(y), z_(z){}template <typename T>bool operator()(const T* const params, T* residual) const{const T A1 = params[0];const T mu1_x = params[1];const T mu1_y = params[2];const T sigma1_x = params[3];const T sigma1_y = params[4];const T A2 = params[5];const T mu2_x = params[6];const T mu2_y = params[7];const T sigma2_x = params[8];const T sigma2_y = params[9];T f1 = A1 * ceres::exp(-(T(x_ - mu1_x) * T(x_ - mu1_x) / (T(2) * sigma1_x * sigma1_x) +T(y_ - mu1_y) * T(y_ - mu1_y) / (T(2) * sigma1_y * sigma1_y)));T f2 = A2 * ceres::exp(-(T(x_ - mu2_x) * T(x_ - mu2_x) / (T(2) * sigma2_x * sigma2_x) +T(y_ - mu2_y) * T(y_ - mu2_y) / (T(2) * sigma2_y * sigma2_y)));residual[0] = z_ - (f1 + f2);return true;}private:const double x_;const double y_;const double z_;
};// 二维双峰高斯数据生成
void GenerateBimodalGaussianData(std::vector<double>& x_data, std::vector<double>& y_data, std::vector<double>& z_data,double A1, double mu1_x, double mu1_y, double sigma1_x, double sigma1_y,double A2, double mu2_x, double mu2_y, double sigma2_x, double sigma2_y,int num_points, double noise_level)
{std::default_random_engine generator;std::normal_distribution<double> distribution(0.0, noise_level);for (int i = 0; i < num_points; ++i){double x = static_cast<double>(rand()) / RAND_MAX * 10.0 - 5.0; // 范围 [-5, 5]double y = static_cast<double>(rand()) / RAND_MAX * 10.0 - 5.0; // 范围 [-5, 5]double f1 = A1 * exp(-(pow(x - mu1_x, 2) / (2 * pow(sigma1_x, 2)) +pow(y - mu1_y, 2) / (2 * pow(sigma1_y, 2))));double f2 = A2 * exp(-(pow(x - mu2_x, 2) / (2 * pow(sigma2_x, 2)) +pow(y - mu2_y, 2) / (2 * pow(sigma2_y, 2))));double z = f1 + f2 + distribution(generator); // 添加噪声x_data.push_back(x);y_data.push_back(y);z_data.push_back(z);}
}// 二维双峰高斯函数拟合
bool FitBimodalGaussian2D(const std::vector<double>& x_data, const std::vector<double>& y_data, const std::vector<double>& z_data, double* params)
{ceres::Problem problem;for (size_t i = 0; i < x_data.size(); ++i){problem.AddResidualBlock(new ceres::AutoDiffCostFunction<BimodalGaussianResidual2D, 1, 10>(new BimodalGaussianResidual2D(x_data[i], y_data[i], z_data[i])),nullptr, params);}ceres::Solver::Options options;options.linear_solver_type = ceres::DENSE_QR;options.minimizer_progress_to_stdout = false;ceres::Solver::Summary summary;ceres::Solve(options, &problem, &summary);return true;
}int main()
{//==========================数据生成==========================std::vector<double> x_data, y_data, z_data;int num_points = 1000;double noise_level = 0.2;double A1_true = -3.0, mu1_x_true = -2.0, mu1_y_true = 2.0, sigma1_x_true = 1.0, sigma1_y_true = 1.0;double A2_true = 5.0, mu2_x_true = 2.0, mu2_y_true = -2.0, sigma2_x_true = 1.5, sigma2_y_true = 1.5;GenerateBimodalGaussianData(x_data, y_data, z_data,A1_true, mu1_x_true, mu1_y_true, sigma1_x_true, sigma1_y_true,A2_true, mu2_x_true, mu2_y_true, sigma2_x_true, sigma2_y_true,num_points, noise_level);//==========================确定初值==========================double initial_params[] = { A1_true, mu1_x_true, mu1_y_true,  sigma1_x_true, sigma1_y_true, A2_true , mu2_x_true ,  mu2_y_true, sigma2_x_true, sigma2_y_true };// 添加随机噪声std::random_device rd;std::mt19937 gen(rd());std::normal_distribution<> d(0, 2*noise_level);for (int i = 0; i < 10;++i)initial_params[i]+= d(gen);//==========================优化求解==========================std::cout << "\nInitial Parameters:\n";std::cout << "A1: " << initial_params[0] << ", mu1_x: " << initial_params[1] << ", mu1_y: " << initial_params[2] << ", sigma1_x: " << initial_params[3] << ", sigma1_y: " << initial_params[4] << "\n";std::cout << "A2: " << initial_params[5] << ", mu2_x: " << initial_params[6] << ", mu2_y: " << initial_params[7] << ", sigma2_x: " << initial_params[8] << ", sigma2_y: " << initial_params[9] << "\n";FitBimodalGaussian2D(x_data, y_data, z_data, initial_params);std::cout << "\nFinal Fitted Parameters:\n";std::cout << "A1: " << initial_params[0] << ", mu1_x: " << initial_params[1] << ", mu1_y: " << initial_params[2] << ", sigma1_x: " << initial_params[3] << ", sigma1_y: " << initial_params[4] << "\n";std::cout << "A2: " << initial_params[5] << ", mu2_x: " << initial_params[6] << ", mu2_y: " << initial_params[7] << ", sigma2_x: " << initial_params[8] << ", sigma2_y: " << initial_params[9] << "\n";//==========================图形绘制==========================// 将数据传递给 Gnuplot 并绘制std::vector<std::tuple<double, double, double>> data;for (size_t i = 0; i < x_data.size(); ++i){data.emplace_back(x_data[i], y_data[i], z_data[i]);}auto format_param = [](double value){std::ostringstream oss;oss << std::fixed << std::setprecision(3) << value;return oss.str();};Gnuplot gp;// 设置绘图标题和坐标轴标签gp << "set terminal wxt enhanced\n";    // 使用支持交互的终端gp << "set mouse\n";                    // 启用鼠标交互gp << "set title 'BimodalGaussian2D Fit'\n";gp << "set xlabel 'X'\n";gp << "set ylabel 'Y'\n";gp << "set zlabel 'Z'\n";// 设置绘图为三维模式gp << "set style data points\n";    // 使用点样式gp << "set ticslevel 0\n";          // 调整 Z 轴刻度位置gp << "set isosamples 100,100\n";   // 设置采样数量gp << "set hidden3d\n";             // 隐藏被遮挡区域gp << "splot '-' with points pointtype 7 pointsize 0.5 title 'Scatter Points', "<< format_param(initial_params[0]) << "*exp(-((x - " << format_param(initial_params[1])<< ")**2/(2*" << format_param(initial_params[3]) << "**2) + (y - " << format_param(initial_params[2])<< ")**2/(2*" << format_param(initial_params[4]) << "**2))) "<< "+ " << format_param(initial_params[5]) << "*exp(-((x - " << format_param(initial_params[6])<< ")**2/(2*" << format_param(initial_params[8]) << "**2) + (y - " << format_param(initial_params[7])<< ")**2/(2*" << format_param(initial_params[9]) << "**2)))\n";gp.send(data);gp.flush();return 0;
}

在真值基础上添加2倍随机噪声的初值及拟合结果如下所示，可以看到其拟合得相当准确的。

Initial Parameters:
A1: -3.11395, mu1_x: -2.43491, mu1_y: 1.9397, sigma1_x: 1.21568, sigma1_y: 1.3008
A2: 4.93584, mu2_x: 2.1473, mu2_y: -1.74558, sigma2_x: 1.88884, sigma2_y: 1.73732Final Fitted Parameters:
A1: -3.00556, mu1_x: -2.01663, mu1_y: 1.98635, sigma1_x: 0.96315, sigma1_y: 1.01506
A2: 4.95913, mu2_x: 2.01769, mu2_y: -2.01035, sigma2_x: 1.49838, sigma2_y: 1.49798

---

多维多峰高斯函数

同样，理论上也可以发散到多维多峰高斯函数的情形，不过一方面这种情况比较少见，另一方面高维的初值一般不好自动确定，可能应用得很少，这里不再作说明和演示。

---

补充
对ceres拟合感兴趣可移步C++ 带约束的Ceres形状拟合。
对gnuplot绘制感兴趣可移步如何在C++中优雅地绘制图表。

---

打完收工。