python 上的瑞利分布 Curve_fit

Rayleigh distribution Curve_fit on python

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

我正在尝试 curve_fit 将其转换为直方图。但是,每当我绘制我的 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 and 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 是一条具有正斜率的线而不是我感兴趣的任何远程问题的问题。现在我对问题是什么感到更加困惑。

这里有几件事。我不认为 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()