使用lmfit进行穆斯堡尔谱曲线拟合

huangapple
go评论75阅读模式
英文:

Curve fitting with lmfit for Mossbauer spectroscopy

问题

我对lmfit还比较新，经常遇到这个错误。我的数据集相当大，有大约27000个数据点，我收集到了这些数据：

[我的错误消息](https://i.stack.imgur.com/uh87r.png)

```python

import numpy as np
import matplotlib.pyplot as plt
from lmfit import Model
import pandas as pd

plt.rcParams['figure.dpi'] = 300
plt.rcParams['savefig.dpi'] = 300

# 使用负号定义洛伦兹函数
def lorentzian(x, amp, cen, wid):
    return (-amp * wid ** 2 / ((x - cen) ** 2 + wid ** 2))

# 定义两个洛伦兹峰的总和
def multi_lorentzian(x, a1, c1, w1, a2, c2, w2):
    return (lorentzian(x, a1, c1, w1) + lorentzian(x, a2, c2, w2))

csv_names = ['Channel', 'Energy', 'Counts']
fp2 = "C:\\Users\\mfern\\Downloads\\Data_3_8.csv"

df = pd.read_csv(fp2, skiprows=23, names=csv_names).dropna()  # 删除包含NaN值的行

xdata = df['Channel']
ydata = df['Counts']

model = Model(multi_lorentzian)

params = model.make_params(a1=0, c1=0, w1=0, a2=0, c2=0, w2=0)
params['a1'].min = 0
params['a2'].min = 0
params['w1'].min = 0
params['w2'].min = 0

result = model.fit(ydata, params=params, x=xdata)

# 打印拟合报告
print(result.fit_report())

# 绘制数据和拟合曲线
plt.scatter(xdata, ydata, label='data', s=2, c='b')
plt.plot(xdata, result.best_fit, 'r-', label='fit')
plt.legend()

所以如果有人可以给我一些启示，我将不胜感激 [这是我的数据的样子]
(https://i.stack.imgur.com/Kg0yq.png)

这是我的数据集链接（使用谷歌驱动器链接）：
1: https://drive.google.com/file/d/1ho5b1iGV9_wft9h_e6xP9KRtQgGCufcD/view?usp=sharing


<details>
<summary>英文:</summary>
Im pretty new to lmfit and I keep running into this error. My dataset is pretty large ~27000 points of data from I gather:
[my error msg](https://i.stack.imgur.com/uh87r.png)

import numpy as np
import matplotlib.pyplot as plt
from lmfit import Model
import pandas as pd

plt.rcParams['figure.dpi'] = 300
plt.rcParams['savefig.dpi'] = 300

Define the Lorentzian function with a negative sign

def lorentzian(x, amp, cen, wid):
return (-amp * wid ** 2 / ((x - cen) ** 2 + wid ** 2))

Define the sum of two Lorentzian peaks

def multi_lorentzian(x, a1, c1, w1, a2, c2, w2):
return (lorentzian(x, a1, c1, w1) + lorentzian(x, a2, c2, w2))

csv_names = ['Channel', 'Energy', 'Counts']
fp2 = "C:\Users\mfern\Downloads\Data_3_8.csv"

df = pd.read_csv(fp2, skiprows=23, names=csv_names).dropna() # Remove rows with NaN values

xdata = df['Channel']
ydata = df['Counts']

model = Model(multi_lorentzian)

params = model.make_params(a1=0, c1=0, w1=0, a2=0, c2=0, w2=0)
params['a1'].min = 0
params['a2'].min = 0
params['w1'].min = 0
params['w2'].min = 0

result = model.fit(ydata, params=params, x=xdata)

Print the fit report

print(result.fit_report())

Plot the data and the fit curve

plt.scatter(xdata, ydata, label='data', s=2, c='b')
plt.plot(xdata, result.best_fit, 'r-', label='fit')
plt.legend()


So if there is anyone that can enlighten me, I would greatly appreciate it [this is what my data looks like btw]
(https://i.stack.imgur.com/Kg0yq.png)
AND this is link to my data set (with google drive link): 
[1]: https://drive.google.com/file/d/1ho5b1iGV9_wft9h_e6xP9KRtQgGCufcD/view?usp=sharing
</details>
# 答案1
**得分**: 1
以下是您提供的代码的中文翻译：
"我不确定，但我认为您所看到的错误可能来自您正在使用的数据集。您没有提供数据，所以很难确定。
我认为您对您的模型的一般方法是正确的，但缺少一些功能。首先，您的洛伦兹分布在远离它们的中心点时将具有零值，而您的莫斯堡尔数据在大约为25000附近。
更重要的是，您绝对需要为参数提供一些相对现实的初始值，而这些值不应位于您设置的边界上。如果不明显，将所有值都设置为零将会产生一个与您的数据完全不同的模型曲线。我不能再强调这一点了：您绝对必须提供初始值，这些初始值应该足够接近，以便它们可以通过小幅改变而在“数据-模型”结果中产生明显变化。
由于您没有提供数据，我制作了一些看起来有点接近的数据，并使用了内置的lmfit模型，类似于以下内容，这可能是一个不错的起点：
```python
#!/usr/bin/env python
import matplotlib.pyplot as plt
import numpy as np
from lmfit.models import ConstantModel, LorentzianModel
from lmfit.lineshapes import lorentzian
npts = 2001
velo = np.linspace(0, 1000.0, npts)
signal = (26900 + np.random.normal(size=npts, scale=100)
- 80000. * lorentzian(velo, center=240, sigma=6)
- 79000. * lorentzian(velo, center=790, sigma=7) )
# 模型 = 常数模型 + 2个洛伦兹分布
model = ConstantModel() + LorentzianModel(prefix='p1_') + LorentzianModel(prefix='p2_')
# 猜测初始值
params = model.make_params(c=25000,
p1_center=250, p1_sigma=10, p1_amplitude=-10000,
p2_center=750, p2_sigma=10, p2_amplitude=-10000)
params['p1_amplitude'].max = 0  # 强制p1_amplitude < 0
params['p2_amplitude'].max = 0
out = model.fit(signal, params, x=velo)
print(out.fit_report())
plt.plot(velo, signal, 'o', label='data')
plt.plot(velo, out.best_fit, label='fit')
plt.legend()
plt.show()

这将生成以下报告：

[[模型]]
((模型(常数) + 模型(洛伦兹分布, 前缀='p1_')) + 模型(洛伦兹分布, 前缀='p2_'))
[[拟合统计信息]]
# 拟合方法   = leastsq
# 函数评估次数   = 109
# 数据点数   = 2001
# 变量数   = 7
卡方值         = 19912468.9
减小的卡方值 = 9986.19301
阿卡伊克信息准则 = 18434.1141
贝叶斯信息准则 = 18473.3239
R方 = 0.96651410
[[变量]]
c:             26897.4478 +/- 2.43041838 (0.01%) (init = 25000)
p1_amplitude: -80171.1149 +/- 638.944430 (0.80%) (init = -10000)
p1_center:     239.908167 +/- 0.04587653 (0.02%) (init = 250)
p1_sigma:      5.99624015 +/- 0.06621609 (1.10%) (init = 10)
p2_amplitude: -78808.0930 +/- 694.381394 (0.88%) (init = -10000)
p2_center:     789.912879 +/- 0.05883249 (0.01%) (init = 750)
p2_sigma:      6.99732305 +/- 0.08516581 (1.22%) (init = 10)
p1_fwhm:       11.9924803 +/- 0.13243219 (1.10%) == '2.0000000*p1_sigma'
p1_height:    -4255.87684 +/- 32.5617804 (0.77%) == '0.3183099*p1_amplitude/max(1e-15, p1_sigma)'
p2_fwhm:       13.9946461 +/- 0.17033163 (1.22%) == '2.0000000*p2_sigma'
p2_height:    -3584.99901 +/- 30.1429487 (0.84%) == '0.3183099*p2_amplitude/max(1e-15, p2_sigma)'
[[相关性]] (未报告的相关性 < 0.100)
C(p2_amplitude, p2_sigma) = -0.7230
C(p1_amplitude, p1_sigma) = -0.7211
C(c, p2_amplitude)        = -0.2989
C(c, p1_amplitude)        = -0.2800
C(c, p2_sigma)            = +0.2133
C(c, p1_sigma)            = +0.1998

以及如下的绘图：

使用lmfit进行穆斯堡尔谱曲线拟合

英文:

I'm not certain, but I think the error you are seeing may come from the dataset you are using. You didn't provide data, so it's hard to know for sure.

I think the general approach to the model you have is OK, but missing some features. For one thing, your Lorentzian's will have a value of zero far from their centroids, whereas your Mossbauer data is around 25000.

More importantly, you definitely need to give somewhat realistic initial values for the parameters, and those values should not at the boundary you set. If it is not obvious, setting all values to zero will give a model curve that looks nothing like your data. I cannot stress this strongly enough: you absolutely must provide initial values that are reasonably close such that they can be changed by a small amount and give a noticeable change in the result of "data-model".

Since you didn't provide data, I made some up that looks sort of close ;), and used the builtin lmfit models, something like this, which might be a good start:

#!/usr/bin/env python
import matplotlib.pyplot as plt
import numpy as np
from lmfit.models import ConstantModel, LorentzianModel
from lmfit.lineshapes import lorentzian
npts = 2001
velo = np.linspace(0, 1000.0, npts)
signal = (26900 + np.random.normal(size=npts, scale=100)
- 80000. * lorentzian(velo, center=240, sigma=6)
- 79000. * lorentzian(velo, center=790, sigma=7) )
# model = constant + 2 Lorentzians
model = ConstantModel() + LorentzianModel(prefix=&#39;p1_&#39;) + LorentzianModel(prefix=&#39;p2_&#39;)
# guess initial values
params = model.make_params(c=25000,
p1_center=250, p1_sigma=10, p1_amplitude=-10000,
p2_center=750, p2_sigma=10, p2_amplitude=-10000)
params[&#39;p1_amplitude&#39;].max = 0  # enforce p1_amplitude &lt; 0
params[&#39;p2_amplitude&#39;].max = 0
out = model.fit(signal, params, x=velo)
print(out.fit_report())
plt.plot(velo, signal, &#39;o&#39;, label=&#39;data&#39;)
plt.plot(velo, out.best_fit, label=&#39;fit&#39;)
plt.legend()
plt.show()

This will give a report of

[[Model]]
((Model(constant) + Model(lorentzian, prefix=&#39;p1_&#39;)) + Model(lorentzian, prefix=&#39;p2_&#39;))
[[Fit Statistics]]
# fitting method   = leastsq
# function evals   = 109
# data points      = 2001
# variables        = 7
chi-square         = 19912468.9
reduced chi-square = 9986.19301
Akaike info crit   = 18434.1141
Bayesian info crit = 18473.3239
R-squared          = 0.96651410
[[Variables]]
c:             26897.4478 +/- 2.43041838 (0.01%) (init = 25000)
p1_amplitude: -80171.1149 +/- 638.944430 (0.80%) (init = -10000)
p1_center:     239.908167 +/- 0.04587653 (0.02%) (init = 250)
p1_sigma:      5.99624015 +/- 0.06621609 (1.10%) (init = 10)
p2_amplitude: -78808.0930 +/- 694.381394 (0.88%) (init = -10000)
p2_center:     789.912879 +/- 0.05883249 (0.01%) (init = 750)
p2_sigma:      6.99732305 +/- 0.08516581 (1.22%) (init = 10)
p1_fwhm:       11.9924803 +/- 0.13243219 (1.10%) == &#39;2.0000000*p1_sigma&#39;
p1_height:    -4255.87684 +/- 32.5617804 (0.77%) == &#39;0.3183099*p1_amplitude/max(1e-15, p1_sigma)&#39;
p2_fwhm:       13.9946461 +/- 0.17033163 (1.22%) == &#39;2.0000000*p2_sigma&#39;
p2_height:    -3584.99901 +/- 30.1429487 (0.84%) == &#39;0.3183099*p2_amplitude/max(1e-15, p2_sigma)&#39;
[[Correlations]] (unreported correlations are &lt; 0.100)
C(p2_amplitude, p2_sigma) = -0.7230
C(p1_amplitude, p1_sigma) = -0.7211
C(c, p2_amplitude)        = -0.2989
C(c, p1_amplitude)        = -0.2800
C(c, p2_sigma)            = +0.2133
C(c, p1_sigma)            = +0.1998

and a plot of

使用lmfit进行穆斯堡尔谱曲线拟合

huangapple
  • 本文由 发表于 2023年5月25日 09:04:24
  • 转载请务必保留本文链接：https://go.coder-hub.com/76328262.html
  • curve-fitting
  • physics
  • python
  • python-2.7
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定