NonlinearModelFit with some shared parameters

Question

I have a big set of data (more than 800 spectra) which need to be fit with a Maxwell-Boltzmann function. Let's start with only two spectra

    data1={1.53108, -0.00184597}, {1.53503, -0.00131884}, {1.539, 
  0.00185461}, {1.54299, -0.00419277}, {1.54701, -0.00639723}, \
{1.55104, -0.00551708}, {1.5551, -0.00800965}, {1.55917, \
-0.00677264}, {1.56327, -0.0110927}, {1.56739, -0.0103286}, {1.57153, 
-0.00973164}, {1.5757, -0.0111152}, {1.57988, -0.0113855}, {1.58409, 
-0.0117398}, {1.58832, -0.0127097}, {1.59257, -0.0139916}, {1.59685, 
-0.014101}, {1.60115, -0.0139898}, {1.60547, -0.0118756}, {1.60981, 
-0.01029}, {1.61418, -0.00805921}, {1.61858, -0.000821386}, {1.62299, 
  0.00466453}, {1.62743, 0.00821456}, {1.6319, 0.0171329}, {1.63639, 
  0.0310873}, {1.6409, 0.0413745}, {1.64544, 0.053401}, {1.65001, 
  0.0616692}, {1.6546, 0.0702816}, {1.65921, 0.0752915}, {1.66385, 
  0.0843292}, {1.66852, 0.0806191}, {1.67322, 0.0865337}, {1.67794, 
  0.0820241}, {1.68268, 0.0869062}, {1.68746, 0.0808511}, {1.69226, 
  0.0864214}, {1.69709, 0.0748011}, {1.70194, 0.0822971}, {1.70683, 
  0.0768256}, {1.71174, 0.0783893}, {1.71668, 0.0702391}, {1.72165, 
  0.0791383}, {1.72665, 0.0711978}, {1.73168, 0.0772164}, {1.73673, 
  0.069743}, {1.74182, 0.0742699}, {1.74693, 0.0635676}, {1.75208, 
  0.0705165}, {1.75726, 0.0580984}, {1.76246, 0.0649913}, {1.7677, 
  0.0537584}, {1.77297, 0.0561502}, {1.77827, 0.0471154}, {1.7836, 
  0.0636209}, {1.78897, 0.0396738}, {1.79437, 0.0417097}, {1.7998, 
  0.036672}, {1.80526, 0.0483584}, {1.81076, 0.0359746}, {1.81629, 
  0.0478157}, {1.82185, 0.0356168}, {1.82745, 0.0382669}, {1.83308, 
  0.0287583}, {1.83875, 0.0391059}, {1.84445, 0.0300186}, {1.85019, 
  0.0283772}, {1.85596, 0.0245703}, {1.86177, 0.0263533}, {1.86761, 
  0.0170249}, {1.8735, 0.0268241}, {1.87942, 0.0134921}, {1.88537, 
  0.0148601}, {1.89137, 0.00781962}, {1.8974, 0.0119102}, {1.90348, 
  0.00445442}, {1.90959, 0.0074297}, {1.91574, 0.00365792}, {1.92193, 
  0.0108381}, {1.92816, -0.00316959}, {1.93443, -0.000539258}, 
{1.94074, 0.0043318}, {1.9471, -0.00430929}, {1.95349, 
  0.00441513}, {1.95993, 
  0.00595201}, {1.96641, -0.00438416}, {1.97293, 
  0.00600895}, {1.9795, -0.0021935}, {1.98611, -0.0121693}, {1.99276, 
-0.0180589}, {1.99946, -0.000531422}, {2.00621, -0.0103922}, {2.013, \
-0.0134376}, {2.01983, -0.00937648}, {2.02671, -0.0123256}, {2.03364, \
-0.0171005}, {2.04062, -0.012415}, {2.04765, -0.00842077}, {2.05472, 
-0.0222206}, {2.06184, -0.0173176}, {2.06902, -0.00691508}, {2.07624, 
-0.0136626}, {2.08351, -0.00384277}, {2.09084, -0.0133483}, {2.09821, 
-0.0214737}, {2.10564, -0.0107265}, {2.11312, -0.0168493}, {2.12066, 
-0.0148875}, {2.12824, -0.0204143}, {2.13589, -0.00490166}, {2.14358, 
-0.0163762}, {2.15134, -0.0132914}, {2.15915, -0.0119276}, {2.16702, 
-0.00761709}, {2.17494, -0.00623985}, {2.18292, -0.0139031}, 
{2.19096, -0.00510826}, {2.19906, -0.00934085}}

and

  data2={1.53108, 0.00187441}, {1.53503, 0.00891081}, {1.539, 
  0.0111849}, {1.54299, 0.004364}, {1.54701, 0.00628025}, {1.55104, 
  0.00680789}, {1.5551, 0.00853707}, {1.55917, 0.00920825}, {1.56327, 
  0.00928874}, {1.56739, 0.0117898}, {1.57153, 0.0114837}, {1.5757, 
  0.0139937}, {1.57988, 0.0161987}, {1.58409, 0.0164855}, {1.58832, 
  0.012975}, {1.59257, 0.0112959}, {1.59685, 0.0197188}, {1.60115, 
  0.021546}, {1.60547, 0.024632}, {1.60981, 0.0281438}, {1.61418, 
  0.0310039}, {1.61858, 0.0405104}, {1.62299, 0.047554}, {1.62743, 
  0.0507924}, {1.6319, 0.0578233}, {1.63639, 0.0714587}, {1.6409, 
  0.0806594}, {1.64544, 0.092654}, {1.65001, 0.101177}, {1.6546, 
  0.110297}, {1.65921, 0.114115}, {1.66385, 0.123051}, {1.66852, 
  0.116146}, {1.67322, 0.121149}, {1.67794, 0.115578}, {1.68268, 
  0.119243}, {1.68746, 0.108077}, {1.69226, 0.117215}, {1.69709, 
  0.0994373}, {1.70194, 0.102059}, {1.70683, 0.0951892}, {1.71174, 
  0.0951288}, {1.71668, 0.0837428}, {1.72165, 0.0940461}, {1.72665, 
  0.0814388}, {1.73168, 0.084036}, {1.73673, 0.0717605}, {1.74182, 
  0.0774831}, {1.74693, 0.0627288}, {1.75208, 0.0652405}, {1.75726, 
  0.0613616}, {1.76246, 0.0609079}, {1.7677, 0.0501096}, {1.77297, 
  0.0564122}, {1.77827, 0.0373842}, {1.7836, 0.0519737}, {1.78897, 
  0.0257087}, {1.79437, 0.0319259}, {1.7998, 0.0243793}, {1.80526, 
  0.0278366}, {1.81076, 0.021786}, {1.81629, 0.0273007}, {1.82185, 
  0.0157447}, {1.82745, 0.0144138}, {1.83308, 0.0149415}, {1.83875, 
  0.00669958}, {1.84445, 0.00744026}, {1.85019, 0.00570051}, {1.85596,
   0.00109972}, {1.86177, 
  0.00250912}, {1.86761, -0.00782543}, {1.8735, -0.00982973}, 
{1.87942, -0.00647511}, {1.88537, -0.00960983}, {1.89137, 
-0.00156248}, {1.8974, -0.0115654}, {1.90348, -0.0117897}, {1.90959, 
-0.0053029}, {1.91574, -0.0144718}, {1.92193, -0.00995421}, {1.92816, 
-0.0120168}, {1.93443, -0.0233383}, {1.94074, -0.0168577}, {1.9471, 
-0.0257402}, {1.95349, -0.0185922}, {1.95993, -0.0176849}, {1.96641, 
-0.0223636}, {1.97293, -0.0259812}, {1.9795, -0.0218552}, {1.98611, 
-0.0207055}, {1.99276, -0.0189714}, {1.99946, -0.0136022}, {2.00621, 
-0.0232702}, {2.013, -0.0199806}, {2.01983, -0.0220758}, {2.02671, 
-0.0211568}, {2.03364, -0.0200752}, {2.04062, -0.0334419}, {2.04765, 
-0.0262612}, {2.05472, -0.0285296}, {2.06184, -0.0262901}, {2.06902, 
-0.0148341}, {2.07624, -0.0239107}, {2.08351, -0.0229732}, {2.09084, 
-0.0159384}, {2.09821, -0.0150739}, {2.10564, -0.013385}, {2.11312, 
-0.014881}, {2.12066, -0.0108471}, {2.12824, -0.0180245}, {2.13589, 
-0.0149589}, {2.14358, -0.01849}, {2.15134, -0.010917}, {2.15915, 
-0.0111556}, {2.16702, -0.00881686}, {2.17494, -0.00877528}, 
{2.18292, -0.0162301}, {2.19096, -0.00540098}, {2.19906, -0.00319884}}

I first define the function (k is a constant)

f[x_, A_, Ef_, T_, d_] := A*Exp[-(x - Ef)/(k*T)] + d;

Only Ef is the shared parameter among all of these 800 spectra. As previous answer, I started from here

    MultiNonlinearModelFit[datasets_, expressions_, params_, 
   constraints : _ : True, independents_Symbol, 
   opts : OptionsPattern[]] := 
  MultiNonlinearModelFit[datasets, expressions, constraints, 
   params, {independents}, opts];

MultiNonlinearModelFit[datasets : {__?(MatrixQ[#, NumericQ] &)}, 
    expressions_List, 
    constraints : _ : 
     True, {fitParams__Symbol}, {independents__Symbol}, 
    opts : OptionsPattern[]] /; 
   Length[expressions] === Length[datasets] := 
  Module[{fitfun, numSets = Length[expressions], 
    augmentedData = 
     Catenate@
      MapIndexed[
       Join[(*Attach indices to the data*)
         ConstantArray[N[#2], Length[#1]], #1, 2] &, datasets]}, 
   fitfun = 
    With[{conditions = 
       Map[{\[FormalN] == #, expressions[[#]]} &, N@Range[numSets]]}, 
     Which @@ Catenate[conditions]];
   NonlinearModelFit[augmentedData, 
    If[TrueQ[constraints], 
     fitfun, {fitfun, constraints}], {fitParams}, {\[FormalN], 
     independents},(*use dataset index as extra independent variable*)
    opts]];

Then, I applied

k= 0.0000651

fit = MultiNonlinearModelFit[{data1, 
    data2}, {amp1 *Exp[(-x + sharedOffset)/(k*T1)]+d1, 
    amp2 *Exp[(-x + sharedOffset)/(k*T2)]+d2}, {amp1,amp2, T1, T2, sharedOffset, d1, d2}, {x}];
fit["BestFitParameters"];

I chose just two spectra to start. The result of the fit is totally a mess. What I am doing wrong? Few question on this approach (I am new, sorry)

1) Is it possible to choose starting parameters for all of them? This would improve the fit a lot. 2) Is it possible to extend it to 800 spectra without defining 800 times A1,A2,A3... and so on?

What is the value of $k$? Setting it to any value gives you results. But your model is overparameterized as some parameter estimators are perfectly correlated with others. In short, the model is too complex for the data. Specifically the estimators for amp1 and d1 cannot be estimated separately. Same thing for amp2 and d2. — JimB, Jun 11 '20 at 14:14
And to show my ignorance of physics, do you not want (-x+sharedOffset)^2 rather than (-x+sharedOffset)? — JimB, Jun 11 '20 at 14:55
Add ClearAll[MultiNonlinearModelFit before your definitions, then in the definition of the fitter function change fitParams__Symbol to simply fitParams__. Finally, use the form: {{amp1, 10}, {amp2, 20}, T1, T2, sharedOffset, d1, d2} to give initial values for certain parameters. You can mix and match, giving initial values for some but not all. — MarcoB, Jun 11 '20 at 14:55
For the second question, you will want to automatically generate lists of parameters. Those would be better formatted as indexed variables, rather than string symbols (i.e. amp[1] rather than amp1 etc). To generate those, look into Array: for instance, consider Array[amp,10]. — MarcoB, Jun 11 '20 at 14:59
From a physical meaning perspective, what do you mean by spectra here? Your data looks like time decays. In fact, your two samples look like linear time decays. you will have a bad time trying to fit those to exponential decays as in your model though. As @JimB hinted at as well, your model may not only be statistically overdetermined, but it may crucially be overly complex for the features in your data. Perhaps you could show a more representative data set? — MarcoB, Jun 11 '20 at 15:02

JimB · Accepted Answer · 2020-06-12 04:17:04Z

I want to try to convince you that what you want to do can't be done.

The basic model is

amp*Exp[(-x+sharedOffset)/(k*T)] + d

For a single dataset both amp and sharedOffset cannot be uniquely identified. That is because the model can be written in the following forms that result in identical predictions:

amp*Exp[(-x + sharedOffset)/(k*T)] + d
amp*Exp[sharedOffset/(k*T)]*Exp[-x/(k*T)] + d
ampShared*Exp[-x/(k*T)] + d

There are only three parameters that can be estimated rather than four. For example, one can estimate d, T, and the product amp*Exp[sharedOffset/(k*T)] but one cannot obtain separate estimates of amp and sharedOffset.

But what about with 2 data sets with a common value of sharedOffset? The answer is "No" here, too. Why? We can estimate amp1*Exp[sharedOffset/(k T1)] and amp2*Exp[sharedOffset/(k T2). That gives us two equations and three unknowns (with the unknowns being amp1, amp2, and sharedOffset. Again, not enough information to estimate all three parameters.

But you say "Why does MultiNonlinearModelFit give plausible-looking results?" Here are those results:

(* Select subsets of the complete datasets *)
subset1 = Select[data1, 1.7 <= #[[1]] <= 2 &];
subset2 = Select[data2, 1.7 <= #[[1]] <= 2 &];

k = 1 ;  (* Setting k=1 only affects the scale of T and 
in this case results in no need to set starting values for
the parameters other than the default parameter values of 1 *)

fit = MultiNonlinearModelFit[{subset1, subset2},
  {amp1*Exp[(-x + sharedOffset)/(k*T1)] + d1,
   amp2*Exp[(-x + sharedOffset)/(k*T2)] + d2},
  {amp1, amp2, T1, T2, sharedOffset, d1, d2}, x];
fit["ParameterConfidenceIntervalTable"]

We see estimates, standard errors, and "nice" confidence intervals. sharedOffset has a 95% confidence interval of {1.3624, 1.44574}. The root mean square error is small:

fit["EstimatedVariance"]^0.5
(* 0.0058183 *)

But standard diagnostics suggest there might be issues. Consider the correlation matrix:

fit["CorrelationMatrix"] // MatrixForm

We see correlations of almost -1: -0.993208 and -0.991959. The rank of the covariance matrix is only 6 (rather than 7 for the 7 parameters):

MatrixRank[fit["CovarianceMatrix"]]
(* 6 *)

But what about the "nice" confidence interval for sharedOffset? Suppose we plug in a value for sharedOffset way outside the confidence interval. If we choose 3 for that value we actually get a (slightly) smaller root mean square error (admittedly probably because of internal precision losses) and a full rank covariance matrix:

fit2 = MultiNonlinearModelFit[{subset1, subset2},
  {amp1*Exp[(-x + 3)/(k*T1)] + d1,
   amp2*Exp[(-x + 3)/(k*T2)] + d2},
  {amp1, amp2, T1, T2, d1, d2}, x,
  MaxIterations -> 1000];

fit["EstimatedVariance"]^0.5
(* 0.00578914 *)

MatrixRank[fit["CovarianceMatrix"]]
(* 6 *)

Should Mathematica give some warnings for the original model (especially about the rank of the covariance matrix)? Probably. But one needs to look at the fit diagnostics and prior to that take a careful look at the structure of the model to be fit.

In short if one can't estimate the shared parameter with an individual data set, one should be very cautious when attempting to fit that shared parameter for multiple data sets. There are just some models for which this is the case and I'm strongly suggesting that yours is one of those models. (The MultiNonlinearModelFit also assumes that there's a common error variance which seems a bit restrictive to me.)

I see what you mean, @JimB. What is your suggestion? I have 800 curves from the same dataset (each of them represents a situation in time) and I need to get the amplitude, temperature out of it, considering that the curve sometimes goes below zero (that's why the d parameter). The sharedoffset is the Fermi energy and it should be common to all of these curves. What's surprise me is that this is a really common model for this stuff, so it should not so crazy to apply — Clarine, Jun 12 '20 at 6:19
If you just need the amplitude and temperature, then I would recommend just running all 800 separately (with setting sharedOffset to zero). If the other 798 curves all have similar numbers of data points, then estimating 800 separate error variances is not a problem. (Again, MultiNonlinearModelFit assumes a single common error variance but with enough data as you seem to have with each curve, that would seem an unnecessary assumption.) — JimB, Jun 12 '20 at 15:26
I believe you when you say this is a common model. Do you have a reference or two for me to read? — JimB, Jun 12 '20 at 20:09

Stack Exchange Network

current community

your communities

more stack exchange communities

NonlinearModelFit with some shared parameters

1 Answer 1

Your Answer

Not the answer you're looking for? Browse other questions tagged calculus-and-analysis fitting nonlinear or ask your own question.

Hot Network Questions

NonlinearModelFit with some shared parameters

1 Answer 1

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged calculus-and-analysis fitting nonlinear or ask your own question.

Related

Hot Network Questions