Analytical Constrained Fitting

[tBuLi]

I am currently working on analytic constrained fitting using generalised Lagrange Multipliers. This method finds the saddle points of the model and therefore allows for analytical expressions for the parameters. These can then be used to calculate the value from the data.

Since multiple solutions are possible, we want to enable the user to select the desired solution as sexy as possible. To give an example of how this would work, I reproduce an example I found here: http://www.asp.ucar.edu/colloquium/1992/notes/part1/node36.html.

The challenge is to calculate the angles of a triangle from a set of measurements on said angles. The constraint is that we know the angles should add up to 180 degrees.

This is how this will look in symfit after vector valued (#15) models have been introduced.

a, b, c = parameters('a, b, c')  # fitted angles
a_i, b_i, c_i = variables('a_i, b_i, c_i')  # variables representing the i'th measured angles.

model = {a_i: a, b_i: b, c_i: c}
constraints = [Eq(a + b + c, 180)]

analytic_fit = LagrangeMultipliers(model, constraints=constraints)
extrema = analytic_fit.extrema
>>> [Extremum(a=2*a_i/3 - b_i/3 - c_i/3 + 60, b=-a_i/3 + 2*b_i/3 - c_i/3 + 60, c=-a_i/3 - b_i/3 + 2*c_i/3 + 60)]

Extremum is a NamedTuple containing an analytical extremum to the system of equations.The beauty of this is that the user can now provide this solution to his favourite fitting object as the one to use to calculate the parameters by averring over all the data:

data = np.array([
    [10., 10., 10., 10., 10., 10., 10.],
    [100., 100., 100., 100., 100., 100., 100.],
    [70., 70., 70., 70., 70., 70., 70.],
])

fit = Fit(model, a_i=data[0], b_i=data[0], c_i=data[0], extremum=extrema[0])
fit.execute()

I also plan on building a ConstrainedFit object to automate the entire process, if it is possible to determine the desired solution from the bounds provided on parameters. In case of only one solution this is easy of course. In case no solution can be found or multiple, the object should (?) raise an Exception.

data = np.array([
    [10., 10., 10., 10., 10., 10., 10.],
    [100., 100., 100., 100., 100., 100., 100.],
    [70., 70., 70., 70., 70., 70., 70.],
])

fit = ConstrainedFit(model, a_i=data[0], b_i=data[0], c_i=data[0], constraints=constraints)
fit.execute()