CODEWITHSUNDEEP
pythonoptimizationscipystatisticsminimize
J. Grünenwald
J. Grünenwald

Reputation: 115

scipy.optimize.minimize choosing parameters that defy constraints

I am running scipy.optimize.minimize trying to maximize the likelihood for left-truncated data on a Gompertz distribution. Since the data is left-truncated at 1, I get this likelihood:

# for a single point x_i, the left-truncated log-likelihood is:
#   ln(tau) + tau*(ln(theta) - ln(x_i)) - (theta / x_i) ** tau - ln(x_i) - ln(1 - exp(-(theta / d) ** tau))
def to_minimize(args, data, d=1):
  theta, tau = args
  if tau <= 0 or theta <= 0 or theta / d < 0 or np.exp(-(theta / d) ** tau) >= 1:
    print('ERROR')
  term1 = len(data) * (np.log(tau) + tau * np.log(theta) - np.log(1 - np.exp(-(theta / d) ** tau)))
  term2 = 0
  for x in data:
    term2 += (-(tau + 1) * np.log(x)) - (theta / x) ** tau
  return term1 + term2

This will fail in all instances where the if statement is true. In other words, tau and theta have to be strictly positive, and theta ** tau must be sufficiently far away from 0 so that np.exp(-theta ** tau) is "far enough away" from 1, since otherwise the logarithm will be undefined.

These are the constraints which I thus defined. I used the notation with a dict instead of a NonlinearConstraints object since it seems that this methods accepts strict inequality (np.exp(-x[0] ** x[1]) must be strictly less than 1). Maybe I have misunderstood the documentation on this.

def constraints(x):
  return [1 - np.exp(-(x[0]) ** x[1])]

To maximize the likelihood, I minimize the negative likelihood.

opt = minimize(lambda args: -to_minimize(args, data),
               x0=np.array((1, 1)), 
               constraints={'type': 'ineq', 'fun': constraints},
               bounds=np.array([(1e-15, 10), (1e-15, 10)]))

As I take it, the two arguments should then never be chosen in a way such that my code fails. Yet, the algorithm tries to move theta very close to its lower bound and tau very close to its upper bound so that the logarithm becomes undefined.

What makes my code fail?

Upvotes: 0

Views: 666

Answers (1)

joni
joni

Reputation: 7157

Both forms of constraints, i.e. NonlinearConstraint and dict constraints don't support strict inequalities. Typically, one therefore uses g(x) >= c + Ɛ to model the strict inequality g(x) > c, where Ɛ is a sufficiently small number.

Note also that it is not guaranteed that each iteration lies inside the feasible region. Internally, most of the methods try to bring it back into the feasible region by a simple clipping of the bounds. In cases where this doesn't work, you can try NonlinearConstraints keep_feasible option and then use the trust-constr method:

import numpy as np
from scipy.optimize import NonlinearConstraint, minimize

def con_fun(x):
    return 1 - np.exp(-(x[0]) ** x[1])

# 1.0e-8 <= con_fun <= np.inf
con = NonlinearConstraint(con_fun, 1.0e-8, np.inf, keep_feasible=True)

x0 = np.array((1., 1.))
bounds = np.array([(1e-5, 10), (1e-5, 10)])

opt = minimize(lambda args: -to_minimize(args, data),
               x0=x0, constraints=(con,),
               bounds=bounds, method="trust-constr")

Upvotes: 1

Related Questions