高斯适合嘈杂和"有趣"的数据集

Question

我有一些看起来像这样的数据(X射线衍射):

我想在这个数据集中加入Gaussian来获得"更宽"部分的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]间隔内的eveything .