Combined fits

[fjosw]

@PiaLJP and I started discussing a possible idea for an interface for combined fits without breaking backward compatibility. We thought about extending fits.least_squares in such a way that inputting lists for x and y and a single function for func results in the known behavior. Alternatively one can pass (non-nested) dictionaries for x, y and func which have to have identical keys. Following this proposal one would define a separate fit function for every data set, but the fit functions can share parameters. The keys of the dictionaries are only used for identification and would in no way propagate into the fit result.

A simple example would be fitting two correlations functions with the same mass but different amplitudes:

def func_a(a, x):
    return a[1] * anp.exp(-a[0] * x)

def func_b(a, x):
    return a[2] * anp.exp(-a[0] * x)

func_dict = {"a": func_a,
             "b": func_b}

The check for the number of fit parameters could be done by performing the usual check on the sum of all fit functions. This should require very little modification of the existing code.

The biggest changes would be in the part where chisquare is constructed, the part for the error propagation probably only requires minimal changes.

A bit more care is required in the correlated case. We could think about the option of passing another dictionary containing the block diagonal covariances for each data set.

Any thoughts?

[s-kuberski]

Thanks for starting to think about it. I like the idea of the implementation but have to think about it and if there is a better way of doing things (from the perspective of the user). Why using dictionaries with (arbitrary?) keys, when lists of lists would suffice? I think we should spend some time to think about the smoothest way to integrate this feature (as it is possible already to do such fits, but not very easily). In the end, also for the chisquared etc., you just need N tuples (x, y, func), where x may be an array, y is your data and func the corresponding function for this pair of (x, y). In that sense, you could even have the same interface as before, but just add another list containing the fit functions (instead of a single function). However, this is probably also not the user-friendliest way (or is it?).

In the correlated case, you would most likely need a full covariance matrix, right? If you, e.g., fit 2- and 3-point functions from the same configurations, and you would like to do this in a correlated fashion, the "off-block-diagonal" terms should be taken into account. I would like to enable the user to provide a covariance matrix, regardless of the other changes. When scanning over fit ranges, the repeated computation of (smaller and smaller) covariance matrices with basically the same entries generates some overhead. In my opinion, it would be enough to ask for a full covariance matrix and the user may decide how to construct it. It could be very helpful to provide a function to construct block-diagonal covariance matrices from several smaller ones.
[This discussion could change a bit if we implement the new way of computing chi_exp.]

As there are some other open questions regarding improvements of the fit functions: Is this the right place or should we separate the different issues?

[fjosw]

We thought about using dictionaries as we had the impression that this gives the user a better idea which sets of data enter a fit. It would also make adding or removing a data set simpler (simply add or pop a key in all three dicts).
As we can simply loop over all keys one would also not have to be so careful with the order in which one initializes the three inputs (and have an explicit check here).
Another benefit for the implementation would be that there is a clear difference between the standard mode of operation (via list input) and the new mode (via dict input). I can imagine situations with multiple x or y values where a nested list can be somewhat ambiguous. We are, however, of course open to other suggestions!

I agree on the suggestion to add a way of passing a covariance. Concerning efficiency it might be beneficial to not have a block diagonal matrix with explicit zeros but use a more efficient way of computing the product, for example via numpy.einsum (maybe not the most pressing issue, though).

If your other suggestions are unrelated to these type of combined fits I would suggest to open a separate discussion.

[s-kuberski]

I see your point concerning the dictionaries. I had the impression that it is unnecessary to explicitly assign some keys to the data, but maybe it even contributes to more readable workflows. Yet another implementation could be done via an explicit fit-parameter class that contains the x, y, and the function. Here, one could overload the addition of objects to create a new, combined fit object. But this might also be an overkill.

Maybe a sample implementation would help to see whether your suggestion needs any improvements for practical purposes?

[PiaLJP]

I can get started on that sample implementation based on dictionaries (based on my current implementation which I use to fit the sigma terms) within the next few days.

[fjosw]

The OOP solution sounds like an interesting idea, we could also consider that (I see some issues with backward compatibility, though). Anyways, what is happening under the hood hardly depends on the interface so we could sort that out later. Thanks for getting started @PiaLJP, just open a PR once you have a first version and then we can discuss in more detail.

[PiaLJP]

In preparation for a PR I just uploaded a 'combined_fits.py' to my fork and also an example so in case you would like to look at that already... In the course of next week I will then include this in the existing function fits.least_squares

[fjosw]

Looks good already! One simple test could be to fit a constant to a data set in the traditional way and then split the data set in two and perform a combined fit to a constant. Both fits should give the same result. Testing all new exceptions related to invalid input would also be a good idea.