时间戳

首页 > 解决方案 > 指数族删失分布的 MLE

python-3.x - 指数族删失分布的 MLE

问题描述

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

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

,

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

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

然而，我发现即使对于简单的分布，例如nelder-mead单纯形算法并不总是收敛，或者它确实收敛但估计的参数与真实参数相去甚远。我在下面附上了我的代码。请注意，可以选择一个分布，并且代码足够通用以适合loc和scale参数，以及可选的形状参数（例如 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)

标签： python-3.xscipy-optimizemle

解决方案

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

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

和

func.logpdf(x, **func_kwargs)

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

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

推荐阅读