问题

我有一组数据，显示辐射的吸收与速度的函数关系。数据显示一个非常明显的低谷，或者如果我们绘制数据的倒数，即吸收度，我们会得到一个明显的峰值。我没有理由不相信这个峰值是一个高斯分布，想要进行拟合以获取这个峰值的方差。因此，我尝试使用scipy.optimize.curve_fit来实现这一目标，无论初始猜测如何，得到的拟合结果都相差很远。

我附上了代码和结果图：

import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import numpy as np

# 以下是数据的读取和处理部分，未翻译
# ...

def func(x, mu, sigma):
    return (1 / (np.sqrt(np.pi * 2) * sigma)) * np.exp(-0.5 * ((x - mu) / sigma) ** 2)

data = np.array(cts)
cts_norm = (data - data.min()) / (data.max() - data.min())
cts_inv = 1 - cts_norm
fit, error = curve_fit(func, vel, cts_inv, p0=[0.2, 0.2])
print(fit)

plt.plot(vel, cts_inv, 'bo')
plt.plot(vel, func(vel, fit[0], fit[1]), 'r')

请注意，以上是代码的翻译，没有包含问题的回答或其他内容。

import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import numpy as np
cts = []
vel = []
file = open("Heisenberg/Mössbauer/Final.lst", "r")
linesArr = file.readlines()
for i in range(210, 260):
    lineList1 = linesArr[i].split()
    cts.append(int(lineList1[1]))
    
    chn = (int(lineList1[0]))
    tempVel = -0.04 * chn + 9.3
    vel.append(tempVel) 

def func (x, mu,sigma):
    return (1 / (np.sqrt(np.pi * 2) *sigma)) * np.exp(-0.5*((x-mu)/sigma)**2)

data = np.array(cts)
cts_norm = (data - data.min())/ (data.max() - data.min())
cts_inv = 1 - cts_norm
fit, error = curve_fit(func, vel, cts_inv, p0=[0.2, 0.2])
print(fit)


plt.plot(vel, cts_inv, 'bo')
plt.plot(vel, func(vel, fit[0],fit[1]),"r")

答案1

得分: 2

问题在于您试图拟合一个正态分布到不是概率分布的数据上！概率分布的积分等于1，但对于您的数据来说并非如此，它可以具有任何幅度。很难将您的数据归一化以满足这一要求。相反，您可以简单地添加一个新的参数来控制正态分布的“幅度”，如下所示：

import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import numpy as np

cts = [0, 0, 0, 0, -1, -2, -5, -10, -5, -2, -1, 0, 0, 0, 0]
vel = np.linspace(-0.75, 1.25, 15)

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

data = np.array(cts)
cts_norm = (data - data.min()) / (data.max() - data.min())
cts_inv = 1 - cts_norm
fit, error = curve_fit(func, vel, cts_inv, p0=[0.2, 0.2, 1])
print(fit)

plt.plot(vel, cts_inv, 'bo')
plt.plot(vel, func(vel, fit[0], fit[1], fit[2]), "r")
plt.show()

（我使用了一些虚拟数据，因为我无法访问您的文件，但这并不重要）

import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import numpy as np

cts = [0, 0, 0, 0, -1, -2, -5, -10, -5, -2, -1, 0, 0, 0, 0]
vel = np.linspace(-0.75, 1.25, 15)


def func(x, mu, sigma, a):
    return a * np.exp(-0.5 * ((x - mu) / sigma) ** 2)  #      << here

data = np.array(cts)
cts_norm = (data - data.min()) / (data.max() - data.min())
cts_inv = 1 - cts_norm
fit, error = curve_fit(func, vel, cts_inv, p0=[0.2, 0.2, 1]) # << here
print(fit)

plt.plot(vel, cts_inv, 'bo')
plt.plot(vel, func(vel, fit[0], fit[1], fit[2]), "r")  #      << and here
plt.show()

答案2

得分: 1

我会为你的模型添加一些灵活性，具体如下。我通过拍摄图像并使用 [此][1] 免费网络服务来检索了你的数据。

```python
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
from scipy.stats import norm

data = np.loadtxt(r"C:\Users\Cristiano\Desktop\data.txt", delimiter=",")
x = data[:, 0]
y = data[:, 1]

def f(x, a, b, mu, sigma):
    return a + b * np.exp(-(x - mu) ** 2 / (2 * sigma ** 2))

popt, pcov = curve_fit(f, x, y)

mean, std = norm.fit(data)

plt.scatter(x, y)
xx = np.linspace(-0.75, 1.25, 1000)
plt.plot(xx, f(xx, *popt))
plt.show()

<details>
<summary>英文:</summary>

I would add a little more flexibility to your model as follows. I retrieved your data by taking a screenshot of the image and using [this][1] free web service.

```python
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
from scipy.stats import norm

data = np.loadtxt(r&quot;C:\Users\Cristiano\Desktop\data.txt&quot;, delimiter=&quot;,&quot;)
x = data[:, 0]
y = data[:, 1]

def f(x, a, b, mu, sigma):
    return a + b * np.exp(-(x - mu) ** 2 / (2 * sigma ** 2))

popt, pcov = curve_fit(f, x, y)

mean, std = norm.fit(data)

plt.scatter(x, y)
xx = np.linspace(-0.75, 1.25, 1000)
plt.plot(xx, f(xx, *popt))
plt.show()