# A case of complicated fitting of multiple datasets

I've stumbled upon a difficulty while trying to do some fitting which may be called sorta complicated. I'll try to introduce my problem step by step.

To begin with, I have a number of datasets datasets. Each dataset consists of a number of points in 2D space. Those points look like {x,y}. Then, I have a set of functions funcs. For each dataset I need to pick a function that fits the best and find fitting parameters. The good thing is that I have some kind of generalized parameters pars which are the same for all functions (yet the functions are significantly inequivalent), so I can accomplish what I need with a code like this:

Table[
fits = Table[NonlinearModelFit[datasets[[i]],funcs[[j]],pars,x],{j,Length@funcs}];
(First@SortBy[#["EstimatedVariance"]&]@fits)["BestFitParameters"]
,{i,Length@datasets}]


Simple enough. The code above will give me best fitting parameters for each dataset with best function used for fitting. I can also easily get name (or number) of a function that provide best fitting for each dataset, if I want to.

Now, to what makes fitting complicated. One of pars, let it be the first one, is not independent. I still have to find the best value for it but this value must be the same for all datasets and all funcs. I tried to do the following:

NMinimize[
Total@Table[
fits = Table[NonlinearModelFit[datasets[[i]],funcs[[j]],pars[[2;;]],x],{j,Length@funcs}];
(First@SortBy[#["EstimatedVariance"]&]@fits)["EstimatedVariance"]
,{i,Length@datasets}]
,pars[[1]]]


I also tried to use FindMinimum instead of NMinimize. With this I expected to get the best value of pars[[1]] in order to fix that value and execute the first code with pars[[2;;]] instead of pars to get best functions and their fitting parameters for each dataset. But that didn't work. Inside of NMinimize or FindMinimum NonlinearModelFit sees pars[[1]] as a symbol instead of a number so it can't do the fitting. Moreover, most fitting algorithms are based on derivatives but calculating a derivative of NonlinearModelFit doesn't seems possible.

Just to be precise, pars are symbols but I expected NMinimize to assign some numerical values to pars[[1]] and try fitting with it, then try another value for it and somehow find the best value. But unfortunately pars[[1]] doesn't becomes a number inside of NMinimize and that ruins my plan.

So, any clue on how to get this to work?

• You need to fit a single model that covers all datasets. I'd first fit each dataset separately to get initial estimates of all of the parameters to go into the all-encompassing model and use the average value of the common parameter as its starting value. Despite there being potentially many parameters to estimate, if the starting values are close enough, it shouldn't take a lot of computing time. – JimB Apr 8 '16 at 20:49
• There isn't a built in function for this kind of model fitting. There is some independently written code for this however. Please take a look at this question here : community.wolfram.com/groups/-/m/t/135933 – Searke Apr 8 '16 at 21:10
• Upon reading this again, maybe things are more complex that I first thought. But for this "common" parameter, is it involved in the same way in all functions being considered? In other words, suppose the parameter is $a$. Is it always $a x$ ? Or might one function be $a + b x$ or $b x^a$ ? If the functions are "significantly inequivalent" as you state and each dataset might be best fit by a different function, how can you possibly state that there is a common parameter with the same true value for all datasets? Have I misunderstood something? If the functions were made available.... – JimB Apr 8 '16 at 22:12
• @JimBaldwin, the "common parameter" is involved in all functions in the same way. Basically, pars are something like {phase, shift, amplitude}. Conditions in my task have physical origin so I assure you the common parameter is common. This is as much as I can tell without clearly identifying my task. And about computing time – it's not a problem at all. – Helyrk Apr 9 '16 at 7:50
• @Searke, I will definitely have a closer look at the topic you suggested. – Helyrk Apr 9 '16 at 7:54

I'm assuming the objective is to find "the best" set of models (all with a single common parameter) with potentially different models for each of many data sets.

Using NonlinearModelFit one could fit a combined model with all data sets and a specific set of models. One would obtain the AIC value for each of the all possible specific sets of models (one model assigned to each data set). One could choose the set of models that results in the smallest AIC value.

But this assumes a common error variance for all data set/model combinations. That might be hard to justify.

One way to get each model with a separate estimate of error variance is to use the LogLikelihood function for each data set/model, add up the contributions from each data set/model, use FindMaximum to obtain the parameter estimates, and finally the associated AIC value. Rank the model combinations by the AIC value to choose "the best" model.

The could involve many, many overall models. If there are $d$ data sets and $m$ models, then there are $m^d$ overall models. When this resulting number is large, then one can easily be accused of data mining and trust that "the best" model is the best model is likely to be low. And one must remember that AIC ranks the models. It is only a relative assessment of fit and not an absolute assessment of fit. One might end up with the best of a bunch of bad models or (which is at least theoretically possible) the best of a bunch of great models.

Here is an outline of the approach with a specific set of just two potential models with two data sets.

Suppose we have the following: x1 and x2 are lists of the values of the predictor variables and y1 and y2 are the lists of observed values. The fixed effects of the two models are defined by

model1[x_, a_, b_, θ_] := a + b Sin[x + θ]
model2[x_, c_, d_, θ_] := c + d Cos[x^2 + θ]


where θ is the common parameter.

We assume that the square root of the error variances of the two models are σ1 and σ2, respectively.

If we wanted to determine the estimates for data set 1 with model 1 and data set 2 with model 2, we could write

logL1 = Simplify[LogLikelihood[NormalDistribution[0, σ1], y1 - model1[x1, a1, b1, θ]]];
logL2 = Simplify[LogLikelihood[NormalDistribution[0, σ2], y2 - model2[x2, c2, d2, θ]]];
sol12 = FindMaximum[{logL1 + logL2, σ1 > 0 && σ2 > 0},
{{a1, 0}, {b1, 1}, {c2, 0.5}, {d2, 0.06}, {θ, 1}, {σ1, 0.1}, {σ2, 2}}]

aic12 = -2 sol12[[1]] + 2*7


This would be repeated for all model assignments to data sets to be considered. (And note that one would probably want to fit each data set/model combination separately at first to obtain starting values for the LogLikelihood approach.)

Update

Just a short example to show that NonlinearModelFit and the use of LogLikelihood gives one the same estimates...

(* Define model *)
model[x_, a_, b_, θ_] := a + b Sin[x + θ]

(* Generate some samples *)
n = 100;
x = Table[2 π i/n, {i, n}];
y = model[x, 0, 1, 0.5] + RandomVariate[NormalDistribution[0, 0.1], n];

(* NonlinearModelFit estimates *)
nlm = NonlinearModelFit[Transpose[{x, y}], model[z, a, b, θ], {a, b, θ}, z];

nlm["BestFitParameters"]
(* {a->-0.00252912,b->1.00503,θ->0.50748} *)
nlm["EstimatedVariance"]^0.5
(* 0.114323 *)

(* LogLikelihood function estimates *)
logL = Simplify[LogLikelihood[NormalDistribution[0, σ], y - model[x, a, b, θ]]];
sol = FindMaximum[{logL, σ > 0}, {{a, 0}, {b, 1}, {θ, 0.5}, {σ, 0.1}}]
(* {76.5017,{a->-0.00252912,b->1.00503,θ->0.507478,σ->0.112595}} *)
(* Adjust estimate of the square root of the variance to match what
NonlinearModelFit produces. 3 is subtracted from n below because there
are 3 fixed effects parameters being estimated *)
(σ /. sol[[2]]) Sqrt[n/(n - 3)]
(* 0.114323 *)

• Wait, I didn't get something in your approach. You have built a function that gives fitting error for two models and two datasets logL1 + logL2 (or it gives some kind of "goodness"). This is good. But each logL includes only one model, right? So you have chosen models for datasets manually. Can I possibly choose models automatically? – Helyrk Apr 11 '16 at 12:10
• If you want the best fit from the models you have, then, yes, that can be automated. If you want an adequate fit, then, no, that can't be automated. (At minimum one has to look at residual diagnostics for that.) I didn't include an automatic process because I didn't want to bury the main issue of simultaneously getting a fit for a selection of models with a single "relative" goodness-of-fit criterion. It's still not clear to me if you want to first screen for the "best" model for each dataset and then fit all together or do all possible combinations. How many datasets and how many models? – JimB Apr 11 '16 at 14:54
• The main purpose of my fitting is to figure out which model corresponds to each dataset. To classify datasets by models. The connection between datasets that makes phase a common parameter is a hint for me. By the way, I have 20 datasets and 5 models. If possible, I want to estimate the value of phase, but the main goal is to get datasets classification. And one more thing. Because the common parameter is phase I was able to check its values with small step and get some result which should be close to "best". – Helyrk Apr 11 '16 at 17:11
• Thanks. Got it. I'll write something up when I get home tonight. Just note that the best model for a data set might not be the best model when the common parameter is forced to match all data sets. If the variability of the estimated values of phase is low, then this likely won't be a problem. – JimB Apr 11 '16 at 17:16

Well, I have finally found a satisfactory solution. In such a case I can shift my datasets along X axis to make their dataranges non-overlapping and then use Piecewise to apply any set of models with any number of common parameters. Basically, It would look somewhat like this:

NonlinearModelFit[Join[datasets], Piecewise[Table[{selectedFuncs[[i]][x], Min[datasets[[i, All, 1]]] <= x <= Max[datasets[[i, All, 1]]]}, {i, Length@datasets}], pars, x]


This would require to check all possible combinations of funcs and do the previous code with selectedFuncs as the combination.