高斯曲线拟合原理以及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实现代码,希望对大家有帮助,欢迎大家点赞收藏!

  • 25
    点赞
  • 140
    收藏
    觉得还不错? 一键收藏
  • 13
    评论
评论 13
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值