С пакетом scipy.stats легко подогнать распределение к данным, например scipy.stats.expon.fit () можно использовать для подгонки данных к экспоненциальному распределению.
Однако я пытаюсь подогнать данные к цензурированному / условному распределению в экспоненциальном семействе. Другими словами, используя MLE, я пытаюсь найти максимум
,
где представляет собой PDF-файл распределения в экспоненциальном семействе, а
- соответствующий ему CDF.
Математически я обнаружил, что функция логарифмического правдоподобия является выпуклой в пространстве параметров , поэтому мой Предполагалось, что применить функцию scipy.optimize.minimize должно быть относительно просто. Обратите внимание на вероятность того, что, взяв
мы получаем традиционную проблему MLE без цензуры.
Однако я считаю, что даже для простых дистрибутивов, например симплексный алгоритм Нелдера-Мида не всегда сходится, или что он сходится, но оценочные параметры далеки от истинных. Я прикрепил свой код ниже. Обратите внимание, что можно выбрать распределение, и что код достаточно общий, чтобы соответствовать параметрам loc и scale, а также необязательной shape параметры (например, бета- или гамма-распределение).
У меня вопрос: что я делаю не так, чтобы получить эти плохие оценки, или иногда возникают проблемы с конвергенцией? Я пробовал несколько алгоритмов, но, к моему удивлению, нет ни одного, который бы легко работал, поскольку проблема выпуклая. Есть ли проблемы с гладкостью, и что мне нужно найти способ использовать якобиан и гессиан в общем виде для этой проблемы?
Есть ли другие методы решения этой проблемы? Первоначально я думал переопределить функцию fit () в классе scipy.stats.rv, чтобы позаботиться о цензуре с помощью CDF, но это казалось довольно громоздким. Но поскольку проблема выпуклая, я бы предположил, что с помощью функции минимизации scipy я смогу легко получить результаты ...
Комментарии и помощь приветствуются!
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import expon, gamma, beta, norm
from scipy.optimize import minimize
from scipy.stats import rv_continuous as rv
def log_likelihood(func: rv, delays, max_delays=10**8, **func_pars)->float:
return np.sum(np.log(func.pdf(delays, **func_pars)+1) - np.log(func.cdf(max_delays, **func_pars)))
def minimize_log_likelihood(func: rv, delays, max_delays):
# Determine number of parameters to estimate (always 'loc', 'scale', sometimes shape parameters)
n_pars = 2 + func.numargs
# Initialize guess (loc, scale, [shapes])
x0 = np.ones(n_pars)
def wrapper(params, *args):
func = args[0]
delays = args[1]
max_delays = args[2]
loc, scale = params[0], params[1]
# Set 'loc' and 'scale' parameters
kwargs = {'loc': loc, 'scale': scale}
# Add shape parameters if existing to kwargs
if func.shapes is not None:
for i, s in enumerate(func.shapes.split(', ')):
kwargs[s] = params[2+i]
return -log_likelihood(func=func, delays=delays, max_delays=max_delays, **kwargs)
# Find maximum of log-likelihood (thus minimum of minus log-likelihood; see wrapper function)
return minimize(wrapper, x0, args=(func, delays, max_delays), options={'disp': True},
method='nelder-mead', tol=1e-8)
# Test code with by sampling from known distribution, and retrieve parameters
distribution = expon
dist_pars = {'loc': 0, 'scale': 4}
x = np.linspace(distribution.ppf(0.0001, **dist_pars), distribution.ppf(0.9999, **dist_pars), 1000)
res = minimize_log_likelihood(distribution, x, 10**8)
print(res)
