CODEWITHSUNDEEP
pythonpython-3.xscipyminimize
Jack.O.
Jack.O.

Reputation: 171

"scipy.optimize.minimize" How to force the coefficients to be not zero

I am trying to fit a function p which depends on two variables x,T. The data for the p,T,x are provided via an excel sheet with pandas. The following code works quite well.

import pandas as pd
import os
from scipy.optimize import minimize
import numpy as np

df = pd.read_excel(os.path.join(os.path.dirname(__file__), "./DataTest.xlsx"))
df = df.sort_values('x')

T = np.array(df['T'], dtype=float)
x = np.array(df['x'], dtype=float)
p = np.array(df['p'], dtype=float)
p0 = 67.17

def cav2(pars, T, x): # function p(T,x)
    a,b,c,d,e,f = pars
    return x * p0 + x * (1 - x) * (a + b * T + c * T ** 2 + d * x + e * x * T + f * x * T ** 2) * p0

def resid(pars, T, x):
    return ((p - cav2(pars, T, x)) ** 2).sum()

def constr(pars):
    return np.gradient(cav2(pars, T, x))

con1 = {'type': 'ineq', 'fun': constr}
pars0 = np.array([0.1, 0.1, 0.1, 0.1, 0.1, 0.1], dtype=float)
res = minimize(resid, pars0, args=(T, x), method='cobyla', options={'maxiter': 50000}, constraints=con1)

print("a = %f , b = %f, c = %f, d = %f, e = %f, f = %f" % (res.x[0], res.x[1], res.x[2], res.x[3], res.x[4], res.x[5]))

The last printgives me the coefficients of my function:

a = 2.891584 , b = 0.000000, c = -0.000000, d = 0.792256, e = -0.000000, f = 0.000000

That brings me to my actual problem. Because some of the coefficients become zero, it makes the function p(T,x) independent from T, which i do not want. To be clear, at the moment cav2(res.x, 300, 0.1) gives the same result as for example cav2(res.x, 500, 0.1).

Is there an (easy) way in scipy.optimize.minimize to force all of the coefficients to accept a value greater than zero ?

Thanks

Upvotes: 2

Views: 8651

Answers (1)

Jonathan Guymont
Jonathan Guymont

Reputation: 497

Some optimizer support bound constraint (e.g. L-BFGS-B) on coefficients.

import pandas as pd
import os
from scipy.optimize import minimize
import numpy as np

T = np.random.normal(10)
x = np.random.normal(10)

p0 = 67.17

# Fake true parameters
a, b, c, d, e, f = np.random.uniform(-1, 1, size=6)

# targets
p = x * p0 + x * (1 - x) * (a + b * T + c * T ** 2 + d * x + e * x * T + f * x * T ** 2) * p0


def cav2(pars, T, x): # function p(T,x)
    a, b, c, d, e, f = pars
    return x * p0 + x * (1 - x) * (a + b * T + c * T ** 2 + d * x + e * x * T + f * x * T ** 2) * p0


def resid(pars, T, x):
    return ((p - cav2(pars, T, x)) ** 2).sum()


def constr(pars):
    return np.gradient(cav2(pars, T, x))

# this will force all parameters to be positive
bounds = [(0, None), (0, None), (0, None), (0, None), (0, None), (0, None)]
pars0 = np.array([0.1, 0.1, 0.1, 0.1, 0.1, 0.1], dtype=float)

res = minimize(resid, pars0, args=(T, x), method='L-BFGS-B', options={'maxiter': 50000}, bounds=bounds)

print("a = %f , b = %f, c = %f, d = %f, e = %f, f = %f" % (res.x[0], res.x[1], res.x[2], res.x[3], res.x[4], res.x[5]))

The way the bounds works is (lower, upper) and putting None means no bound is applied. So if you don't want a bound on the first parameter, for example, you can replace bounds with:

[(None, None), (0, None), (0, None), (0, None), (0, None), (0, None)]

Upvotes: 1

Related Questions