# Fitting two sets of equations and two data sets with NonlinearModelFit

I have 2 nonlinear equations of the type:

model_shi_real=(((β - γ)/2)*(Cos[π*α/2]/(Cosh[(1 - α)*Log[x*τ]] + Sin[α*π/2])))


and

model_shi_imag=(γ + ((β - γ)/2)*(1 - Sinh[(1 - α)*Log[x*τ]]/(Cosh[(1 - α)*Log[x*τ]] + Sin[α*π/2])))


with the real value parameters α, β, γ, and τ for both equations.

I further have 2 data sets:

data_shi_real={{10, 2.08228}, {12.6899, 2.04886}, {16.103, 2.00934}, {20.4341, 1.96466}, {25.9292, 1.93907}}


and

 data_shi_imag={{10, 0.37359}, {12.6899, 0.34438}, {16.103, 0.327364}, {20.4341, 0.307265}, {25.9292, 0.29121}}.


I have restrictions to the parameters as

1<α<0, β>γ and τ>0.


I can use NonlinearModelFit in order to fit each data[i] set to its function and get good fits. But all the parameters should be the same for both fits.

I used:

real = NonlinearModelFit[Data_shi_real, {model_shi_real[\[Alpha], \[Beta], \[Gamma], \[Tau],x], {1 > \[Alpha] > 0, \[Beta] > 5, \[Tau] > 0, 2 > \[Gamma] > 1}}, {\[Alpha], \[Beta], \[Gamma], \[Tau]}, x, MaxIterations -> 5000]


and

imag = NonlinearModelFit[data_shi_imag, {model_shi_imag[\[Alpha], \[Beta], \[Tau], \[Gamma], x], {1 > \[Alpha] > 0, \[Beta] > 5, \[Tau] > 0, 2 > \[Gamma] > 1}}, {\[Alpha], \[Beta], \[Tau], \[Gamma]}, x, MaxIterations -> 5000]


Obtaining the following when fitting the real part and using to simulate the imag part:

{\[Alpha] -> 0.801163, \[Beta] -> 10.9567, \[Gamma] -> 1.15813, \[Tau] -> 7046.43}


Obtaining the following when fitting the imag part and using to simulate the real part:

{\[Alpha] -> 0.722099, \[Beta] -> 24.1831, \[Tau] -> 8529.64, \[Gamma] -> 2.


How can I create a simple code with Mathematica, so that I can do a simultaneous fit resulting in the parameters to be the same?

• Welcome to Mathematica Stack Exchange. If you haven't done so already, please take the tour at mathematica.stackexchange.com/help. Also, if you post a minimally sufficient example (including data or an appropriately-sized subset of the data), you'll quicker and better help.
– JimB
Mar 24, 2017 at 3:10
• You also might want to look at mathematica.stackexchange.com/questions/866/… and mathematica.stackexchange.com/questions/138174/…. But don't forget you need to decide on whether the error variances should be the same or different between the two models sharing some common parameters.
– JimB
Mar 24, 2017 at 4:32
• @JimBaldwin, Thank you for respond. An example of data is data1(shi_real)= {{10, 2.08228}, {12.6899, 2.04886}, {16.103, 2.00934}, {20.4341, 1.96466}, {25.9292, 1.93907}} and data2(shi_imag)= {{10, 0.37359}, {12.6899, 0.34438}, {16.103, 0.327364}, {20.4341, 0.307265}, {25.9292, 0.29121}}. I have some restrictions to the parameters as 1>α>0, τ>0 and β>γ. I didn't think about error variances and I don't know what´s the best option in this case, sorry. Mar 24, 2017 at 19:58
• Please edit your question and include the data in the question. That will also increase the chances of getting help (as some might not notice the data buried in a comment). And while it probably seems that I can't be satisfied, 10 data points for 6 parameters ($\alpha$, $\beta$, $\gamma$, $\tau$, $\sigma^2_{real}$, and $\sigma^2_{imag}$), is pretty minimal. I hope your real data set has many more points.
– JimB
Mar 24, 2017 at 20:00
• @JimBaldwin, Sorry to correct, but I have four parameters. In this case, I only have this data set. Mar 24, 2017 at 20:24

If one can assume that the error variances are identical, then one can use NonlinearModelFit. (This means one is "able to assume" rather than "willing to assume".) If one can't assume identical error variances, then consider using the LogLikelihood function described at this LogLikelihood example.

Assuming that the error variances are identical, one can use the following:

modelReal = (((β - γ)/2)*(Cos[π*α/2]/(Cosh[(1 - α)*Log[x*τ]] + Sin[α*π/2])));
modelImag = (γ + ((β - γ)/2)*(1 - Sinh[(1 - α)*Log[x*τ]]/(Cosh[(1 - α)*Log[x*τ]] + Sin[α*π/2])));
dataReal = {{10, 2.08228}, {12.6899, 2.04886}, {16.103, 2.00934}, {20.4341, 1.96466}, {25.9292, 1.93907}}
dataImag = {{10, 0.37359}, {12.6899, 0.34438}, {16.103, 0.327364}, {20.4341, 0.307265}, {25.9292, 0.29121}}
data = Flatten[{{1, #[[1]], #[[2]]} & /@ dataReal, {0, #[[1]], #[[2]]} & /@ dataImag}, 1]

fit = NonlinearModelFit[data,
{a modelReal + (1 - a) modelImag,
0 < α < 1 && β > γ > 0 && τ > 0},
{{α, 0.5}, {β, 12}, {γ, 10}, {τ, 7046.43}}, {a, x}];
fit["BestFitParameters"]
(* {α -> 0.136504850521697,β -> 12808872012093087,γ -> 2.3584229991919354*^-7,τ -> 25602.266522044003} *)


This approach creates a dummy variable to identify the two models. However, the fit is poor and the model is overparameterized. Here is the fit:

And here is the estimated parameter correlation matrix:

We see that the estimators of β and τ are perfectly correlated (implying overparameterization).

In this case having more data points between the existing data points is unlikely to fix the issue of near perfect correlation among the parameter estimators. But having more data between x = 0 and x = 10 might help.

Update

To attempt to decide what the problem actually is (bad theory, bad data, not-so-hot fitting algorithm, not enough data, not a wide enough range of predictor variables, etc.) one can set the parameters to known values, select some predictor variables, and try some fits. Here's the code that sets the parameters to what was found for the real data and uses the original predictor values. I've also added in a bit more noise than what appears to be the case.

parameters = {α -> 0.801163, β -> 10.9567, γ -> 1.15813, τ -> 7046.43, σ -> 0.05}
modelReal = (((β - γ)/2)*(Cos[π*α/2]/(Cosh[(1 - α)*Log[x*τ]] +       Sin[α*π/2])))
modelImag = (γ + ((β - γ)/2)*(1 - Sinh[(1 - α)*Log[x*τ]]/(Cosh[(1 - α)*Log[x*τ]] + Sin[α*π/2])))

(* Generate some data *)
SeedRandom[12345];
dataReal = Table[{x, modelReal /. parameters}, {x, {10, 12.6899, 16.103, 20.4341, 25.9292}}] +
RandomVariate[NormalDistribution[0, σ /. parameters], 5]
dataImag = Table[{x, modelImag /. parameters}, {x, {10, 12.6899, 16.103, 20.4341, 25.9292}}] +
RandomVariate[NormalDistribution[0, σ /. parameters], 5]

(* Combine into a single dataset *)
data = Flatten[{{1, #[[1]], #[[2]]} & /@ dataReal, {0, #[[1]], #[[2]]} & /@ dataImag}, 1];

(* Estimate parameters *)
fit = NonlinearModelFit[
data, {a modelReal + (1 - a) modelImag,
0 < α < 1 && β > γ > 0 && τ > 0},
{{α, 0.8}, {β, 11}, {γ, 1.2}, {τ, 7000}}, {a, x}];
fit["BestFitParameters"]
(* {α -> 0.19081295486508426, β -> 2.5544925230021898, γ -> 1.7927656092570712, t -> 0.14173596480881764} *)

Show[ListPlot[{dataReal, dataImag}],
Plot[{fit[1, x], fit[0, x]}, {x, 0, 26},
PlotLegends -> {"Real", "Imaginary"}]]


We see that (1) the predictions are fine for both curves and (2) the estimated coefficients are wildly different from the parameters generating the data. The parameter estimation correlation matrix still contains most entries very near -1 and +1. (Note that my choice of parameter values seems to have reversed the values for the real and imaginary data. Maybe I've slipped up somewhere.)

My conclusion would be that the two models with common parameters do not fit the data (but I don't have the subject matter background to declare if the problem is with the theory or the way the data was collected) and there is just not enough data nor are the predictor values widespread enough to get good estimates of the parameters. The fitting algorithm seems to work just fine.

Wood Allen's line about a restaurant comes to mind: "The food was bad and the portions too small."

• Thanks a lot for your help! The first code is functioning very well! Yes, the problem is in the model to describe the data. I was trying to reproduce only part of my results, but this is not good. I will modify the model and increase all the data to try to reproduce my results if you like I can show you my results then? But I have a question, why you write this: "σ /. parameters" when you generate Data, what are you saying in this case? Finally, I need to say thank you, again!!!!! Mar 29, 2017 at 1:04
• /. is ReplaceAll. It allows me to set a value for (in this instance) $\sigma$ just in this one part of the command and not affect the value of $\sigma$ elsewhere.
– JimB
Mar 29, 2017 at 2:32

I just add that fitting under the assumptions applied in Jim Baldwin's answer was a large part of the motivation for me to write the TransformedFit package (but please get the updated code from GitHub). As Jim mentions, these assumptions are far from generally applicable; but sometimes they are -- and if not, but one just wishes to see whether the model is reasonable or check if the fit converges, it's usually good enough even when the answer will only be a starting point for further analysis.

Now, this probably wouldn't be worth mentioning, except that by combining the real and imaginary parts of the model in one expression, we get a considerable simplification:

complexModel = realModel + I imagModel // FullSimplify
(* -> I β + (x (β - γ) τ)/(I x τ + E^((I/2) Pi α) (x τ)^α) *)


(or, if we wish to avoid using assumptions that are only correct for generic values of the variables, we can use

complexModel = realModel + I imagModel // TrigToExp // ExpandAll // Simplify


Forming the complex-valued data:

realData = {{10, 2.08228}, {12.6899, 2.04886}, {16.103, 2.00934},
{20.4341, 1.96466}, {25.9292, 1.93907}};
imagData = {{10, 0.37359}, {12.6899, 0.34438}, {16.103, 0.327364},
{20.4341, 0.307265}, {25.9292, 0.29121}};
complexData = realData;
complexData[[All, 2]] += I imagData[[All, 2]];


Now we can write:

Needs["TransformedFit"]
fit = ComplexFit[complexData, {
complexModel, {
0 < TransformedParameter[Re, α] < 1,
0 < TransformedParameter[Re, γ] < TransformedParameter[Re, β],
0 < TransformedParameter[Re, τ]
}
}, {{α, 0.5}, {β, 12}, {γ, 10}, τ}, x,
"CoordinateSystem" -> "Real" (* this means that parameter values are strictly real *),
"FitFunction" -> NonlinearModelFit
]


and get just the same answer that Jim got. In fact, if you add the option Hold -> True, you'll see what's passed to NonlinearModelFit, which will show that what the package does is just the same as what Jim did.

Obviously this doesn't really help as it is, but perhaps it'll be easier for you to see/fix the problem with your model if you don't have to worry about the fact that FindFit et al. usually don't like complex models/data unless split up in this way. At least you don't have to do the splitting manually.

• +1 Very nice! Does your function provide any estimates of precision? Parameter estimator covariance matrix? I ask because at least a few of the questions similar to the above involve what ends up to be essentially a linear relationship in the range of the observed predictor values. That makes getting separate estimates of the parameters difficult if not impossible. So while I constantly complain of too few data points, the issue for this question is more about getting a better set of predictor values.
– JimB
Mar 25, 2017 at 16:11
• @JimBaldwin it's simply a thin wrapper for either FindFit or NonlinearModelFit, so you can do fit["CorrelationMatrix"] as usual. The main potential difficulty is that the parameters are split into their real/imaginary, modulus/argument (etc.) components and these are treated as separate (possibly correlated) values. There is also the problem that by fitting effectively twice as many data points (one real and one imaginary for each actual point) then the error variance will not be estimated correctly by default. But as a starting point for a log-likelihood maximization, it can be helpful. Mar 25, 2017 at 16:22
• @OleksandrR. Thank you for your response! I was trying to use your answer but I have a problem, when I write this: "complexData = realData; complexData[[All, 2]] += I imagData[[All, 2]];" the program said: "Symbol called with 0 arguments; 1 argument is expected". And I was trying to find the reason for my own, but I could not find the answer. And sorry to bother you, but how can I plot the real part vs the imag part? Well, I want to say thank you for your help! Mar 29, 2017 at 1:20
• @Charlie your variable names contain underscores, which are not valid in symbol names in Mathematica. So I just renamed them, data_shi_real to realData and so on. Fix this and you should get the same results I did. Plotting would be a matter of ListPlot[Transpose[{Re[complexData], Im[complexData]}]]`. Hope that helps Mar 31, 2017 at 20:12