python - 噪声和 'interesting' 数据集上的高斯拟合

Question

我有一些数据(X 射线衍射)如下所示:

我想对这个数据集进行高斯拟合，以获得“更宽”部分的 FWHM。 7 度左右的双峰不是重要信息，来自不需要的来源。

为了让自己更清楚，我想要这样的东西(我用油漆画的:)):

我尝试使用以下代码在 python 中编写脚本:

import math
from pylab import *
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
data2=np.loadtxt('FWHM.spc')    
x2,y2=data2[:,0],data2[:,7]
plt.title('Full Width Half Max of 002 Peak')
plt.plot(x2, y2, color='b')
plt.xlabel('$\\theta$', fontsize=10)
plt.ylabel('Intensity', fontsize=10)
plt.xlim([3,11])
plt.xticks(np.arange(3, 12, 1), fontsize=10)
plt.yticks(fontsize=10)

def func(x, a, x0, sigma):
    return a*np.exp(-(x-x0)**2/(2*sigma**2))

mean = sum(x2*y2)/sum(y2)
sigma2 = sqrt(abs(sum((x2-mean)**2*y2)/sum(y2)))
popt, pcov = curve_fit(func, x2, y2, p0 = [1, mean, sigma2])
ym = func(x2, popt[0], popt[1], popt[2])
plt.plot(x2, ym, c='r', label='Best fit')
FWHM = round(2*np.sqrt(2*np.log(2))*popt[2],4)
axvspan(popt[1]-FWHM/2, popt[1]+FWHM/2, facecolor='g', alpha=0.3, label='FWHM = %s'%(FWHM))
plt.legend(fontsize=10)
plt.show()

我得到以下输出:

显然，这与我们想要的相去甚远。有没有人对我如何实现这一点有一些建议？

(我在这里附上了数据:https://justpaste.it/109qp)

Answer 1

如 OP 评论中所述，在存在不需要的数据时限制信号的方法之一是将其与所需信号一起建模。当然，这种方法仅在那些污染数据随时可用的有效模型时才有效。对于您提供的数据，可以考虑对以下组件求和的复合模型:

1. 线性基线，因为您的所有点都不断偏离零。
2. 两个窄高斯分量，用于模拟光谱中心部分的双峰特征。
3. 窄高斯分量。这是您实际上要限制的对象。

一旦您将合理的起始猜测传递给 curve_fit，就可以同时拟合所有四个分量(双峰计数两次):

def composite_spectrum(x, # data
                       a, b, # linear baseline
                       a1, x01, sigma1, # 1st line
                       a2, x02, sigma2, # 2nd line
                       a3, x03, sigma3 ): # 3rd line
    return (x*a + b + func(x, a1, x01, sigma1)
                    + func(x, a2, x02, sigma2)
                    + func(x, a3, x03, sigma3))

guess = [1, 200, 1000, 7, 0.05, 1000, 6.85, 0.05, 400, 7, 0.6]

popt, pcov = curve_fit(composite_spectrum, x2, y2, p0 = guess)
plt.plot(x2, composite_spectrum(x2, *popt), 'k', label='Total fit')
plt.plot(x2, func(x2, *popt[-3:])+x2*popt[0]+popt[1], c='r', label='Broad component')
FWHM = round(2*np.sqrt(2*np.log(2))*popt[10],4)
plt.axvspan(popt[9]-FWHM/2, popt[9]+FWHM/2, facecolor='g', alpha=0.3, label='FWHM = %s'%(FWHM))
plt.legend(fontsize=10)
plt.show()

在无法正确建模不需要的源的情况下，可以按照 Mad Physicist 的建议屏蔽不需要的不连续性.对于最简单的情况，您甚至可以简单地屏蔽掉 [6.5; 中的所有内容； 7.4]间隔。