2

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

Jack.O.
  • 171
  • 1
  • 14
  • Do you have a range of values you wish to bound `T` or `x` over? If so, specifying this may produce a local minimum with non-zero coefficients. – iacob Feb 25 '19 at 14:25
  • Actualy I don't realy understand how bounds work. If I have something like this `bnds = ((1, 0.1, 1),(300, 1, 67))` and `bounds=bnds` in the `minimize` line. Does this stand for `bnds = ((Tmin, xmin, pmin),(Tmax, xmax, pmax))` ? If so it unfortunately does not solve my problem. – Jack.O. Feb 25 '19 at 14:52

1 Answers1

1

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)]
Jonathan Guymont
  • 497
  • 2
  • 11
  • Thanks! Your solution with `method='L-BFGS-B'` works in way but unfortunately the method does not support constraints, which I need for my condition where the derivation of the function is positive everywhere (see `def constr()`). On the other hand the 'method='cobyla'` does not support bounds. =( – Jack.O. Feb 25 '19 at 16:00
  • But you gave me a good hint to look how to formulate bounds as constraints https://stackoverflow.com/questions/12781622/does-scipys-minimize-function-with-method-cobyla-accept-bounds So thanks for that! – Jack.O. Feb 25 '19 at 16:04