1

私はPythonでMLEの実装を行っています。私の対数尤度関数には、推定する 5 つのパラメーターがあり、そのうちの 2 つには、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

One common approach to get an interior solution that satisfies the constraints is to transform the parameters so the optimization is unconstrained.

For example: A constraint to be in the open interval (0, 1) can be transformed using the Logit function, used for example here:

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

We can use a mulinomial logit for probabilities, parameters that are in (0, 1) and add up to one.

In generalized linear models we use the link functions to impose similar restriction for the predicted mean, see statsmodels/genmod/families/links.py.

If constraints can be binding, then this does not work. Scipy has constrained optimizers but those are not yet connected to the statsmodels LikelihoodModel classes.

Related aside: scipy has a global optimizer, basinhopping, that works pretty well if there are multiple local minima, and it is connected to the LikelihoodModels and can be chosen using the method argument in fit.

于 2016-06-14T13:01:18.670 に答える
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 に答える
0

IMHOこれは数学の質問であり、簡単な答えは、質問を再構築する必要があるということです。

問題に具体的に対処するには、元のモデルから派生した特別なケースのモデルを作成し、その中で固有に定義された制約を使用する必要があります。次に、特別なケース モデルの計算された MLE により、探している推定値が得られます。

しかし、元のモデルでは2つのパラメーターが制約されていなかったため、一般的なケースモデルではなく、制約のある派生モデルの推定は膨らみます。

実際には、MCMC、ANN、ニュートンベースの反復法などのパラメータ推定に使用する方法はすべて、派生モデルと制約付きモデルの推定値を提供します。

于 2016-06-14T10:31:01.833 に答える