指数族删失分布的 MLE答案

【问题标题】：MLE for censored distributions of the exponential family指数族删失分布的 MLE
【发布时间】：2020-10-06 15:42:00
【问题描述】：

使用 scipy.stats 包，可以直接将分布拟合到数据，例如scipy.stats.expon.fit() 可用于将数据拟合成指数分布。

但是，我正在尝试将数据拟合到指数族中的删失/条件分布。换句话说，使用 MLE，我试图找到最大值

其中是指数族分布的 PDF，是其对应的 CDF。

在数学上，我发现对数似然函数在参数空间中是凸的，所以我的假设是应用 scipy.optimize.minimize 函数应该相对简单.请注意，在上述对数似然中，通过采用，我们得到了传统/未经审查的 MLE 问题。

但是，我发现即使对于简单的发行版，例如nelder-mead 单纯形算法并不总是收敛，或者它确实收敛但估计的参数与真实参数相去甚远。我在下面附上了我的代码。请注意，可以选择分布，并且代码足够通用以适应 loc 和 scale 参数，以及可选的 shape参数（例如 Beta 或 Gamma 分布）。

我的问题是：我做错了什么来获得这些错误的估计，或者有时会出现收敛问题？我已经尝试了一些算法，但没有一个可以轻松工作，令我惊讶的是，问题是凸的。是否存在平滑问题，并且我需要找到一种方法以通用方式使用 Jacobian 和 Hessian 来解决此问题？

还有其他方法可以解决这个问题吗？最初我想重写 scipy.stats.rv 类中的 fit() 函数来处理 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)

【问题讨论】：

补充说明：一般来说，使用 MLE 拟合位置参数非常困难。保持 'loc' 参数固定解决了我所有的问题，并且由于上述问题是凸的，我现在确实得到了几乎任意初始值 x0 的正确解决方案。

标签： python-3.x scipy-optimize mle

【解决方案1】：

我发现由于数值不准确，收敛性很差。最好是更换

np.log(func.pdf(x, **func_kwargs))

与

func.logpdf(x, **func_kwargs)

这导致参数的正确估计。 CDF 也是如此。 scipy的文档也表明后者的数值精度表现更好。

这一切都适用于指数、正态、伽玛、chi2 分布。 Beta 分布仍然给我带来问题，但我认为这又是一些（其他）数字不准确的问题，我将单独分析。

【讨论】：