Constrained fit

[fjosw]

Mathias made me take a closer look at constrained fits and I believe the constrained fits triggered by the keyword argument const_par do not work as expected. I attached a small example in which one can achieve the constraint also by subtracting a constant from the data. The errors for the unconstrained parameter differ by a relevant factor between the two methods. I have the assumption that the two should agree in the limit in which the error of the constraint is small compared to the error on the data.

import numpy as np
import pyerrors as pe

dim = 10
x = np.arange(dim)
y = -0.06 * x + 5 + np.random.normal(0.0, 0.3, dim)
yerr = [0.3] * dim

oy = []
for i, item in enumerate(x):
    oy.append(pe.pseudo_Obs(y[i], yerr[i], 'test'))

def func(a, x):
    y = a[0] * x + a[1]
    return y

my_constant = pe.pseudo_Obs(5, 0.0001, 'test')

# Fit with constrained parameter
out = pe.least_squares(x, oy, func, const_par=[my_constant])

def alt_func(a, x):
    y = a[0] * x
    return y

alt_y = np.array(oy) - my_constant
[o.gamma_method() for o in alt_y]

# Fit with the constant subtracted from the data
alt_out = pe.least_squares(x, alt_y, alt_func)

print(out, alt_out)
print(out[0] == alt_out[0])

[s-kuberski]

Hm. You are right, this looks suspicious. I have made the following observations:

• If you perform a two-parameter fit and use the result for one of the resulting parameters to constrain a new one-parameter fits, the error on the free parameter does not change. I thought, that this was supposed to happen, since the fit parameters are not independent of each other.
• Whether you add/remove an offset using an Obs (see my code) exactly or constrain the fit should not make any difference for the mean value of the free parameter (this is the case). However, somehow one could expect the error to be independent of the procedure, as well (this is certainly not the case).
• The dependence of the error of the free fit parameter on the error of the constrained one is not entirely clear to me. My current implementation seems to lead to a result, where the error of the free parameter is insensitive to the constraint. This makes sense, when I think about how the error propagation is implemented. You are right: If the constraint has a very small error, the constrained-two-parameter fit should lead to a similar result as the one-parameter fit to the shifted data set.

I have to find the conceptual error in the current implementation...

dim = 10
x = np.arange(dim)
y = -0.06 * x +  np.random.normal(0.0, 0.25, dim)
yerr = [0.3] * dim

oy = []
for i, item in enumerate(x):
    oy.append(pe.pseudo_Obs(y[i], yerr[i], 'test'))

shift = pe.pseudo_Obs(1, .0000000001, 'shift')

oy = [o + shift for o in oy]
[o.gamma_method() for o in oy]

def func(a, x):
    y = a[0] * x + a[1]
    return y

# Fit with constrained parameter
out = pe.least_squares(x, oy, func, const_par=[shift])

def alt_func(a, x):
    y = a[0] * x
    return y

alt_y = np.array(oy) - shift
[o.gamma_method() for o in alt_y]

# Fit with the constant subtracted from the data
alt_out = pe.least_squares(x, alt_y, alt_func)

# Fit to the data with two free parameters
twop_out = pe.least_squares(x, oy, func)

print(out, '\n', alt_out, '\n', twop_out)

[s-kuberski]

Maybe the implementation is too simple and one should use Lagrange multipliers (together with MINUIT). There seems to be some discussion (and a fit package https://www.desy.de/~sschmitt/blobel/wwwcondl.html ) in the experimental community.

[fjosw]

As the issue seems to be rather complicated I would propose to remove the feature from the develop version in order to not delay the 2.0 release any further and put it on our agenda for the next minor release 2.1 for which I wanted to rework parts of the fitting routines anyway.

[s-kuberski]

This seems to be the best way to handle this right now.

[fjosw]

I temporarily removed the feature from the develop branch in commit 2a925ba.