4
$\begingroup$

This question already has an answer here:

dataTC = {{2.2993`,49.`},{2.3106`,47.553`},{2.3202`,48.`},{2.3233`,48.`},{2.3233`,47.4`},{2.3353`,46.2`},{2.353`,47.49`},{2.3534`,48.3`},{2.3548`,47.`},{2.3805`,47.455`},{2.3903`,47.4`},{2.3917`,46.8`},{2.4093`,47.357`},{2.4135`,47.409`},{2.4268`,47.`},{2.4319`,47.5`},{2.4396`,49.4`},{2.4477`,45.3`},{2.4578`,47.224`},{2.5039`,46.985`},{2.5135`,47.2`},{2.5273`,46.669`},{2.5748`,46.13`},{2.5853`,45.827`},{2.5921`,45.1`},{2.6332`,45.533`},{2.6625`,45.331`},{2.6705`,45.174`},{2.6971`,41.4`},{2.7085`,45.008`},{2.721`,44.928`},{2.7317`,45.54`},{2.7539`,44.651`},{2.7578`,44.5`},{2.7656`,44.466`},{2.785`,44.401`},{2.803`,44.188`},{2.8097`,44.156`},{2.8132`,44.114`},{2.8539`,47.1`},{2.8561`,43.61`},{2.8643`,43.669`},{2.8983`,41.6`},{2.9084`,43.138`},{2.9401`,42.978`},{2.9505`,43.2`},{2.9782`,42.1`},{2.9967`,42.68`},{3.05`,42.316`},{3.0774`,41.6`},{3.0774`,43.`},{3.0884`,42.136`},{3.1217`,39.7`},{3.1363`,42.2`},{3.1552`,41.765`},{3.2257`,42.1`},{3.2376`,41.457`},{3.3026`,41.377`},{3.3533`,41.165`},{3.4227`,41.171`},{3.5025`,41.6`},{3.5041`,40.878`},{3.582`,40.848`},{3.5836`,41.6`},{3.6183`,38.3`},{3.6272`,40.6`},{3.8748`,40.6`},{4.051`,41.6`},{4.0704`,40.075`},{4.0739`,40.`},{4.108`,40.`},{4.1307`,40.1`},{4.5191`,39.4`},{4.5397`,39.9`},{4.5397`,40.2`},{4.5397`,39.6`},{4.5417`,41.1`},{4.5622`,40.`},{4.9344`,39.6`},{4.9344`,39.4`},{5.0096`,39.`},{5.3001`,39.1`},{5.4739`,39.29`},{5.6091`,38.7`},{5.6424`,38.7`},{5.9179`,39.7`},{5.9653`,38.7`},{6.1707`,38.9`},{6.1813`,39.7`},{6.2416`,38.6`},{6.2716`,38.4`},{6.2716`,38.85`},{6.2716`,39.06`},{6.4774`,39.4`},{6.5637`,38.3`},{6.705`,39.39`},{6.7745`,39.7`},{6.8434`,38.9`},{6.8707`,38.7`},{6.9116`,39.3`},{6.9791`,38.8`},{7.1673`,38.8`},{7.4218`,39.9`},{7.6213`,38.55`},{7.6213`,38.59`},{7.8755`,38.4`},{8.2136`,38.46`},{8.2136`,38.49`},{8.766`,38.85`},{8.766`,38.5`},{9.0296`,38.42`},{9.2857`,38.45`},{9.7777`,38.2`},{9.7777`,38.14`},{9.7777`,38.42`},{9.7777`,37.74`},{9.7777`,38.46`},{9.9866`,37.87`},{10.246`,38.43`},{10.694`,38.44`},{11.124`,38.65`},{11.457`,36.68`},{11.538`,38.28`},{11.779`,37.6`},{13.626`,37.6`},{13.763`,38.7`},{13.763`,38.46`},{13.763`,38.39`},{13.899`,38.9`},{15.065`,38.58`},{15.189`,37.`},{16.662`,37.6`},{16.662`,38.47`},{16.83`,38.69`},{16.83`,38.62`},{17.91`,38.83`},{18.015`,37.4`},{18.17`,39.6`},{19.224`,38.5`},{19.418`,38.98`},{19.418`,38.9`},{19.658`,39.`},{21.263`,39.24`},{22.961`,39.42`},{23.5`,39.65`},{23.5`,39.4`},{23.5`,39.13`},{23.5`,38.8`},{23.5`,38.9`},{23.6`,38.7`},{23.764`,39.21`},{23.764`,40.68`},{23.882`,39.`},{24.156`,39.59`},{25.294`,39.69`},{26.383`,39.77`},{27.6`,40.6`},{30.6`,40.22`},{30.6`,40.11`},{30.6`,39.91`},{30.6`,40.07`},{30.6`,40.2`},{30.6`,40.1`},{30.7`,40.1`},{30.8`,40.`},{35.2`,40.42`},{44.699`,41.7`},{44.7`,41.9`},{44.9`,41.89`},{45.2`,42.5`},{52.8`,43.01`},{52.8`,42.38`},{52.8`,42.1`},{52.8`,42.85`},{52.8`,42.71`},{52.9`,42.4`},{53.2`,42.9`},{62.3`,43.55`},{62.401`,43.1`},{62.5`,44.`},{62.7`,43.82`},{62.7`,42.96`},{62.8`,44.1`},{63.`,42.2`},{6164.1`,93.`},{7000.`,98.3`},{8000.`,101.7`},{8125.8`,101.`},{10712.`,117.`},{14121.`,104.`},{18615.`,93.`},{24539.`,124.`},{30000.`,120.`}};

dataRho = {{2.7675`,-0.426`},{3.1363`,-0.39`},{3.3208`,-0.38`},{3.3626`,-0.389`},{3.8446`,-0.3`},{3.8748`,-0.339`},{3.8748`,-0.4`},{4.0647`,-0.331`},{4.0739`,-0.29`},{4.3051`,-0.33`},{4.4208`,-0.351`},{4.5108`,-0.343`},{4.5108`,-0.345`},{4.5397`,-0.33`},{4.5397`,-0.31`},{4.5622`,-0.43`},{4.7209`,-0.26`},{4.923`,-0.29`},{5.3054`,-0.272`},{6.1554`,-0.245`},{6.1707`,-0.33`},{6.3073`,-0.205`},{6.8434`,-0.19`},{6.8598`,-0.157`},{7.1279`,-0.154`},{7.1673`,-0.32`},{7.3072`,-0.23`},{7.6213`,-0.183`},{8.7037`,-0.176`},{9.0296`,-0.194`},{9.7777`,-0.068`},{9.8995`,-0.157`},{9.967`,-0.153`},{9.9866`,-0.176`},{10.281`,-0.154`},{10.641`,-0.122`},{11.525`,-0.106`},{11.538`,-0.104`},{11.538`,-0.115`},{12.324`,-0.096`},{12.625`,-0.1216`},{12.625`,-0.1194`},{13.342`,-0.098`},{13.763`,-0.092`},{13.763`,-0.074`},{13.899`,-0.1029`},{13.966`,-0.0987`},{13.966`,-0.1024`},{15.373`,-0.024`},{16.523`,-0.064`},{16.83`,-0.04`},{16.83`,0.008`},{17.54`,-0.048`},{18.138`,-0.039`},{18.17`,-0.011`},{18.682`,-0.038`},{19.37`,-0.034`},{19.418`,0.019`},{19.848`,-0.0247`},{20.125`,-0.02`},{20.818`,-0.0176`},{21.424`,-0.013`},{21.7`,-0.041`},{22.171`,-0.009`},{22.494`,0.022`},{23.002`,0.0118`},{23.205`,0.0099`},{23.5`,0.02`},{23.5`,0.022`},{23.764`,-0.028`},{23.882`,-0.011`},{24.3`,0.009`},{25.612`,0.025`},{27.185`,0.039`},{27.361`,0.012`},{30.6`,0.03`},{30.6`,0.034`},{30.6`,0.042`},{44.7`,0.062`},{52.8`,0.077`},{52.9`,0.078`},{62.4`,0.095`},{7000,0.141},{8000,0.12}};

ReA = 
  -aP*(E^(bP*α0P))*(Cos[Pi*α0P/2])*s^(2*(α0P)) - af*(E^(bf*α0f))*
    (Cos[Pi*α0f/2])*s^(2*(α0f));

ImA = 
   aP*(E^(bP*α0P))*(Sin[Pi*α0P/2])*s^(2*(α0P)) + af*(E^(bf*α0f))* 
     (Sin[Pi*α0f/2])*s^(2*(α0f));

modelTC = 4*0.389*(Pi/s^2)*ImA;

modelRho = ReA/ImA;

fit = FindFit[dataTC, modelTC, {{α0P, 1.0808}, {α0f, 0.5}, aP,bP, af, bf}, s]

fit = FindFit[dataRho, modelRho, {{α0P, 1.0808}, {α0f, 0.5}, aP, bP, af, bf}, s]

I have two different datasets (dataTC and dataRho) and two different models (modelTC and modelRho) which have shared parameters (α0P, α0f, aP, bP, af, bf). I have to obtain a fit for dataTC using modelTC and a fit for dataRho using modelRho. so that the values of fitted parameters be the same in case of both fits.

I am able to fit the data separately, but the result is not satisfying for both data sets simultaneously.

I obtain the following when fitting dataTC using modelTC:

{α0P -> 1.08414, α0f -> 0.550223, aP -> 1.25628,bP -> 1.16465, af -> 1.5452, bf -> 3.84271}

enter image description here

enter image description here

I obtain the following when fitting dataRho using modelRho:

{α0P -> 1.09486, α0f -> 0.583092, aP -> 0.718297,bP -> 0.903191, af -> 1.99279, bf -> 3.05292}

enter image description here

enter image description here

These models would describe these datasets properly if they had the same parameter values for both models. Can I fit both models for the corresponding data simultaneously? How can I solve this problem?

$\endgroup$

marked as duplicate by MarcoB, m_goldberg, Anton Antonov, Sascha, happy fish Feb 20 '17 at 10:53

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

  • $\begingroup$ Note to question poser: Do not include irrelevant or assumed points in your title, such as "How do I solve this?" or "I need help" as these make it harder to find relevant questions in the future. $\endgroup$ – David G. Stork Feb 20 '17 at 0:00
  • 1
    $\begingroup$ You can formulate an appropriate minimization problem, in which the minimized function is derived from both datasets and both models. E.g. NMinimize[ Norm[Map[modelTC[params], dataTC[[All,1]] ] - dataTC] + Norm[Map[modelRho[params], dataRho[[All,1]] ] - dataRho], params, paramsRanges] $\endgroup$ – Anton Antonov Feb 20 '17 at 1:58
  • $\begingroup$ @MarcoB Agreed -- it is a duplicate. My proposal to use NMinimize is basically the same as Heike's answer in Simultaneously fitting multiple datasets. $\endgroup$ – Anton Antonov Feb 20 '17 at 3:58
  • $\begingroup$ @MarcoB For whatever it's worth, I don't agree about closing because there are bigger issues with the data and models that the potential duplicate does not address. $\endgroup$ – JimB Feb 20 '17 at 4:03
10
$\begingroup$

You have several issues with the data and the models. First, if you are going to simultaneously estimate all of the parameters, you need to account for the fact that both models involve different variances about the curves. (Whether or not those variances are constant across the predictor variable is another question but we won't need to worry about that because of bigger problems to come.)

To account for two different variances you (currently) need to use the LogLikelihood function and maximize the sum of the two likelihood functions (close to what @AntonAntonov suggested).

(* Get initial values for the standard deviations *)
fitTC = NonlinearModelFit[dataTC, 
   modelTC, {{α0P, 1.0808}, {α0f, 0.5}, aP, bP, af, bf},
    s];
σTCInitial = fitTC["EstimatedVariance"]^0.5;
fitRho = NonlinearModelFit[dataRho, 
   modelRho, {{α0P, 1.0808}, {α0f, 0.5}, aP, bP, af, 
    bf}, s];
σRhoInitial = fitRho["EstimatedVariance"]^0.5;

(* Construct log of the likelihood by adding the two log likelihoods *)
logL = LogLikelihood[NormalDistribution[0, σTC], 
    dataTC[[All, 2]] - (modelTC /. s -> dataTC[[All, 1]])] +
   LogLikelihood[NormalDistribution[0, σRho], 
    dataRho[[All, 2]] - (modelRho /. s -> dataRho[[All, 1]])]; 

(* Maximize the sum of the log likelihoods *)
sol = FindMaximum[{logL, σTC > 0 && σRho > 0},
  {{α0P, 1.0808}, {α0f, 0.4}, {aP, 1.3}, {bP, 1.2}, {af, 2}, {bf, 5},
  {σTC, σTCInitial}, {σRho, σRhoInitial}}]
(* {-339.98440796408477`,{α0P -> 1.0814972832958802`, α0f -> 0.41732447545678597`,
    aP -> 1.2881667623357842`,bP -> 1.1878643359708285`,af -> 2.1619418128981676`,
    bf -> 5.367892449067913`,σTC -> 2.1820904238485257`,σRho -> 0.08852423436184981`}} *)

(* Show results *)
Show[ListPlot[dataTC],
  Plot[{fitTC[s], modelTC /. sol[[2]]}, {s, 0, 80}, PlotLegends -> {"Single", "Combined"}]]

Show[ListPlot[dataRho],
Plot[{fitRho[s], modelRho /. sol[[2]]}, {s, 0, 80}, PlotLegends -> {"Single", "Combined"}]]

This results in the following fits:

Data and fits

It doesn't look good for the combined fit. That could be because the data do not support common parameters or other things going on.

Looking at the parameter estimator correlation matrices for the individual models suggests a big problem:

fitTC["CorrelationMatrix"] // MatrixForm
fitRho["CorrelationMatrix"] // MatrixForm

Parameter correlation matrices

We see that all of the correlations are very high in magnitude and many are either +1 or -1. This means that the model is overparameterized (and even more so for the Rho model). This issue is what needs to be fixed.

Update: Pinning down the overparameterization

We see in the that in ReA and ImA that af and bf are never separated from each other as are ap and af so there are really only two parameters to estimate rather than four parameters. We could define c1 and c2 in the following manner:

modelTC = 4*0.389*(Pi/s^2)*ImA /. {aP*(E^(bP*α0P)) -> c1, af*(E^(bf*α0f)) -> c2}
(* (4.88832 (c2 s^(2 α0f) Sin[(π α0f)/2] + c1 s^(2 α0P) Sin[(π α0P)/2]))/s^2 *)
fitTC = NonlinearModelFit[dataTC, modelTC, {α0P, α0f, c1, c2}, s];
fitTC["CorrelationMatrix"] // MatrixForm

Correlation matrix for fitTC

We see that while the correlations are still very large in magnitude, there are no longer any correlations of +1 or -1.

For fitRho, there is an additional dependency in that we can only estimate c2/c1. However, when the datasets are combined c1 and c2 can be estimated because of the model structure for fitTC.

In short, the models have built-in dependencies that need re-thinking. But even with removing the dependencies, the theory with common parameters doesn't seem plausible with the available data. Or your instruments need a tune-up. The best advice I can give is "Get thee to a statistician." (Yes, I'm a bit biased about that.)

$\endgroup$
  • $\begingroup$ Jim, thank you for contributing this answer. Indeed, I had not identified the correlation issue at all. After trying to provide a combined fit by constructing a joint sum of squares and minimizing that, I obtained bad results (as you did) and looked no deeper. Your approach is much superior. $\endgroup$ – MarcoB Feb 20 '17 at 4:34
  • $\begingroup$ @Jim Baldwin I can't tell you how much I appreciate your in depth answers. I have learned a great deal from the answers you have contributed to fitting. $\endgroup$ – Jack LaVigne Feb 20 '17 at 15:59
  • 1
    $\begingroup$ @JackLaVigne. Good. I'm glad I can occasionally contribute. While I've had a Mathematica license since version 1 on a MacPlus, it's only been in the last few years that I've started to take advantage of great functions available. This forum has helped me greatly. $\endgroup$ – JimB Feb 20 '17 at 18:19
  • $\begingroup$ @Jim Baldwin Thank you for your detailed answer. It was very helpful. It did not solve my problem completely, but I have obtained an acceptable fit using your advices, reducing the interval of the datasets and modifying the model. $\endgroup$ – István Szanyi Feb 25 '17 at 19:10

Not the answer you're looking for? Browse other questions tagged or ask your own question.