Fitting two sets of equations and two data sets with NonlinearModelFit

Question

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?

Answer 1

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."

Answer 2

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

instead.)

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.