5
$\begingroup$

I am trying to fit a set of data with a specific model using the following code. As the figure shows, the fitting procedure is not great. First there is an error message, and as is obvious the fitting curve lacks the required number of peaks. Any suggestions as to how to best fit the data with the given model equation?

dataS = {{1.5040631016957404, 
    22.05446614925076}, {1.5214948122848406, 
    22.05446614925076}, {1.5377673786834203, 
    22.486895879419606}, {1.5535919973280015, 
    22.703110744504034}, {1.5771249660714466, 22.919325609588455}, 
       {1.5920936236176717, 23.351755339757304}, {1.6073491419697115, 
    23.784185069926153}, {1.620287200783405, 
    23.892292502468365}, {1.6360904891320525, 
    24.432829665179423}, {1.6494973106811417, 24.540937097721635}, 
       {1.66495981834507, 25.189581692974905}, {1.676981103969422, 
    25.730118855685962}, {1.689177243413443, 
    25.946333720770383}, {1.7063600475229623, 
    26.81119318110808}, {1.7219298010236765, 27.67605264144578}, 
       {1.737786304720355, 28.432804669241257}, {1.7508863658226992, 
    29.18955669703674}, {1.762126288255509, 
    30.5949533200855}, {1.7818844164451941, 
    31.24359791533877}, {1.7913990048481847, 32.4327796733031}, 
       {1.7956604034971266, 34.16249859397848}, {1.8107362979093506, 
    35.02735805431618}, {1.8238614375442779, 
    34.16249859397848}, {1.8316060123558566, 
    32.865209403471944}, {1.833830843046453, 31.4598127804232}, 
       {1.842784479227779, 30.270631022458865}, {1.8484250367218138, 
    28.64901953432569}, {1.8495572943867873, 
    27.02740804619251}, {1.858665538478947, 
    25.40579655805933}, {1.8690201396360508, 24.432829665179423}, 
       {1.8771538278817868, 23.459862772299516}, {1.8900793538510194, 
    23.027433042130667}, {1.9043844754728239, 
    22.91932560958846}, {1.9164718706122645, 
    23.459862772299516}, {1.9311808411251234, 24.108507367552786}, 
       {1.9448638139104728, 25.081474260432692}, {1.9536725713860417, 
    25.946333720770383}, {1.9754011663459174, 
    26.703085748565876}, {1.988409955436317, 
    27.24362291127693}, {1.9976185247944775, 26.919300613650297}, 
       {2.008247598975348, 25.946333720770383}, {2.017641260008336, 
    25.081474260432692}, {2.0298487330170984, 
    23.784185069926153}, {2.0380694470618033, 
    22.486895879419606}, {2.0463570182521065, 21.29771412145528}, 
       {2.0575125305410404, 20.432854661117588}, {2.070208760448426, 
    20.000424930948736}, {2.085940778620852, 
    20.000424930948736}, {2.1019137305460767, 
    20.21663979603316}, {2.1211091182481505, 20.973391823828642}, 
       {2.1482735570135274, 21.5139289865397}, {2.1652192677575655, 
    22.05446614925076}, {2.198323879071299, 
    22.91932560958846}, {2.222596642792378, 
    23.459862772299516}, {2.2357625855648067, 23.784185069926153}, 
       {2.262568026172226, 24.216614800095}, {2.283097746582417, 
    24.757151962806056}, {2.3022466842110143, 
    25.297689125517113}, {2.3270876381088414, 
    25.730118855685962}, {2.3543047675106723, 26.27065601839702}, 
       {2.382166082861296, 27.02740804619251}, {2.4106947283310736, 
    27.892267506530203}, {2.439914970220427, 
    28.865234399410102}, {2.465818249278822, 
    29.405771562121167}, {2.492277432741851, 29.946308724832228}, 
       {2.529864794115782, 29.946308724832228}, {2.55989253914495, 
    29.405771562121167}, {2.5751753139224256, 
    28.432804669241257}, {2.604047165517507, 
    27.67605264144578}, {2.6244175831132344, 26.486870883481448}, 
       {2.6404829252930715, 24.86525939534827}, {2.6567461665260383, 
    23.459862772299516}, {2.6684859646174215, 
    22.16257358179297}, {2.682711406542016, 
    20.75717695874422}, {2.70918918918919, 18.919350605526617}, 
       {2.7263125359799636, 17.18963168485123}, {2.738676622554785, 
    16.10855735942911}, {2.7586941268239316, 
    14.919375601464775}, {2.776451030803175, 
    13.513978978416027}, {2.799620029080939, 12.21668978790948}, 
       {2.8284682290473735, 10.27075600214967}, {2.849824976799784, 
    9.83832627198082}, {2.9074519307569133, 
    9.405896541811972}, {2.9357204381207054, 
    8.216714783847644}, {2.973301814567886, 7.351855323509947}, 
       {3.0028718216398964, 6.486995863172255}, {3.045271888519732, 
    5.405921537750136}, {3.0794356379551346, 
    4.432954644870231}, {3.1079632323585713, 
    3.6762026170747464}, {3.130519433275781, 2.162698561483778}, 
       {3.1699585216345945, 1.0816242360616586}, {3.2104040165031127, 
    0.32487220826618085}, {3.2589145824326042, -0.215664954444882}, \
{3.323482156888569, -0.10755752190266817}, {3.4136601387916805, 
    0.5410870733506021}, 
       {3.4805550917279287, 1.5140539662305077}, {3.575346757407733, 
    2.4870208591104133}, {3.648796861127818, 
    3.351880319448111}, {3.734543366932933, 
    3.78431004961696}, {3.7908067754888184, 4.000524914701382}, 
       {3.853703633481055, 4.000524914701382}, {3.9391724852514995, 
    4.108632347243596}, {4.017797352084466, 
    4.108632347243596}, {4.077480015018497, 
    4.108632347243596}, {4.133296702688033, 4.216739779785803}}; 


ep = 4.44; 
w0 = 0; 
w1 = 1.88; 
w2 = 2.03; 
w3 = 2.78; 
w4 = 2.91; 
w5 = 4.31; 
model = ep + 
   10^5*Re[((a0*wp^2)/(w0^2 - \[Omega]^2 - I*\[Omega]*b0))*
       r0 + ((a1*wp^2)/(w1^2 - \[Omega]^2 - I*\[Omega]*b1))*
       r1 + ((a2*wp^2)/(w2^2 - \[Omega]^2 - I*\[Omega]*b2))*
       r2 + ((a3*wp^2)/(w3^2 - \[Omega]^2 - I*\[Omega]*b3))*r3 + 
             ((a4*wp^2)/(w4^2 - \[Omega]^2 - I*\[Omega]*b4))*
       r4 + ((a5*wp^2)/(w5^2 - \[Omega]^2 - I*\[Omega]*b5))*r5]; 
result = NonlinearModelFit[
  dataS, {model, {0.03 > wp > 0.02, 5 > a0 > 0.5, 5 > a1 > 0.5, 
    5 > a2 > 0.5, 5 > a3 > 0.5, 5 > a4 > 0.5, 10 > a5 > 0.5, 
    5 > b0 > 0, 5 > b1 > 0, 5 > b2 > 0, 5 > b3 > 0, 5 > b4 > 0, 
    5 > b5 > 0, 
         1 > r0 > 0.1, 1 > r1 > 0.1, 1 > r2 > 0.1, 1 > r3 > 0.1, 
    1 > r4 > 0.1, 1 > r5 > 0.1}}, {wp, a0, a1, a2, a3, a4, a5, b0, 
   b1, b2, b3, b4, b5, r0, r1, r2, r3, r4, r5}, \[Omega], 
  MaxIterations -> 1500, Method -> {NMinimize}]
result["BestFitParameters"]
result["AdjustedRSquared"]
result["AIC"]
fitplot1 = 
 Show[ListPlot[dataS], Plot[result[f], {f, 1, 4}, PlotRange -> Full]]

enter image description here

Improve this question
$\endgroup$
1
  • $\begingroup$ The Newton method was not employed as it gives a negative value for the frequency wp, and also negative values for some other parameters which take on only positive values. This method can only be employed for unconstrained situations. $\endgroup$
    – thils
    Commented Aug 5, 2016 at 2:23

3 Answers 3

Reset to default
8
$\begingroup$

Brute force method:

result = NonlinearModelFit[
  dataS, {model, {0.03 > wp > 0.02, 5 > a0 > 0.5, 5 > a1 > 0.5, 
    5 > a2 > 0.5, 5 > a3 > 0.5, 5 > a4 > 0.5, 10 > a5 > 0.5, 
    5 > b0 > 0, 5 > b1 > 0, 5 > b2 > 0, 5 > b3 > 0, 5 > b4 > 0, 
    5 > b5 > 0, 1 > r0 > 0.1, 1 > r1 > 0.1, 1 > r2 > 0.1, 
    1 > r3 > 0.1, 1 > r4 > 0.1, 1 > r5 > 0.1}}, {wp, a0, a1, a2, a3, 
   a4, a5, b0, b1, b2, b3, b4, b5, r0, r1, r2, r3, r4, r5}, ω, 
  MaxIterations -> 20000, 
  Method -> {NMinimize, 
    Method -> {"RandomSearch", "SearchPoints" -> 10000}}]

enter image description here

result["BestFitParameters"]

 {wp -> 0.02, a0 -> 0.698604, a1 -> 0.503243, a2 -> 0.5, a3 -> 3.33348,
  a4 -> 2.85747, a5 -> 5.88554, b0 -> 0.000967767, b1 -> 0.143442,   b2 -> 0.102361, 
  b3 -> 0.564017, b4 -> 1.23437, b5 -> 1.12135,   r0 -> 0.778439, r1 -> 0.145413, 
  r2 -> 0.1, r3 -> 0.124165,   r4 -> 0.766377, r5 -> 0.504192}

result["AdjustedRSquared"]

0.993176

Update 2:

Fast and more elegant

result = NonlinearModelFit[
  dataS, {model, {0.03 > wp > 0.02, 5 > a0 > 0.5, 5 > a1 > 0.5, 
    5 > a2 > 0.5, 5 > a3 > 0.5, 5 > a4 > 0.5, 10 > a5 > 0.5, 
    5 > b0 > 0, 5 > b1 > 0, 5 > b2 > 0, 5 > b3 > 0, 5 > b4 > 0, 
    5 > b5 > 0, 1 > r0 > 0.1, 1 > r1 > 0.1, 1 > r2 > 0.1, 
    1 > r3 > 0.1, 1 > r4 > 0.1, 1 > r5 > 0.1}}, {wp, a0, a1, a2, a3, 
   a4, a5, b0, b1, b2, b3, b4, b5, r0, r1, r2, r3, r4, r5}, ω, 
  MaxIterations -> 2000, 
  Method -> {NMinimize, 
   Method -> {"DifferentialEvolution", "SearchPoints" -> 40, 
    "ScalingFactor" -> 0.95, "CrossProbability" -> 0.05, 
    "PostProcess" -> {FindMinimum, Method -> "QuasiNewton"}}}]

enter image description here

result["BestFitParameters"]

{wp -> 0.02, a0 -> 4.99856, a1 -> 0.656895, a2 -> 0.5, a3 -> 0.738923,
  a4 -> 3.57501, a5 -> 5.85331, b0 -> 4.2197*10^-7, b1 -> 0.14334, 
 b2 -> 0.102336, b3 -> 0.564501, b4 -> 1.23393, b5 -> 1.12183, 
 r0 -> 0.108618, r1 -> 0.111336, r2 -> 0.1, r3 -> 0.562351, 
 r4 -> 0.611452, r5 -> 0.506908}

result["AdjustedRSquared"]

0.993176

Reference:

http://reference.wolfram.com/language/tutorial/ConstrainedOptimizationOverview.html http://reference.wolfram.com/language/tutorial/ConstrainedOptimizationGlobalNumerical.html#24713453

Improve this answer
$\endgroup$
3
  • $\begingroup$ Certainly faster than the brute force method. What is the meaning of "CrossProbability" $\endgroup$
    – thils
    Commented Aug 5, 2016 at 3:25
  • $\begingroup$ @thils look specifically here: reference.wolfram.com/language/tutorial/… $\endgroup$
    – Young
    Commented Aug 5, 2016 at 3:30
  • $\begingroup$ Useful, increasing searchpts helps $\endgroup$
    – thils
    Commented Aug 5, 2016 at 3:51
5
$\begingroup$

I suspect that you need to make predictions from the fitting function outside of Mathematica. If not, then just using the Mathematica interpolation functions will provide a perfect fit and (by definition) allow for interpolation. This wouldn't give you standard errors for any parameters but because there is no a priori model form and you have no repeated observations to obtain a measure of pure fit error, I don't see that such standard errors or prediction intervals would have much meaning.

But if you just need a fit that looks good to the eye and can be "easily" programmed in whatever outside package is to be used, you might consider kernel regression.

n = Length[dataS[[All, 1]]];
{xMin, xMax} = MinMax[dataS[[All, 1]]];

(* Kernel regression function *)
f[x_, data_, n_, h_] := 
 Sum[dataS[[i, 2]] Exp[-(x - dataS[[i, 1]])^2/h^2], {i, n}]/
  Sum[Exp[-(x - dataS[[i, 1]])^2/h^2], {i, n}]

h = 0.01; (* Initial bandwidth *)

ListPlot[{dataS, 
  Table[{x, f[x, dataS, n, h (1 + 5 (x - xMin)/(xMax - xMin))]},
   {x, xMin, xMax, 0.005}]}, Joined -> {False, True}]

Fit with kernel regression

Note that I've adjusted the bandwidth to depend on x because the bandwidth needs to be small at the beginning of the sequence to fit the sharp peaks and larger towards the end of the sequence.

Improve this answer
$\endgroup$
6
  • $\begingroup$ It is interesting that the data doesn't look noisy at all... OTOH, OP did specify a model, and maybe he knows something we don't. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Aug 5, 2016 at 4:46
  • $\begingroup$ @J.M. Agree. Little or no noise makes me suspicious. But then again, I work mainly with biological systems where everything has lots of noise. (And I've heard rumors that engineers and physicists (i.e., what seems to be the majority of Mathematica users) don't believe in noise. $\endgroup$
    – JimB
    Commented Aug 5, 2016 at 4:51
  • 1
    $\begingroup$ Nice fit, just too good! The original model incorporates the five experimentally observed absorption peaks as specified by w0 to w5, hence the equation is a summation of several oscillators. $\endgroup$
    – thils
    Commented Aug 5, 2016 at 5:35
  • $\begingroup$ I say "no apriori model" because you've looked at the data to count peaks prior to constructing the model. (If you saw just 3 "peaks", I'm sure you'd start with a less complicated model. But please correct me if I'm wrong.) I also think that there's a potential identifiability issue for such tasks. Knowing the sum is 35 does not tell you if that's from 10 + 25 or 7.3 + 27.6 or 1 + 7 + 3 + 24. The linear part between 2.2 and 2.5 could be the result of two "peaks"/"kernels" that are close to each other that when summed together don't show a peak. Do you need a good fit everywhere? $\endgroup$
    – JimB
    Commented Aug 5, 2016 at 15:54
  • 1
    $\begingroup$ @Baldwin Can't hide peaks in the spectrum, so its an empirical model, with one of the best fit parameters (@Young above) nicely matching the expt. value. Perhaps the interpolation procedure will be useful for other systems? $\endgroup$
    – thils
    Commented Aug 5, 2016 at 22:26
3
$\begingroup$

Suggest to use Newton's method as follows:

NonlinearModelFit[dataS, model, {wp, a0, a1, a2, a3, a4, a5, b0, b1, b2, b3, b4, b5, r0, r1, r2, r3, r4, r5}, ω ,  Method -> "Newton"]

enter image description here

Improve this answer
$\endgroup$
2
  • $\begingroup$ The Newton method gives a negative value for the frequency wp, and also negative values for some other parameters which take on only positive values. $\endgroup$
    – thils
    Commented Aug 5, 2016 at 2:23
  • $\begingroup$ Okay, I see @Young has probably helped you out already. $\endgroup$
    – Alexander B.
    Commented Aug 5, 2016 at 3:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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