MLE用于检查指数族的分布

问题描述

使用 scipy.stats 包，可以很容易地将分布分配到数据中，例如 scipy.stats.expon.fit（）可用于将数据拟合为指数分布。

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

$\text{log likelihood}=\sum_{i=1}^{N}\left\[\log(f(x_{i}\mid\theta))-\log(F(c_{i}\mid\theta)) \right\]$

，

其中

$f(x_{i}\mid \theta)$

是指数族分布的PDF，而

$F(c_{i}\mid \theta)$

是其对应的CDF。

在数学上，我发现对数似然函数在参数空间

$\theta$

中是凸的，因此我的假设是应用 scipy.optimize.minimize 应该相对简单>功能。请注意，在上述对数可能性中，通过采用

$c_{i}\rightarrow\infty$

，我们可以获得传统/未经审查的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,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,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,1000)

res = minimize_log_likelihood(distribution,x,10**8)
print(res)

解决方法

我发现由于数值不正确，收敛性不好。最好是替换

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

使用

func.logpdf(x,**func_kwargs)

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

所有这些都可以很好地与指数分布，正态分布，伽玛分布，chi2分布一起使用。 Beta分布仍然给我带来了问题，但是我认为这还是一些（其他）数值不正确的地方，我将分别进行分析。

mle python-3.x scipy-optimize