高斯曲线拟合原理以及Python源码

高斯函数曲线拟合数学基础

为了更好的对实验数据更好的拟合使用高斯函数曲线进行拟合。
使用高斯函数拟合比多项式拟合更加合适，多项式拟合必须把曲线分为两段，高斯函数拟合是对所有数据进行整体拟合，更能够反映出数据的总体变化情况，而多项式拟合只能对数据进行分段拟合，对数据的变化趋势进行割裂。
一下给出高斯函数拟合的数值基础：


证明X^T X矩阵非奇异是很有意义的需要证明在你的数据集上最小二乘法优化方式是可行的，但是我的证明也是代入真实的坐标值点进行计算，没有办法给出lambda1在实数集上都成立不等于0，这个证明留给后来人了。能力有限没有给出证明。





为使峰值点与Geant 4模拟出的最大值点相计算，需要求解出高斯函数曲线的峰值横坐标得到峰值（公式5.21目的）。由于公式的复杂性，无法求解出解析解进而求解数值解。高斯函数曲线先增加后减少的特性，对公式5.13进行求导后得出的函数表达式在数据横坐标区间上一定有一个零点，通过二分法迭代求解法控制绝对值误差小于0.0001的条件求解峰值点的横坐标，如公式5.21。

Python求解高斯函数代码

高斯函数：

//因为实验数据明显不是单高斯函数，其实更像朗道分布，但是这个函数是拟合不出来的
//所以使用二次高斯拟合，更加符合分布
def gauss(x, *param):
    return param[0] * np.exp(-np.power(x - param[2], 2.) / (2 * np.power(param[4], 2.))) + \
           param[1] * np.exp(-np.power(x - param[3], 2.) / (2 * np.power(param[5], 2.)))

拟合高斯曲线函数

def read_data(filename):
    ll = 0
    with open(filename, 'r') as f:
        lines = f.readlines()
        line_1 = len(lines)
        try:
            data = np.zeros((line_1, 2))
            for line1 in lines:
                value = [float(s) for s in line1.split()]
                data[ll, 0] = value[0]
                data[ll, 1] = value[1]
                ll = ll + 1
        except IndexError:
            data = np.zeros((line_1, 1))
            for line1 in lines:
                value = [float(s) for s in line1.split()]
                data[ll, 0] = value[0]
                ll = ll + 1
    return data


def fit_gauss_nor():
	filename3 = 'Tdc_gau.txt'
	//读取数据，横纵坐标分为两列
    data1 = read_data(filename3)
    plt.scatter(data1[:, 0], data1[:, 1], label='Raw Data', linewidth=2, color='blue')
    font1 = {'family': 'Times New Roman',
             'weight': 'normal',
             'size': 13,
             }
    x_axis = range(580, 740, 20)
    plt.xticks(x_axis, ('2.9', '3', '3.1', '3.2', '3.3', '3.4', '3.5', '3.6'))
    plt.xlabel("Times(ns)", font1)
    plt.ylabel("Counts", font1)
    plt.legend(loc='upper left', prop=font1)
    plt.show()
    x_sort_exp = my_sort(data1[:, 0])
    mid1 = data1[0, 1]
    for i in range(0, len(data1)):   #normalization
        if data1[i, 1] > mid1:
            mid1 = data1[i, 1]
    data2 = data1
    for i in range(0, len(data2)):
        data2[i, 1] = data2[i, 1]
    //使用优化器，最小二次优化方向，和理论推导非常相似，P0为参数预估值
    popt, pcov = optimize.curve_fit(gauss, data2[:, 0], data2[:, 1], p0=[110, 125, 620, 680, 4, 11])
    print('popt = ', popt)
    mse = cal_MSE(gauss(x_sort_exp, *popt), data1[:, 1])
    print('cal_mse = ', mse)
    plt.plot(x_sort_exp, gauss(x_sort_exp, *popt), label='Fitting Function',
             linewidth=2, color='red')
    plt.scatter(data1[:, 0], data1[:, 1], label='Raw Data', linewidth=2, color='blue')
    font1 = {'family': 'Times New Roman',
             'weight': 'normal',
             'size': 13,
             }
    x_axis = range(580, 740, 20)
    plt.xticks(x_axis, ('2.9', '3', '3.1', '3.2', '3.3', '3.4', '3.5', '3.6'))
    plt.xlabel("Times(ns)", font1)
    plt.legend(loc='upper left', prop=font1)
    plt.show()
    return popt

//这个函数为传输峰值数据，进行数值填充算法，进行卷积计算
def padding_data(data):
    params = fit_gauss_nor()
    //符号求导
    x = sp.Symbol("x")
    fun1 = sp.diff(gauss_1(x, *params), x)
    print('fun1 = ', fun1)
    x_min = 620
    x_max = 660
    x_mid = (x_max + x_min) / 2
    //迭代计算峰值横坐标
    while abs(fun1.subs(x, x_mid)) > 0.0001:
        if fun1.subs(x, x_mid) < 0:
            x_max = x_mid
        elif fun1.subs(x, x_mid) > 0:
            x_min = x_mid
        x_mid = (x_max + x_min) / 2
    x_peak = x_mid
    x_plot_min = 590
    x_plot_max = 700
    point_number = 50
    x_point1 = frange(x_plot_min, x_peak, (x_peak - x_plot_min)/point_number)
    x_point2 = frange(x_peak, x_plot_max, (x_plot_max - x_peak)/point_number)
    y_point1 = gauss(x_point1, *params)
    y_point2 = gauss(x_point2, *params)
    x_point1[len(x_point1):len(x_point1)] = x_point2
    y_point1 = np.append(y_point1, y_point2)
    plt.plot(x_point1, y_point1, label='Fitting Function',
             linewidth=2, color='red')
    plt.show()
    data_array = []
    for i in range(0, len(y_point1)):
        data_array.append(y_point1[i] * data)
    return data_array

结论

以上整个博客给出了高斯函数曲线拟合的数值理论基础与Python实现代码，希望对大家有帮助，欢迎大家点赞收藏！