如何使用 scipy.stats.beta 对不对称峰的曲线进行建模？

Question

我正在尝试积分色谱峰。我需要能够对峰形状进行建模，以一致地确定峰的开始和结束时间。我可以对高斯峰值进行建模，但许多峰值是不对称的，如本示例数据所示，并且使用 Beta 分布进行建模会更好。我可以将高斯模型拟合到数据，但我似乎无法获得拟合峰值的 Beta 模型。我在这里做错了什么？我读过很多有类似问题的帖子，但没有一个展示如何在原始数据上绘制派生模型以显示拟合度。 例如：

Scipy - 如何使用 Python Scipy 曲线拟合来拟合此 beta 分布

努力使用 Python 绘制 Beta 分布图

如何在 scipy.stats 中正确绘制 beta 函数的 pdf

这是我的代码：

from scipy.stats import beta
import plotly.graph_objects as go
import peakutils

# raw data 
yy = [128459, 1822448, 10216680, 24042041, 30715114, 29537797, 25022446, 18416199, 14138783, 12116635, 9596337, 7201602, 5668133, 4671416, 3920953, 3259980, 2756295, 2326780, 2095209, 1858894, 1646824, 1375129, 1300799, 1253879, 1086045, 968363, 932041, 793707, 741462, 741593]
xx = list(range(0,len(yy)))
fig = go.Figure()
fig.add_trace(go.Scatter(x = xx, y = yy, name = 'raw', mode = 'lines'))

# Gaussian model
gamp, gcen, gwid = peakutils.peak.gaussian_fit(xx, yy, center_only= False) # this calculates the Amplitude, Center & Width of a guassian peak
print(gamp, gcen, gwid)
gauss_model = [peakutils.peak.gaussian(x, gamp, gcen, gwid) for x in xx]
fig.add_trace(go.Scatter(x = xx, y = gauss_model, mode = 'lines', name = 'Gaussian'))

# Beta distribution model
aa, bb, loc, scale = beta.fit(dat)
print('beta fit:', beta_fit)
beta_model = beta.pdf(xx, aa, bb, loc, scale)
fig.add_trace(go.Scatter(x = xx, y = beta_model, mode = 'lines', name = 'Beta'))

fig.show()

Answer 1

我假设

dat = np.repeat(xx, yy)

-

xx

是观测值，

yy

是频率。

我还假设您必须使用

beta

来拟合离散数据。 （不过，它看起来并不是最好的发行版选择。）

您的代码可能按原样工作，但您正在根据积分为

np.sum(yy)

的直方图绘制 PDF（积分为 1）。它们的领域非常不同，因此不具有可比性。您需要将 PDF 乘以

np.sum(yy)

或标准化直方图。我会做后者，因为我没有你正在使用的库。

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

yy = [128459, 1822448, 10216680, 24042041, 30715114, 29537797, 25022446, 18416199, 14138783, 12116635, 9596337, 7201602, 5668133, 4671416, 3920953, 3259980, 2756295, 2326780, 2095209, 1858894, 1646824, 1375129, 1300799, 1253879, 1086045, 968363, 932041, 793707, 741462, 741593]
xx = list(range(0,len(yy)))
dat = np.repeat(xx, yy)  # asssume this is how `dat` is defined

# estimate `loc` and `scale` manually
loc = np.min(xx) - 1 
scale = (np.max(xx) + 1) - loc

# fit a subset of the data, guessing `loc` and `scale`
a, b, loc, scale = stats.beta.fit(dat[::1000], loc=loc, scale=scale)

# make bins to contain the observations nicely
bins = np.arange(np.min(xx)-0.5, np.max(xx)+0.5)

# plot _normalized_ histogram
plt.hist(dat, bins=bins, density=True)

# plot pdf
pdf = stats.beta(a, b, loc, scale).pdf(xx)
plt.plot(xx, pdf)