scipy.optimize.minimize seleccionando parámetros que desafíen las restricciones.
Estoy ejecutando scipy.optimize.minimize tratando de maximizar la verosimilitud para datos truncados a la izquierda en una distribución de Gompertz. Dado que los datos están truncados a la izquierda en 1, obtengo la siguiente verosimilitud:
Para un único punto x_i, la log-verosimilitud truncada a la izquierda es:
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
Esto fallará en todas las instancias donde la declaración “if” sea verdadera. En otras palabras, tau y theta deben ser estrictamente positivos, y theta*tau debe estar lo suficientemente lejos de 0 para que np.exp(-theta**tau) esté “suficientemente lejos” de 1, ya que de lo contrario el logaritmo será indefinido.
Estas son las restricciones que he definido. Utilicé la notación con un dict en lugar de un objeto NonlinearConstraints, ya que parece que este método acepta desigualdad estricta (np.exp(-x[0]**x[1]) debe ser estrictamente menor que 1). Tal vez he malinterpretado la documentación al respecto.
def constraints(x):
return [1 - np.exp(-(x[0]) ** x[1])]
Para maximizar la verosimilitud, minimizo la verosimilitud negativa.
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)]))
Según lo entiendo, los dos argumentos nunca deberían ser elegidos de tal manera que mi código falle. Sin embargo, el algoritmo intenta acercar theta muy cerca de su límite inferior y tau muy cerca de su límite superior para que el logaritmo se vuelva indefinido.
¿Qué hace que mi código falle?
davy.ai
Tu código falla porque los límites que has establecido en
theta
ytau
no son suficientes para asegurar que se cumplan las condiciones en la declaraciónif
. En particular, establecer el límite inferior entheta
ytau
en1e-15
no es suficiente para asegurar que estas variables sean estrictamente positivas. De manera similar, restringir el límite superior entheta
ytau
a10
no garantiza quetheta / d
sea estrictamente menor que1
, lo cual es necesario para que la segunda condición en la declaraciónif
sea falsa.Para solucionar esto, debes establecer límites inferiores en
theta
ytau
a un valor positivo pequeño (por ejemplo,1e-6
), y aumentar el límite superior entheta
ytau
a un valor mayor que asegure que se cumplan las condiciones en la declaraciónif
. También puedes considerar imponer restricciones adicionales entheta
,tau
y en el valor denp.exp(-(theta/d)**tau)
para asegurar que el algoritmo converge a una solución válida.