Указание и подгонка пользовательского распределения к данным

У меня есть набор данных измерений $Y$, к которому я хочу применить пользовательское распределение, чтобы получить оценку параметров распределения.

Основываясь на знаниях предметной области, я знаю, что процесс, генерирующий $Y$, представляет собой взвешанную смесь двух гамма-распределений:

$Y\sim w\times\Gamma(\alpha_1, \theta_1)+(1-w)\times\Gamma(\alpha_2, \theta_2)$

Я хочу оценить $w$, $\alpha_1$, $\theta_1$, $\alpha_2$ и $\theta_2$.

SciPy имеет различные стандартные распределения, но я хочу более общий способ указать мою пользовательскую смесь.

scipy.optimize.curve_fit не совсем подходит, потому что я хочу оценить распределение по данным, а не напрямую подгонять кривую.

Как это можно сделать с помощью SciPy?

Начните с определения функции, которая строит модель функции плотности вероятности (PDF). В вашем случае это смесь двух гамма-распределений, но это может быть что-то другое:

def two_gamma(y_samples, w, alpha1, scale1, alpha2, scale2):
    #Конструирование текущей оцениваемой PDF
    gamma1_pdf = stats.gamma.pdf(y_samples, alpha1, scale=scale1)
    gamma2_pdf = stats.gamma.pdf(y_samples, alpha2, scale=scale2)
    y_pdf_estimate = w * gamma1_pdf + (1 - w) * gamma2_pdf

    return y_pdf_estimate

Затем вы можете использовать scipy.optimize.minimize, чтобы настроить параметры таким образом, чтобы отрицательная логарифмическая правдоподобность (NLL) была минимизирована. Функция two_gamma_mle ниже использует оцененные параметры для построения PDF и затем возвращает оценку (NLL), насколько хорошо параметры соответствуют данным:

def two_gamma_mle(params):
    #Конструирование текущей оцениваемой PDF (в точках данных Y)
    y_pdf_estimate = two_gamma(Y, *params)

    #Возвращаем среднее отрицательное логарифмическое правдоподобие
    avg_neg_ll = -np.log(y_pdf_estimate + 1e-8).mean()
    return avg_neg_ll

results = optimize.minimize(
  two_gamma_mle,
  x0=initial_guess,
  bounds=bounds
)

w, alpha1, scale1, alpha2, scale2 = results['x']

Технически вы можете использовать curve_fit для подгонки модели смеси к ядерной оценке плотности $Y$ (тем самым оценивая параметры смеси), но подход на основе MLE, описанный выше, более прямой и работает лучше.

Демонстрация

Рассмотрим тестовый набор данных ниже. У нас есть выборки $Y$ и мы хотим оценить параметры смеси $\Gamma$ (зеленый), которая их сгенерировала.

Использование minimize() приводит к хорошей подгонке$^\dagger$ к эталонному распределению:

Converged? True
Estimated params: [ 0.49 19.81  0.51  2.03  0.98]

Код

Импорты, данные для тестирования и визуализация:

import numpy as np
from scipy import stats

import matplotlib.pyplot as plt

#Выборка PDF с использованием выборки обратного преобразования
def sample_pdf(x, pdf_x, n_samples=1000):
    cdf = np.cumsum(pdf_x)

    cdf_vals_sampled = np.random.uniform(cdf.min(), cdf.max(), n_samples)
    cdf_ixs = np.searchsorted(cdf, cdf_vals_sampled)

    return x[cdf_ixs]

def two_gamma(y_samples, w, alpha1, scale1, alpha2, scale2):
    #Конструирование текущей оцениваемой PDF
    gamma1_pdf = stats.gamma.pdf(y_samples, alpha1, scale=scale1)
    gamma2_pdf = stats.gamma.pdf(y_samples, alpha2, scale=scale2)
    y_pdf_estimate = w * gamma1_pdf + (1 - w) * gamma2_pdf

    return y_pdf_estimate

#
# Данные для тестирования: бимодальное Y, состоящее из двух гамма-распределений
#
np.random.seed(0)

W = 0.5
ALPHA1, SCALE1 = 2, 1
ALPHA2, SCALE2 = 20, 0.5

#Создание PDF Y, которая является смесью двух гамма-распределений
y_axis = np.linspace(0, 20, num=5_000)
pdf_y = two_gamma(y_axis, W, ALPHA1, SCALE1, ALPHA2, SCALE2) #w*g1 + (1-w)*g2

#Выборка данных для тестирования
Y = sample_pdf(y_axis, pdf_y, n_samples=5_000)

#
# Визуализация распределения Y и его основных компонентов
#
pdf_gamma1 = stats.gamma(ALPHA1, scale=SCALE1).pdf(y_axis)
pdf_gamma2 = stats.gamma(ALPHA2, scale=SCALE2).pdf(y_axis)

f, axs = plt.subplots(nrows=2, ncols=1, figsize=(8, 4), sharex=True, layout="tight")

ax = axs[0]
ax.plot(y_axis, pdf_gamma1, label="$\gamma_1$")
ax.plot(y_axis, pdf_gamma2, label="$\gamma_2$")
ax.plot(
    y_axis, pdf_y,
    label="Y ~ $w\gamma_1 + (1 - w)\gamma_2$",
    linewidth=7, color="green", alpha=0.25
)

ax.legend(ncols=3, shadow=True)
ax.set(ylabel="PDF", title="PDFs")

ax.spines[['top', 'right']].set_visible(False)
ax.spines.bottom.set_bounds(0, 20)
ax.spines.left.set_bounds(0, 0.2)

#
# Визуализация точек данных
#
ax = axs[1]
ax.plot(Y, range(len(Y)), linestyle="none", marker=".", markersize=1, alpha=0.8, color="black")

ax.set(xlabel="y", ylabel="sample index", title="Samples")

ax.spines[['top', 'right']].set_visible(False)
ax.spines.bottom.set_bounds(0, 20)
ax.spines.left.set_bounds(0, 4000)

Подгонка к данным и вывод результатов:

from scipy import optimize

def two_gamma_mle(params):
    #Конструирование текущей оцениваемой PDF (в точках данных Y)
    y_pdf_estimate = two_gamma(Y, *params)

    #Возвращаем среднее отрицательное логарифмическое правдоподобие
    avg_neg_ll = -np.log(y_pdf_estimate + 1e-8).mean()
    return avg_neg_ll

#
# Подгонка с использованием scipy.optimize.minimize()
#
initial_guess = [1] * 5
bounds = [(0, 1)] + [(0, None)] * 4 #w ограничен (0, 1). Другие параметры (0, inf).

res = optimize.minimize(two_gamma_mle, x0=initial_guess, bounds=bounds)
estimated_params = (w, alpha1, scale1, alpha2, scale2) = res['x']

print('Converged?', res['success'])
print('Estimated params:', estimated_params.round(2))

#
# Просмотр результатов
#
ax = plt.figure(figsize=(8, 3), layout="tight").add_subplot()

ax.plot(y_axis, pdf_y, linewidth=3, label="target PDF")

y_pdf_fit = two_gamma(y_axis, *estimated_params)
ax.plot(y_axis, y_pdf_fit, linestyle="--", color="lawngreen", label="estimated PDF")
ax.legend()
ax.set(xlabel="y", ylabel="PDF", title="Target and MLE-estimated PDF")
ax.spines[['top', 'right']].set_visible(False)
ax.spines.left.set_bounds(0, 0.15)
ax.spines.bottom.set_bounds(0, 20)

$^\dagger$Оцененные параметры для отдельных компонентов поменялись местами по сравнению с истинными значениями, но общая подгонка хорошая. Вы можете ограничить параметры при необходимости и/или задать более точное начальное приближение.