Scipy & Optimize: Minimize example, how to add constraints?

Question

I have a dataset, and I'd like to find a mixed gaussian model by least square error method.

The code is like this:

from sklearn.neighbors import KernelDensity
kde = KernelDensity().fit(sample)

def gaussian_2d(x,y,meanx,meany,sigx,sigy,rho):
    # rho <= 1
    part1 = 1/(2*np.pi*sigx*sigy*sqrt(1-0.5**2))
    part2 = -1/2*(1-rho**2)
    part3 = (((x-meanx)/sigx)**2-2*rho*(x-meanx)*(y-meany)/(sigx*sigy)+((y-meany)/sigy)**2)
    return part1*exp(part2*part3)

def square_error(f1,f2, u1,v1,sigu1,sigv1,rho1, u2,v2,sigu2,sigv2,rho2, u3,v3,sigu3,sigv3,rho3):
    # 1. Generate Mixed Gaussian Model
    def gaussian1(x,y):
        return gaussian_2d(x,y,u1,v1,sigu1,sigv1,rho1)
    def gaussian2(x,y):
        return gaussian_2d(x,y,u2,v2,sigu2,sigv2,rho2)
    def gaussian3(x,y):
        return gaussian_2d(x,y,u3,v3,sigu3,sigv3,rho3)
    mixed_model = f1*gaussian1(x,y)+f2*gaussian2(x,y)+(1-f1-f2)*gaussian3(x,y)
    # 2. Calculate the sum of square error
    sum_error = 0
    for point in sample:
        error = (exp(mixed_model(point)) - exp(kde.score(point)))**2
        sum_error += error
    return sum_error

# How can I add constraints:
# f1+f2 <= 1
# rho1,2,3 <= 1
result = sp.optimize.minimize(square_error)

But I don't know how to add constrains in minimize method. The example in http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html#scipy.optimize.minimize is very hard to understand.

UPDATE: This is what I end up with,

result = sp.optimize.minimize(
    square_error,
    x0 = [0.2,0.5,
          1,1,1,1,0.3,
          1,1,1,1,0.3,
          1,1,1,1,0.3,],
    bounds = [(0., 1.),(0., 1.),
              (None, None),(None, None),(0., None),(0., None),(0., 1.),
              (None, None),(None, None),(0., None),(0., None),(0., 1.),
              (None, None),(None, None),(0., None),(0., None),(0., 1.),])

But it gives me TypeError: square_error() takes exactly 17 arguments (1 given), what's the problem?

possible duplicate of scipy minimize with constraints

Daniel Lenz
– Daniel Lenz

2015-09-11 12:19:31 +00:00
Commented Sep 11, 2015 at 12:19 — Daniel Lenz
– Daniel Lenz, Commented Sep 11, 2015 at 12:19

Daniel Lenz · Accepted Answer · 2015-09-11 11:28:39Z

2

You can add only add bounds if the solver supports them, so only for method='L-BFGS-B', TNC and SLSQP.

The bounds are passed via a sequence of (min, max) tuples whose length corresponds to the number of your parameters. An example for fitting with 3 parameters would be:

result = sp.optimize.minimize(
                        square_error,
                        method='L-BFGS-B',
                        bounds=[(0., 5.), (None, 1.e4), (None, None)])

Here, None corresponds to no bound. I'm afraid that constraints on a combination of parameters such as f1+f2 <= 1 in your example is not possible within the framework of bounds in scipy.minimize.

You can, however, simply return np.inf in your cost function if your bounds are violated. I'm not sure about the stability of this, though:

def square_error(f1,f2, other_args):
    if f1+f2 <= 1:
        return np.inf
    # rest of the cost function

Moreover, I'd suggest to use the python implementation of multivariate Gaussians instead of creating them from scratch. That will speed up your fitting, help to avoid bugs and be more readable.

edited Sep 11, 2015 at 11:28

answered Sep 11, 2015 at 9:33

Daniel Lenz

3,9371 gold badge22 silver badges37 bronze badges

Sign up to request clarification or add additional context in comments.

7 Comments

cel Over a year ago

I'm afraid that constraints on a combination of parameters such as f1+f2 <= 1 in your example is not possible. This is not true. There are solvers especially for constrained optimization that support them.

Daniel Lenz Over a year ago

Could you provide a link? I was referring to the solvers that are included in scipy.minimize.

cel Over a year ago

Just have a look at the constraints keyword. OP is trying to define constraints rather than simple bounds.

Daniel Lenz Over a year ago

Thanks for pointing that out, I'm afraid I misunderstood the OP. I'll update my answer accordingly. Meanwhile, check out stackoverflow.com/questions/20075714/… for a great answer.

ZK Zhao Over a year ago

For the np.random.multivariate_normal, I need the exact probability density value, and the module seems only provide sample?

|

Collectives™ on Stack Overflow

Scipy & Optimize: Minimize example, how to add constraints?

1 Answer 1

7 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

7 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related