1

我正在使用 Python 进行 MLE 实现。我的对数似然函数有 5 个要估计的参数,其中两个具有必须介于 0 和 1 之间的约束。我可以使用 statsmodels 包中的 GenericLikelihoodModel 模块实现 MLE,但我不知道如何使用约束执行此操作。具体来说,我的负对数似然函数是

def ekop_ll(bs,alpha,mu,sigma,epsilon_b,epsilon_s):
    ll=[]
    for bsi in bs:
        b=bsi[0]
        s=bsi[1]
        part1 = (1-alpha)*stats.poisson.pmf(b,epsilon_b)*stats.poisson.pmf(s,epsilon_s)
        part2 = alpha*sigma*stats.poisson.pmf(b,epsilon_b)*stats.poisson.pmf(s,mu+epsilon_s)
        part3 = alpha*(1-sigma)*stats.poisson.pmf(b,mu+epsilon_b)*stats.poisson.pmf(s,epsilon_s)
        li = part1+part2+part3
        if part1+part2+part3 == 0:
            li = 10**(-100)
        lli = np.log(li)
        ll.append(lli)
    llsum = -sum(ll)
    return llsum 

和 MLE 优化类是

class ekop(GenericLikelihoodModel):
    def __init__(self,endog,exog=None,**kwds):
        if exog is None:
            exog = np.zeros_like(endog)
        super(ekop,self).__init__(endog,exog,**kwds)
    def nloglikeobs(self,params):
        alpha = params[0]
        mu = params[1]
        sigma = params[2]
        epsilon_b = params[3]
        epsilon_s = params[4]
        ll = ekop_ll(self.endog,alpha=alpha,mu=mu,sigma=sigma,epsilon_b=epsilon_b,epsilon_s=epsilon_s)
        return ll
    def fit(self, start_params=None, maxiter=10000, maxfun=5000, **kwds):
         if start_params == None:
             # Reasonable starting values
             alpha_default = 0.5
             mu_default = np.mean(self.endog)
             sigma_default = 0.5
             epsilon_b_default = np.mean(self.endog)
             epsilon_s_default = np.mean(self.endog)
             start_params =[alpha_default,mu_default,sigma_default,epsilon_b_default,epsilon_s_default]
         return super(ekop, self).fit(start_params=start_params,
                                      maxiter=maxiter, maxfun=maxfun,
                                      **kwds) 

主要是

if __name__ == '__main__':
    bs = #my data#
    mod = ekop(bs)
    res = mod.fit() 

我不知道如何修改我的代码以合并约束。我希望 alpha 和 sigma 介于 0 和 1 之间。

4

3 回答 3

2

获得满足约束的内部解决方案的一种常见方法是转换参数,以便优化不受约束。

例如:在开区间 (0, 1) 中的约束可以使用 Logit 函数进行转换,例如这里使用:

https://github.com/statsmodels/statsmodels/blob/master/statsmodels/miscmodels/count.py#L243

我们可以对概率使用多项式 logit,参数在 (0, 1) 中并且加起来为 1。

在广义线性模型中,我们使用链接函数对预测均值施加类似的限制,请参见 statsmodels/genmod/families/links.py。

如果约束可以绑定,那么这不起作用。Scipy 限制了优化器,但尚未连接到 statsmodels LikelihoodModel 类。

相关的:scipy 有一个全局优化器,basinhopping,如果有多个局部最小值,它可以很好地工作,并且它连接到 LikelihoodModels 并且可以使用 fit 中的方法参数进行选择。

于 2016-06-14T13:01:18.670 回答
0

恕我直言,这是一个数学问题,简单的答案是你应该重塑你的问题。

为了具体解决这个问题 - 您应该创建一个从原始模型派生的特例模型,其中固有定义的约束。然后为特殊情况模型计算的 MLE 将为您提供您正在寻找的估计。

但是 - 对于具有约束的派生模型的估计会膨胀,而不是对于一般案例模型,因为在原始模型中,两个参数没有受到约束。

实际上,您将用于参数估计的任何方法,如 MCMC、ANN、基于牛顿的迭代方法和其他方法,它们都将为您提供派生和约束模型的估计。

于 2016-06-14T10:31:01.833 回答
0

事实上,这是一道数学题,而不是编程题。我设法通过将带有约束的参数(即 alpha 和 sigma)转换为 alpha_hat 和 sigma_hat 来解决这个问题,

    alpha = 1/(1+np.exp(-alpha_hat))
    sigma = 1/(1+np.exp(-sigma_hat))

这样我们就可以不受约束地估计 alpha_hat 和 sigma_hat。

于 2016-06-16T02:18:47.273 回答