Как я могу смоделировать кривую асимметричных пиков, используя scipy.stats.beta?
Я пытаюсь интегрировать хроматографические пики. Мне нужно иметь возможность моделировать форму пика, чтобы последовательно определять, когда пики начинаются и заканчиваются. Я могу моделировать гауссовы пики, но многие пики асимметричны, как в этом примере данных, и их лучше моделировать с помощью бета-распределения. Я могу подогнать к данным гауссовские модели, но не могу получить бета-модели, соответствующие пику. Что я здесь делаю не так? Я читал много сообщений SO с похожими вопросами, но ни один из них не демонстрирует, как построить производную модель на основе необработанных данных, чтобы показать соответствие. Например:
Scipy — Как подогнать этот бета-дистрибутив с помощью Python Scipy Curve Fit
Пытаюсь построить график бета-распределения с помощью Python
Как правильно построить PDF-файл бета-функции в scipy.stats
Вот мой код:
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()
Ответы 2
Я предполагаю, что dat = np.repeat(xx, yy): xx — это значения наблюдений, а yy — это частоты.
Я также предполагаю, что вы должны использовать beta для подгонки дискретных данных, как указано (хотя это не похоже на лучший выбор распределения), и что вы хотите использовать оптимизацию максимального правдоподобия для подгонки (поскольку это то, что ваш код пытается сделать прямо сейчас).
Ваш код может работать как есть, но вы рисуете PDF-файл, который интегрируется до 1, на гистограмме, которая интегрируется до np.sum(yy). У них будут очень разные области, поэтому их нельзя сравнивать. Вам нужно либо умножить PDF на np.sum(yy), либо нормализовать гистограмму. Я сделаю последнее, поскольку у меня нет библиотек, которые вы используете.
Кроме того, я предоставлю разумные предположения для loc и scale, основываясь на том факте, что в np.min(xx)-1 и np.max(xx)+1 было ноль наблюдений).
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 to speed up optimization,
# 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)
Если у вас возникли проблемы с оптимизацией, которая не сходится для полного набора данных, вы можете исправить loc и scale.
# _fix_ loc and scale this time (floc=loc, fscale=scale)
a, b, loc, scale = stats.beta.fit(dat, floc=loc, fscale=scale)
plt.hist(dat, bins=np.arange(np.min(xx)-0.5, np.max(xx)+0.5), density=True)
pdf = stats.beta(a, b, loc, scale).pdf(xx)
plt.plot(xx, pdf)
Если вы хотите построить график в исходном масштабе, а не в нормализованном PDF-файле:
plt.plot(xx, yy, label='data')
plt.plot(xx, pdf*np.sum(yy), label='beta fit')
plt.legend()
Когда пики асимметричны, особенно в хроматографии, можно использовать классическую модель Экспоненциально модифицированный гауссов Пик.
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats, optimize
def model(x, A, sigma, x0, xG):
z = (x - xG) / sigma - sigma / x0
return A / x0 * np.exp(-0.5 * np.power(sigma / x0, 2) - (x - xG) / x0) * stats.norm.cdf(z)
y = np.array([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])
x = np.arange(y.size)
popt, pcov = optimize.curve_fit(model, x, y, p0=[1e9, 1, 5, 1], bounds=(0., np.inf))
Привет, jlandercy, спасибо за ответ. Ваша модель выглядит лучше, чем бета-дистрибутив. Как замечает Мэтт Хаберланд ниже, бета-распределение — не лучшая модель для этих данных.
Спасибо, Мэтт. Это отличный ответ. Теперь, когда я вижу график, я согласен, что бета-распределение не лучшим образом соответствует данным. Однако мне нужно было это увидеть, чтобы понять это.