Constrained MLE with Python

Question

I'm doing a MLE implementation with Python. My log-likelihood function has 5 parameters to be estimated and two of them has the constraint that they must be between 0 and 1. I'm able to implement a MLE using GenericLikelihoodModel module from the statsmodels package but I don't know how to do this with the constraint. To be specific, my negative log-likelihood function is

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 

and the MLE optimization class is

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) 

And the main is

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

I don't know how to modify my code to incorporate the constraint. I would like alpha and sigma to be between 0 and 1.

Answer 1

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.