Inequality constraints

The usual pattern for constraints within bumps is to set the value for one parameter to be some function of the other parameters. This does not allow contraints of the form \(a < b\) for parameters \(a\) and paramter \(b\).

Instead, along with the fit problem definition, you can supply your own penalty constraints function which adds an artificial value to the probability function for points outside the feasible region. The ideal constraints function will incorporate the distance from the boundary of the feasible region so that if the fitter is started outside forces the fit back into the feasible region.

The soft_limit value can be used in conjunction with the penalty to avoid evaluating the function outside the feasible region. For example, the function \(\log(a-b)\) is only defined for \(a > b\), so setting a constraint such as \(10^6 + (a-b)^2\) for \(a <= b\) and \(0\) along with a soft limit of \(10^6\) will keep the function defined everywhere. With the penalty value sufficiently large, the probability of any evaluation in the infeasible region will be neglible, and will not skew the posterior distribution statistics.

Define the model as usual

from bumps.names import *

def line(x, m, b):
    return m*x + b

x = [1, 2, 3, 4, 5, 6]
y = [2.1, 4.0, 6.3, 8.03, 9.6, 11.9]
dy = [0.05, 0.05, 0.2, 0.05, 0.2, 0.2]

M = Curve(line, x, y, dy, m=2, b=0)
M.m.range(0, 4)
M.b.range(0, 5)

Define the constraints as a function which takes no parameters and returns a floating point value. Note the value 1e6 in the penalty condition: this is the soft limit value which we will use to avoid evaluating the curve in the infeasible region.

def constraints():
    m, b = M.m.value, M.b.value
    return 0 if m < b else 1e6 + (m-b)**6

Attach the constraints to the problem. Give the soft limit value that is used for the constraints. Without the soft limit, the fit would stall since we started it at a deep local minimum near the true solution without constraints.

problem = FitProblem(M, constraints=constraints, soft_limit=1e6)

The constraint relies on the ability for python to access the parameters from the module. Furthermore, the parameters still “boxed”, and so you need to reference the value attribute to get the parameter value at the time the constraint is evaluated. Not an elegant solution, but it works. Eventually we will add constraint expressions such as M.m < M.b or M.m + M.b < 10 using the same infrastructure as equality constraints.

Download: inequality.py.