Scipy minimize / Scipy Curve fit / lmfit

Question

log(VA) = gamma - (1/eta)log[alphaL^(-eta) + betaK^(-eta)]

I'm trying to estimate the above function with nonlinear least squares. I am using 3 different packages (Scipy-minimize, Scipy-curve_fit and lmfit - Model) for this but I find different parameter results in each one. I can't understand why. I would be very grateful if anyone can help with a solution or offer a different solution method.

SCIPY-MINIMIZE

import numpy as np
from scipy.optimize import minimize, curve_fit
from lmfit import Model, Parameters

L  = np.array([0.299, 0.295, 0.290, 0.284, 0.279, 0.273, 0.268, 0.262, 0.256, 0.250])
K  = np.array([2.954, 3.056, 3.119, 3.163, 3.215, 3.274, 3.351, 3.410, 3.446, 3.416])
VA = np.array([0.919, 0.727, 0.928, 0.629, 0.656, 0.854, 0.955, 0.981, 0.908, 0.794])

def f(param):
    gamma = param[0]
    alpha = param[1]
    beta  = param[2]
    eta   = param[3]
    VA_est = gamma - (1/eta)*np.log(alpha*L**-eta + beta*K**-eta)
    
    return np.sum((np.log(VA) - VA_est)**2)

bnds = [(1, np.inf), (0,1),(0,1),(-1, np.inf)]
x0 = (1,0.01,0.98, 1)
con = {"type":"eq", "fun":c}

result = minimize(f, x0, bounds = bnds)

print(result.fun)
print(result.message)
print(result.x[0],result.x[1],result.x[2],result.x[3])

SCIPY-MINIMIZE - OUT

0.30666062040617503
CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL
1.0 0.5587147011643757 0.9371430857380681 5.873041615873815

SCIPY-CURVE_FIT

def f(X, gamma, alpha, beta, eta):
    L,K = X
  
    return gamma - (1/eta) * np.log(alpha*L**-eta + beta*K**-eta)

p0 = 1,0.01,0.98, 1

res, cov = curve_fit(f, (L, K), np.log(VA), p0,  bounds = ((1,0,0,-1),(np.inf,1,1,np.inf)) )
gamma, alpha, beta, eta = res[0],res[1],res[2],res[3] 
gamma, alpha, beta, eta

SCIPY-CURVE_FIT - OUT

(1.000000000062141,
 0.26366547263939205,
 0.9804436474926481,
 13.449747863921704)

LMFIT-MODEL

def f(x, gamma, alpha, beta, eta):
    L = x[0]
    K = x[1]
    
    return gamma - (1/eta)*np.log(alpha*L**-eta + beta*K**-eta)

fmodel = Model(f)
params = Parameters()
params.add('gamma', value = 1,    vary=True, min = 1)
params.add('alpha', value = 0.01, vary=True, max = 1, min = 0)
params.add('beta',  value = 0.98, vary=True, max = 1, min = 0)
params.add('eta',   value = 1,    vary=True, min = -1)

result = fmodel.fit(np.log(VA), params, x=(L,K))
print(result.fit_report())

LMFIT-MODEL - OUT

[[Model]]
    Model(f)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 103
    # data points      = 10
    # variables        = 4
    chi-square         = 0.31749840
    reduced chi-square = 0.05291640
    Akaike info crit   = -26.4986758
    Bayesian info crit = -25.2883354
##  Warning: uncertainties could not be estimated:
    gamma:  at initial value
    gamma:  at boundary
    alpha:  at boundary
[[Variables]]
    gamma:  1.00000000 (init = 1)
    alpha:  1.3245e-13 (init = 0.01)
    beta:   0.20130064 (init = 0.98)
    eta:    447.960413 (init = 1)

score 3 · Answer 1 · answered Sep 04 '21 at 05:53

A fitting algorithm always seeks for a local minimizer of the underlying least-squares problem. Note that your problem is convex but not strictly convex. Consequently, there's no unique global minimizer. But each local minimizer is a global one. By evaluating the first function f for each found solution, we can observe that they all have the same objective function value. Hence, each solution is a global minimizer.

Why does each method find a different minimizer? The reason is simple. Each one uses a different algorithm to solve the underlying nonlinear optimization problem, e.g. scipy.optimize.minimize uses the 'L-BFGS-B' algorithm while scipy.optimize.curve_fit uses scipy.optimize.least_squares with the Trust Region Reflective algorithm ('TRF'). In short, you can only expect the same solution for different algorithms for a strictly convex problem.

Scipy minimize / Scipy Curve fit / lmfit

1 Answers1