Python直方图拟合：多高斯曲线拟合及曲线交点求解

我有一个直方图，分布存在多个众数区间，分别在5-10、110-115和180附近。希望用Python实现一条全局高斯曲线加三条局部高斯曲线拟合该直方图，同时求解三条局部高斯曲线的交点，以此划分地质断层方向的三个类别。尝试了以下代码但拟合效果不佳，恳请协助解决。

原问题代码 #

x_array = np.linspace(0, 180, 36)

amp1 = 700
sigma1 = 40
cen1 = 7

amp2 = 200
sigma2 = 10
cen2 = 65

amp3 = 750
sigma3 = 7
cen3 = 115


amp4 = 700
sigma4 = 4
cen4 = 172

y_array_4gauss = amp1*(1/(sigma1*(np.sqrt(2*np.pi))))*(np.exp((-1.0/2.0)*(((x_array-cen1)/sigma1)**2))) + \
                amp2*(1/(sigma2*(np.sqrt(2*np.pi))))*(np.exp((-1.0/2.0)*(((x_array-cen2)/sigma2)**2))) + \
                amp3*(1/(sigma3*(np.sqrt(2*np.pi))))*(np.exp((-1.0/2.0)*(((x_array-cen3)/sigma3)**2))) + \
                amp4*(1/(sigma4*(np.sqrt(2*np.pi))))*(np.exp((-1.0/2.0)*(((x_array-cen4)/sigma4)**2)))

def _4gaussianB(x, amp1,cen1,sigma1, amp2,cen2,sigma2,amp3,cen3,wid3,amp4,cen4,wid4,base):
    return base + amp1*(1/(sigma1*(np.sqrt(2*np.pi))))*(np.exp((-1.0/2.0)*(((x_array-cen1)/sigma1)**2))) + \
            amp2*(1/(sigma2*(np.sqrt(2*np.pi))))*(np.exp((-1.0/2.0)*(((x_array-cen2)/sigma2)**2)))+ \
            amp3*(1/(sigma3*(np.sqrt(2*np.pi))))*(np.exp((-1.0/2.0)*(((x_array-cen3)/sigma3)**2)))+ \
            amp4*(1/(sigma4*(np.sqrt(2*np.pi))))*(np.exp((-1.0/2.0)*(((x_array-cen4)/sigma4)**2)))
                
popt_4gauss, pcov_4gauss = scipy.optimize.curve_fit(_4gaussianB, x_array, y_array_4gauss, p0=[amp1,cen1,wid1,amp2,cen2,wid2,amp3,cen3,wid3,amp4,cen4,wid4, y_array_4gauss.min()], bounds=([0,x_array.min(),0,0,x_array.min(),0,0,x_array.min(),0,0,x_array.min(),0,0.9999*y_array_4gauss.min()],[np.inf,x_array.max(),np.inf,np.inf,x_array.max(),np.inf,np.inf,x_array.max(),np.inf,np.inf,x_array.max(),np.inf,y_array_4gauss.min()]))

# calcul des erreurs standards
perr_4gauss = np.sqrt(np.diag(pcov_4gauss))
sd_1 = perr_4gauss[0:3]
sd_2 = perr_4gauss[3:6]
sd_3 = perr_4gauss[6:9]
sd_4 = perr_4gauss[9:12]

pars_1 = popt_4gauss[0:3]
pars_2 = popt_4gauss[3:6]
pars_3 = popt_4gauss[6:9]
pars_4 = popt_4gauss[9:12]


def _1gaussian(x, amp1,cen1,sigma1):
    return amp1*(1/(sigma1*(np.sqrt(2*np.pi))))*(np.exp((-1.0/2.0)*(((x_array-cen1)/sigma1)**2)))

gauss_peak_1 = _1gaussian(x_array, *pars_1)
gauss_peak_2 = _1gaussian(x_array, *pars_2)
gauss_peak_3 = _1gaussian(x_array, *pars_3)
gauss_peak_4 = _1gaussian(x_array, *pars_4)

# calcul des différences entre série et fonction
residual_4gauss = (y_array_4gauss) - (_4gaussianB(x_array, *popt_4gauss))

# prints the fitting parameters with their errors
if popt_4gauss.size>12:
    baselinevalue=popt_4gauss[12]
    print("-------------Baseline-------------")
    print("baseline value = %0.3f" % (baselinevalue))
print("-------------Peak 1-------------")
print("amplitude = %0.2f (+/-) %0.2f" % (pars_1[0], sd_1[0]))
print("center = %0.2f (+/-) %0.2f" % (pars_1[1], sd_1[1]))
print("sigma = %0.2f (+/-) %0.2f" % (pars_1[2], sd_1[2]))
print("area = %0.2f" % np.trapz(gauss_peak_1))
print("-------------Peak 2-------------")
print("amplitude = %0.2f (+/-) %0.2f" % (pars_2[0], sd_2[0]))
print("center = %0.2f (+/-) %0.2f" % (pars_2[1], sd_2[1]))
print("sigma = %0.2f (+/-) %0.2f" % (pars_2[2], sd_2[2]))
print("area = %0.2f" % np.trapz(gauss_peak_2))
print("-------------Peak 3-------------")
print("amplitude = %0.2f (+/-) %0.2f" % (pars_3[0], sd_3[0]))
print("center = %0.2f (+/-) %0.2f" % (pars_3[1], sd_3[1]))
print("sigma = %0.2f (+/-) %0.2f" % (pars_3[2], sd_3[2]))
print("area = %0.2f" % np.trapz(gauss_peak_3))
print("-------------Peak 4-------------")
print("amplitude = %0.2f (+/-) %0.2f" % (pars_4[0], sd_4[0]))
print("center = %0.2f (+/-) %0.2f" % (pars_4[1], sd_4[1]))
print("sigma = %0.2f (+/-) %0.2f" % (pars_4[2], sd_4[2]))
print("area = %0.2f" % np.trapz(gauss_peak_4))
print("--------------------------------")

# Création de l'histogramme
hist, bins = np.histogram(valeurs, bins=36, range=(0, 180))



# Définir les graduations sur l'axe des abscisses (x) tous les 5 degrés
ax1.set_xticks(np.arange(0, 181, 5))

fig = plt.figure(figsize=(20,12))
gs = gridspec.GridSpec(2,1, height_ratios=[1,0.25])
ax1 = fig.add_subplot(gs[0])
ax2 = fig.add_subplot(gs[1])
gs.update(hspace=0) 



# Affichage de l'histogramme avec des lignes de séparation entre les bins
ax1.bar(bins[:-1], hist, width=(bins[1]-bins[0]), color='blue', alpha=0.5, edgecolor='black')
ax1.grid(axis='x', alpha=0.5, linestyle='--')

# peak 1
ax1.plot(x_array, gauss_peak_1+baselinevalue, "g")
ax1.fill_between(x_array, baseline, gauss_peak_1+baselinevalue, facecolor="green", alpha=0.5)
  
# peak 2
ax1.plot(x_array, gauss_peak_2+baselinevalue, "y")
ax1.fill_between(x_array, baseline, gauss_peak_2+baselinevalue, facecolor="yellow", alpha=0.5)  

# peak 3
ax1.plot(x_array, gauss_peak_3+baselinevalue, "y")
ax1.fill_between(x_array, baseline, gauss_peak_3+baselinevalue, facecolor="yellow", alpha=0.5)  

# peak 4
ax1.plot(x_array, gauss_peak_4+baselinevalue, "r")
ax1.fill_between(x_array, baseline, gauss_peak_4+baselinevalue, facecolor="red", alpha=0.5)  

# baseline
if baselinevalue!=0:
    baseline=np.ones(np.shape(x_array))*baselinevalue
    ax1.plot(x_array, baseline, "k-.",label='baseline')
    ax1.fill_between(x_array, 0, baseline, facecolor="grey", alpha=0.5)

# distribution
ax1.plot(x_array, y_array_4gauss, "ro")
ax1.plot(x_array, _4gaussianB(x_array, *popt_4gauss), 'k--',label='best fitting curve')

# residual
ax2.plot(x_array, residual_4gauss, "bo", label="Abs. Residual average = %0.3f" % (abs(residual_4gauss)).mean())
    


ax2.set_xlim(-5,105)
ax2.set_ylim(residual_4gauss.min()-(0.6*np.ptp(residual_4gauss)),residual_4gauss.max()+(0.2*np.ptp(residual_4gauss)))

ax2.set_xlabel("x_array",family="serif",  fontsize=12)
ax1.set_ylabel("y_array",family="serif",  fontsize=12)
ax2.set_ylabel("Res.",family="serif",  fontsize=12)

ax1.legend(loc="best")
ax2.legend(loc="best")

ax1.xaxis.set_major_locator(ticker.MultipleLocator(20))
#ax1.yaxis.set_major_locator(ticker.MultipleLocator(50))

ax2.xaxis.set_minor_locator(AutoMinorLocator(2))
ax1.yaxis.set_minor_locator(AutoMinorLocator(2))

ax1.xaxis.set_major_formatter(plt.NullFormatter())

ax1.tick_params(axis='x',which='major', direction="out", top="on", right="on", bottom="off", length=8, labelsize=8)
ax1.tick_params(axis='x',which='minor', direction="out", top="on", right="on", bottom="off", length=5, labelsize=8)
ax1.tick_params(axis='y',which='major', direction="out", top="on", right="on", bottom="off", length=8, labelsize=8)
ax1.tick_params(axis='y',which='minor', direction="out", top="on", right="on", bottom="on", length=5, labelsize=8)

ax2.tick_params(axis='x',which='major', direction="out", top="off", right="on", bottom="on", length=8, labelsize=8)
ax2.tick_params(axis='x',which='minor', direction="out", top="off", right="on", bottom="on", length=5, labelsize=8)
ax2.tick_params(axis='y',which='major', direction="out", top="off", right="on", bottom="on", length=8, labelsize=8)
ax2.tick_params(axis='y',which='minor', direction="out", top="off", right="on", bottom="on", length=5, labelsize=8)

fig.tight_layout()
#fig.savefig("fit2Gaussian_peaks_resid.png", format="png",dpi=1000)

代码核心问题分析 #

1. 参数引用混乱：_4gaussianB函数定义中使用wid3、wid4，但计算时实际调用sigma3、sigma4，且初始参数p0中使用未定义的wid1、wid2，导致拟合参数错位。
2. 硬编码全局数组：拟合函数未使用传入的x参数，而是直接调用全局x_array，导致curve_fit无法正确计算梯度。
3. 基线边界错误：基线的上下界被限制为几乎相等的值，相当于固定基线，失去拟合意义。
4. 数据目标混淆：代码同时使用生成的模拟数据y_array_4gauss和实际直方图数据hist，但拟合目标是模拟数据，而非真实直方图。
5. 无交点求解逻辑：原代码未实现局部高斯曲线交点的计算，无法完成断层类别划分。

修正后的完整代码 #

import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import matplotlib.ticker as ticker

# 模拟原始地质数据（替换为你的真实数据valeurs）
np.random.seed(42)
valeurs = np.concatenate([
    np.random.normal(7, 40, 700),
    np.random.normal(115, 7, 750),
    np.random.normal(172, 4, 700)
])
valeurs = valeurs[(valeurs >= 0) & (valeurs <= 180)]

# 生成直方图数据，用区间中点作为拟合x值
hist, bins = np.histogram(valeurs, bins=36, range=(0, 180))
x_array = (bins[:-1] + bins[1:]) / 2

# 定义1全局+3局部高斯拟合函数
def four_gaussian(x, amp1, cen1, sigma1, amp2, cen2, sigma2, amp3, cen3, sigma3, base):
    # 全局高斯（对应5-10区间）
    g1 = amp1 * np.exp(-(x - cen1)**2 / (2 * sigma1**2))
    # 局部高斯1（110-115区间）
    g2 = amp2 * np.exp(-(x - cen2)**2 / (2 * sigma2**2))
    # 局部高斯2（180附近）
    g3 = amp3 * np.exp(-(x - cen3)**2 / (2 * sigma3**2))
    return base + g1 + g2 + g3

# 初始参数设置（根据直方图峰值调整）
p0 = [700, 7, 40, 750, 115, 7, 700, 172, 4, hist.min()]
# 合理设置参数边界：振幅≥0，中心在0-180，σ>0，基线≥0
bounds = (
    [0, 0, 1e-3, 0, 0, 1e-3, 0, 0, 1e-3, 0],
    [np.inf, 180, 100, np.inf, 180, 50, np.inf, 180, 50, np.inf]
)

# 执行拟合
popt, pcov = opt.curve_fit(four_gaussian, x_array, hist, p0=p0, bounds=bounds)
perr = np.sqrt(np.diag(pcov))

# 提取拟合参数
amp1, cen1, sigma1, amp2, cen2, sigma2, amp3, cen3, sigma3, base = popt

# 单个高斯函数，用于绘制和交点计算
def single_gaussian(x, amp, cen, sigma):
    return amp * np.exp(-(x - cen)**2 / (2 * sigma**2)) + base

# 求解两条高斯曲线的交点
def find_intersection(g1_params, g2_params, x_range):
    amp1, cen1, sigma1, base = g1_params
    amp2, cen2, sigma2, base = g2_params
    def func(x):
        return (amp1 * np.exp(-(x - cen1)**2/(2*sigma1**2)) - 
                amp2 * np.exp(-(x - cen2)**2/(2*sigma2**2)))
    return opt.root_scalar(func, bracket=x_range, method='brentq').root

# 计算三个高斯曲线的交点（按中心顺序划分区间）
intersection1 = find_intersection([amp1, cen1,