高斯拟合噪声和“感兴趣”数据集
高斯拟合噪声和“感兴趣”数据集
提问于 2016-11-10 18:53:18
我有一些数据(X射线衍射)是这样的：
我想要拟合一个高斯到这个数据集，以获得‘更宽’部分的半高宽。双峰在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)
回答 1
Stack Overflow用户
回答已采纳
发布于 2016-11-10 20:33:32
正如OP注释中提到的那样，在存在不想要的数据的情况下约束信号的方法之一是将信号与所需的信号一起建模。当然，这种方法只有在那些受污染的数据有一个有效的模型时才是有效的。对于您提供的数据，可以考虑将以下组件相加的复合模型：
- 一个线性基线，因为你所有的点都不断地从零中偏移。
- 两个窄高斯分量，将模拟双峰特征在您的光谱的中心部分。
- 一个狭窄的高斯分量。这才是你真正想要约束的。
一旦您通过对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]间隔内的每一件事情。
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/40534986
复制相关文章
相似问题
腾讯云开发者
