python 上的瑞利分布曲线拟合

Question

我目前正在使用此 PDF 方程编写布朗运动的实验室报告，目的是评估 D： Brownian PDF equation

我正在尝试将其曲线拟合到直方图。但是，每当我绘制我的 curve_fits 时，它都是一条线，并且在直方图上显示不正确。 Example Histogram with bad curve_fit

这是我的代码：

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

# Variables
eta = 1e-3 
ra = 0.95e-6
T = 296.5
t = 0.5

# Random data
r = np.array(np.random.rayleigh(0.5e-6, 500))

# Histogram
plt.hist(r, bins=10, density=True, label='Counts')

# Curve fit
x,y = np.histogram(r, bins=10, density=True)
x = x[2:]
y = y[2:]
bin_width = y[1] - y[2]
print(bin_width)
bin_centers = (y[1:] + y[:-1])/2
err = x*0 + 0.03

def f(r, a):
    return (((1e-6)3*np.pi*r*eta*ra)/(a*T*t))*np.exp(((-3*(1e-6 * r)**2)*eta*ra*np.pi)/(a*T*t))

print(x) # these are flipped for some reason
print(y)

plt.plot(bin_centers, x, label='Fitting this', color='red')


popt, pcov = optimize.curve_fit(f, bin_centers, x, p0 = (1.38e-23), sigma=err, maxfev=1000)

plt.plot(y, f(y, popt), label='PDF', color='orange')
print(popt)

plt.title('Distance vs Counts')
plt.ylabel('Counts')
plt.xlabel('Distance in micrometers')
plt.legend()

是我的curve_fit 有问题吗？还是有我遗漏的潜在问题？

编辑：我将 D 分解为函数中的玻尔兹曼常数 a，这就是为什么 f 中的数字比上面的等式更多。 D 和 Gamma.

我尝试弄乱初始条件并使用1.38e-23 而不是popt 绘制函数，但这确实this (the purple line). 这告诉我f 的方程有问题，但没有问题跳转当我看着它的时候给我看。我错过了什么吗？

编辑 2：我将函数更改为此以简化它并匹配 numpy.random.rayleigh() 分布：

def f(r, a):
    return ((r)/(a))*np.exp((-1*(r)**2)/(2*a))

但这并不能解决curve_fit 是一条具有正斜率的线而不是我感兴趣的任何东西的问题。现在我对问题是什么感到更加困惑。

Answer 1

这里有几件事。我认为 x 和 y 从未翻转过，或者至少当我认为它们没有翻转时，一切似乎都运行良好。我还清理了代码的一些部分，例如，我不确定为什么要调用两个不同的直方图；而且我认为处理参数的单元素元组可能存在问题。此外，对于曲线拟合，初始参数猜测通常需要在大致范围内，所以我也进行了更改。

这是一个适合我的版本：

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

# Random data
r = np.array(np.random.rayleigh(0.5e-6, 500))

# Histogram
hist_values, bin_edges, patches = plt.hist(r, bins=10, density=True, label='Counts')

bin_centers = (bin_edges[1:] + bin_edges[:-1])/2

x = bin_centers[2:]  # not necessary, and I'm not sure why the OP did this, but I'm doing this here because OP does
y = hist_values[2:]

def f(r, a):
    return (r/(a*a))*np.exp((-1*(r**2))/(2*a*a))

plt.plot(x, y, label='Fitting this', color='red')

err = x*0 + 0.03
popt, pcov = optimize.curve_fit(f, x, y, p0 = (1.38e-6,), sigma=err, maxfev=1000)

plt.plot(x, f(x, *popt), label='PDF', color='orange')

plt.title('Distance vs Counts')
plt.ylabel('Counts')
plt.xlabel('Distance in Meters') # Motion seems to be in micron range, but calculation and plot has been done in meters
plt.legend()