【高斯拟合】不用库手写高斯拟合算法：从最小二乘到拟合参数推导

高斯分布（正态分布）是科学计算和机器学习中最常见的函数之一，拟合一组数据为高斯曲线在信号处理、统计建模、图像处理中都有广泛应用。

市面上很多工具包（如 NumPy、SciPy）都可以快速进行高斯拟合。但你有没有想过，如果不使用任何库函数，只使用原始数学知识和基础 Python 语法，我们也可以完成一次完整的高斯拟合？

本文将一步步讲解如何手写一个 高斯拟合算法，核心思路是：

• 将高斯函数取对数后，转化为一个二次函数
• 使用最小二乘法拟合该二次函数
• 反推出高斯分布的参数：振幅 a，中心 mu，标准差 sigma

一、高斯函数公式

高斯函数的表达式如下：

$$f(x)=a\cdot \exp \left(-\frac{(x-\mu {)}^2}{2{\sigma }^2}\right)$$

其中：

• $a$：峰值（振幅）
• $\mu$：均值（峰值位置）
• $\sigma$：标准差（控制曲线宽度）

二、对高斯函数取对数

我们不能直接拟合指数函数，所以我们先对其取自然对数：

$$\ln (f(x))=\ln (a)-\frac{(x-\mu {)}^2}{2{\sigma }^2}$$

整理为一个关于 $x$ 的二次函数：

$$\ln (f(x))=A{x}^2+Bx+C$$

其中：

• $A=-\frac{1}{2{\sigma }^2}$
• $B=\frac{\mu }{{\sigma }^2}$
• $C=\ln (a)-\frac{{\mu }^2}{2{\sigma }^2}$

所以我们可以将目标转为一个二次曲线拟合问题，只要拟合出 $A,B,C$，就可以反推高斯参数 $a,\mu ,\sigma$。

三、最小二乘法拟合二次函数

我们要拟合下列形式的曲线：

$$z=A{x}^2+Bx+C$$

其中 $z=\ln (y)$。

构造三个正则方程来求解 $A,B,C$：

这是一个 $3\times 3$ 的线性方程组，可以使用克拉默法则（Cramer’s Rule）或高斯消元来求解。

四、从拟合结果还原高斯参数

根据上文推导公式反解：

• $\sigma =\sqrt{-\frac{1}{2A}}$
• $\mu =-\frac{B}{2A}$
• $a=\exp \left(C-\frac{B^2}{4A}\right)$

这就完成了高斯参数的恢复。

五、完整 Python 实现（不依赖库）

import math

def gaussian_fit_log_least_squares(x_vals, y_vals):
    assert len(x_vals) == len(y_vals)
    n = len(x_vals)

    S_x4 = S_x3 = S_x2 = S_x1 = S_1 = 0.0
    S_zx2 = S_zx1 = S_z = 0.0

    valid_points = 0

    for i in range(n):
        x = x_vals[i]
        y = y_vals[i]
        if y <= 0:
            continue

        z = math.log(y)
        x2 = x * x
        x3 = x2 * x
        x4 = x2 * x2

        S_x4 += x4
        S_x3 += x3
        S_x2 += x2
        S_x1 += x
        S_1 += 1

        S_zx2 += x2 * z
        S_zx1 += x * z
        S_z += z

        valid_points += 1

    if valid_points < 3:
        raise ValueError("有效点太少，无法拟合")

    A = [
        [S_x4, S_x3, S_x2],
        [S_x3, S_x2, S_x1],
        [S_x2, S_x1, S_1]
    ]
    b = [S_zx2, S_zx1, S_z]

    def solve_linear_3x3(A, b):
        def det3(m):
            return (
                m[0][0] * (m[1][1]*m[2][2] - m[1][2]*m[2][1]) -
                m[0][1] * (m[1][0]*m[2][2] - m[1][2]*m[2][0]) +
                m[0][2] * (m[1][0]*m[2][1] - m[1][1]*m[2][0])
            )

        D = det3(A)
        if abs(D) < 1e-12:
            raise ValueError("矩阵不可逆")

        def replace_col(m, col_idx, new_col):
            new_m = [row[:] for row in m]
            for i in range(3):
                new_m[i][col_idx] = new_col[i]
            return new_m

        Dx = det3(replace_col(A, 0, b))
        Dy = det3(replace_col(A, 1, b))
        Dz = det3(replace_col(A, 2, b))

        return Dx / D, Dy / D, Dz / D

    A_, B_, C_ = solve_linear_3x3(A, b)

    if A_ >= 0:
        raise ValueError("拟合失败：A 应该为负")

    sigma = math.sqrt(-1 / (2 * A_))
    mu = -B_ / (2 * A_)
    a = math.exp(C_ - (B_ * B_) / (4 * A_))

    return a, mu, sigma

六、示例数据生成 + 拟合验证

import random

def generate_data(a, mu, sigma, num_points=100):
    x_vals = []
    y_vals = []
    for i in range(num_points):
        x = -1 + i * 0.05
        y = a * math.exp(-((x - mu) ** 2) / (2 * sigma ** 2)) + random.uniform(-0.05, 0.05)
        y = max(y, 0.01)
        x_vals.append(x)
        y_vals.append(y)
    return x_vals, y_vals

x_vals, y_vals = generate_data(a=5.0, mu=2.0, sigma=1.5)

a_fit, mu_fit, sigma_fit = gaussian_fit_log_least_squares(x_vals, y_vals)
print(f"拟合结果: a={a_fit:.4f}, mu={mu_fit:.4f}, sigma={sigma_fit:.4f}")

输出示例：

拟合结果: a=5.0123, mu=2.0009, sigma=1.4988

七、小结

本文手把手地实现了一个不用任何库的高斯拟合算法，从数学推导、线性方程组求解到最终还原高斯参数。该方法适用于单峰数据的拟合，能帮我们更深入理解最小二乘法与高斯分布之间的数学联系。

如果你希望进一步了解多峰高斯拟合（比如双高斯拟合）、使用更鲁棒的拟合算法（如最大似然、非线性最小二乘法），可以尝试引入 NumPy / SciPy 进行下一步扩展。

📌 附：该方法的优缺点

如果你觉得本文对你有帮助，欢迎点赞、收藏或转发给有需要的小伙伴！