【问题标题】:Error for my plot and calculation of FWHM我的绘图和 FWHM 计算错误
【发布时间】:2013-04-30 10:43:42
【问题描述】:

我找到并复制了此代码以从 Finding the full width half maximum of a peak 获取 FWHM(第二个到最后一个答案)。我下面的代码使用我自己的数据。生成的图看起来不对,因为我的数据出现在一侧,而绿色框在另一侧。必须进行哪些更改才能在数据上看到高斯分布?

import numpy as np, scipy.optimize as opt
from pylab import *

def gauss(x, p):
    return 1.0/(p[1]*np.sqrt(2*np.pi))*np.exp(-(x-p[0])**2/(2*p[1]**2))

x = [6711.19873047, 6712.74267578, 6714.28710938, 6715.83544922, \
     6717.38037109, 6718.92919922, 6720.47509766, 6722.02490234, \
     6723.57128906, 6725.11767578, 6726.66845703, 6728.21630859, \
     6729.76757812, 6731.31591797, 6732.86816406, 6734.41699219, \
     6735.96630859, 6737.51953125, 6739.06933594, 6740.62353516, \
     6742.17431641, 6743.72900391]

y = [20.86093712, 23.60984612, 23.079916,   18.17703056, 18.24843597, \
     16.70049095, 19.48906136, 16.7509613,  19.09896088, 32.03007889, \
     54.56513977, 58.76417542, 40.93075562, 24.77710915, 17.68757629, \
     17.60736847, 18.89552498, 17.84486008, 17.49455452, 18.29696465, \
     18.55847931, 19.26465797]

# Fit a guassian
p0 = [0,70]
errfunc = lambda p, x, y: gauss(x, p) - y # Distance to the target function
p1, success = opt.leastsq(errfunc, p0[:], args=(x, y))

fit_mu, fit_stdev = p1

FWHM = 2*np.sqrt(2*np.log(2))*fit_stdev
print "FWHM", FWHM

plot(x,y)
plot(x, gauss(x,p1),lw=3,alpha=.5, color='r')
axvspan(fit_mu-FWHM/2, fit_mu+FWHM/2, facecolor='g', alpha=0.5)
show()

【问题讨论】:

    标签: python python-2.7 numpy matplotlib scipy


    【解决方案1】:

    看起来像最初的猜测,p0,你传递给scipy.optimize.leastsq() 是这里的问题。如果它不够接近(例如不在数据范围内),那么返回的解决方案将没有多大意义。

    为了得到合理的输出,可以初始化p0 = [6730,70],会产生如下的情节:

    但这仍然不是正确的答案,因为您已经跳过了original answer 中的 PDF 规范化步骤。这是适用于您的代码的规范化 sn-p(在适合高斯之前添加它):

    # Renormalize to a proper PDF
    y /= ((max(x) - min(x)) / len(x)) * np.sum(y)
    

    添加它会导致以下情节:

    高斯拟合不是很紧,但至少它有一个合理的均值和标准差,正如我最初在下面描述的那样。


    问题源于最小二乘函数p1的输出,与p0相同。

    换句话说,假设得到的平均值是0 而不是np.mean(x) = 6727.45,那么块将被绘制在一个完全不同的位置。此外,块的宽度绘制为329.67,基于70 的标准偏差,而不是基于np.std(x) = 9.8346.28

    【讨论】:

    • @tcaswell:当然,我只是按照链接的答案。随意修改我的答案。
    猜你喜欢
    • 2015-07-05
    • 2018-08-12
    • 1970-01-01
    • 2023-02-15
    • 1970-01-01
    • 2014-08-14
    • 2017-05-29
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多