1
$\begingroup$

I tried to fit the curve of the dataset. I used 'Nonlinearmodelfit' function to get the specific fitted equation. Yes, I can get the result, but I am wondering there is a way to improve the result.

This is the dataset, basically {Angle, Intensity}

data={{0.,4852.},{0.0872665,3128.},{0.174533,2686.},{0.261799,6450.},{0.349066,15532.},{0.436332,28730.},{0.523599,42838.},{0.610865,52648.},{0.698132,51744.},{0.785398,66190.},{0.872665,71260.},{0.959931,73682.},{1.0472,67820.},{1.13446,53786.},{1.22173,36626.},{1.309,26750.},{1.39626,28354.},{1.48353,26258.},{1.5708,19888.},{1.65806,19844.},{1.74533,16350.},{1.8326,15732.},{1.91986,12062.},{2.00713,9996.},{2.0944,15350.},{2.18166,17160.},{2.26893,21986.},{2.35619,19248.},{2.44346,22100.},{2.53073,20096.},{2.61799,16962.},{2.70526,15510.},{2.79253,14280.},{2.87979,11278.},{2.96706,7664.},{3.05433,5452.},{3.14159,3114.},{3.22886,1712.},{3.31613,1860.},{3.40339,3510.},{3.49066,7498.},{3.57792,19636.},{3.66519,38690.},{3.75246,56766.},{3.83972,68802.},{3.92699,64136.},{4.01426,74730.},{4.10152,83330.},{4.18879,63400.},{4.27606,48962.},{4.36332,35622.},{4.45059,22414.},{4.53786,15704.},{4.62512,11696.},{4.71239,10184.},{4.79966,9836.},{4.88692,11946.},{4.97419,11496.},{5.06145,13668.},{5.14872,18054.},{5.23599,31018.},{5.32325,37838.},{5.41052,47124.},{5.49779,55154.},{5.58505,46564.},{5.67232,34936.},{5.75959,27916.},{5.84685,22714.},{5.93412,15792.},{6.02139,12696.},{6.10865,8584.},{6.19592,8166.}}

And, These below are the model and parameters that I used.

Model= a + b*Cos[x] + c*Cos[2x] + d*Cos[3x] + e*Cos[4x] + f*Cos[5x] + g*Cos[6x] + h*Sin[x] + i*Sin[2x] + j*Sin[3x] + k*Sin[4x] + l*Sin[5x] + m*Sin[6x]
coeff = {a, b, c, d, e, f, g, h, i, j, k, l, m};
nlm = Normal[NonlinearModelFit[data, model , coeff , x, MaxIterations -> 1000]];

Thank you

Improve this question
$\endgroup$
3
  • 2
    $\begingroup$ Another approach would be to use a Fourier transform. $\endgroup$ – bill s Jul 9 '20 at 0:43
  • $\begingroup$ bill's approach is the most efficient if your points are equispaced; otherwise, you would really need to do full least squares. $\endgroup$ – J. M.'s torpor♦ Jul 9 '20 at 4:12
  • 2
    $\begingroup$ I think this is a linear model (linear in the unknown coefficients) $\endgroup$ – mikado Jul 9 '20 at 6:09
5
$\begingroup$

Mikado's comment to the original question is on the money: this is a linear model; not a non-linear one. This makes the fit much easier and faster with LinearModelFit, since you don't have to worry about local minima. With linear models you normally don't have to fudge around nearly as much:

basisFun = {
   1, Cos[x], Cos[2 x], Cos[3 x], Cos[4 x], Cos[5 x],
   Cos[6 x], Sin[x], Sin[2 x], Sin[3 x], Sin[4 x], Sin[5 x], Sin[6 x]
};
fit = LinearModelFit[data, basisFun, x];
Show[ListPlot[data], Plot[fit[x], {x, 0, 2 π}]]

enter image description here

Edit

Bonus content:

fourrierFit[data_, n_Integer, k : _?NumericQ : 1] := LinearModelFit[
   data,
   Flatten[{Sin[# k \[FormalX]], Cos[# k \[FormalX]]} & /@ Range[n]],
   \[FormalX]
];
Manipulate[
 Show[ListPlot[data],
  Plot[
   Evaluate @ Map[
     Legended[fourrierFit[data, #, k]["BestFit"], #] &,
     Range[6]
     ],
   {\[FormalX], 0, 2 π}
   ]
 ],
 {{k, 1}, 0.1, 5}
]

If you want nice, Bayesian, error bars take a look at my ResourceFunction:

k = 1;
n = 6
basisFuns = Flatten[{Sin[# k \[FormalX]], Cos[# k \[FormalX]]} & /@ Range[n]];
fit = ResourceFunction["BayesianLinearRegression"][data, basisFuns, \[FormalX]];
Show[
 ListPlot[data],
 Plot[
  Evaluate @ InverseCDF[fit["Posterior", "PredictiveDistribution"], {0.95, 0.5, 0.05}],
  {\[FormalX], 0, 2 π},
  Filling -> {1 -> {2}, 3 -> {2}},
  PlotLegends -> Table[Quantity[i, "Percent"], {i, {95, 50, 5}}]
  ],
 Plot[
  Evaluate @ InverseCDF[fit["Posterior", "UnderlyingValueDistribution"], {0.95, 0.5, 0.05}],
  {\[FormalX], 0, 2 π},
  Filling -> {1 -> {2}, 3 -> {2}},
  PlotStyle -> Dashed
 ]
]

enter image description here

Improve this answer
$\endgroup$
1
  • $\begingroup$ +1 Mathematica needs more Bayesian exposure. $\endgroup$ – JimB Jul 9 '20 at 16:14
2
$\begingroup$

If you're going to stick with pairs of sines and cosines, then using AICc is a reasonable approach to determine an appropriate number of pairs.

Rewriting your code to allow for any number of pairs of sines and cosines:

data = {{0., 4852.}, {0.0872665, 3128.}, {0.174533, 2686.}, {0.261799, 6450.}, {0.349066, 15532.}, {0.436332, 28730.}, {0.523599, 42838.}, {0.610865, 52648.}, {0.698132, 51744.}, {0.785398, 66190.}, {0.872665, 71260.}, {0.959931, 73682.}, {1.0472, 67820.}, {1.13446, 53786.}, {1.22173, 36626.}, {1.309, 26750.}, {1.39626, 28354.}, {1.48353, 26258.}, {1.5708, 19888.}, {1.65806, 19844.}, {1.74533, 16350.}, {1.8326, 15732.}, {1.91986, 12062.}, {2.00713, 9996.}, {2.0944, 15350.}, {2.18166, 17160.}, {2.26893, 21986.}, {2.35619, 19248.}, {2.44346, 22100.}, {2.53073, 20096.}, {2.61799, 16962.}, {2.70526, 15510.}, {2.79253, 14280.}, {2.87979, 11278.}, {2.96706, 7664.}, {3.05433, 5452.}, {3.14159, 3114.}, {3.22886, 1712.}, {3.31613, 1860.}, {3.40339, 3510.}, {3.49066, 7498.}, {3.57792, 19636.}, {3.66519, 38690.}, {3.75246, 56766.}, {3.83972, 68802.}, {3.92699, 64136.}, {4.01426, 74730.}, {4.10152, 83330.}, {4.18879, 63400.}, {4.27606, 48962.}, {4.36332, 35622.}, {4.45059, 22414.}, {4.53786, 15704.}, {4.62512, 11696.}, {4.71239, 10184.}, {4.79966, 9836.}, {4.88692, 11946.}, {4.97419, 11496.}, {5.06145, 13668.}, {5.14872, 18054.}, {5.23599, 31018.}, {5.32325, 37838.}, {5.41052, 47124.}, {5.49779, 55154.}, {5.58505, 46564.}, {5.67232, 34936.}, {5.75959, 27916.}, {5.84685, 22714.}, {5.93412, 15792.}, {6.02139, 12696.}, {6.10865, 8584.}, {6.19592, 8166.}};

n = 6;
model = cs0 + Sum[c[i] Cos[i x] + s[i] Sin[i x], {i, n}];
coeff = Flatten[{cs0, Table[{c[i], s[i]}, {i, n}]}];
nlm = NonlinearModelFit[data, model, coeff, x]
Show[ListPlot[data], Plot[nlm[x], {x, 0, 2 π}]]

Data and initial fit

I assume this is the fit you found.

We can loop through and find the number of pairs with the smallest AICc value.

nmax = 25;
aicc = ConstantArray[{0, 0}, nmax];
Do[
 model = cs0 + Sum[c[i] Cos[i x] + s[i] Sin[i x], {i, n}];
 coeff = Flatten[{cs0, Table[{c[i], s[i]}, {i, n}]}];
 aicc[[n]] = {n, NonlinearModelFit[data, model, coeff, x]["AICc"]},
 {n, nmax}]

ListPlot[aicc, Frame -> True,
  FrameLabel -> (Style[#, Bold, 16] &) /@ {"Number of sine and cosine pairs", 
  "\!\(\*SubscriptBox[\(AIC\), \(c\)]\)"}]

Number of pairs of sines and cosines vs AICc

So we see that n = 10 gives the lowest AICc value so we should not have any more than 10 pairs of sines and cosines. That fit is as follows:

n = 10;
model = cs0 + Sum[c[i] Cos[i x] + s[i] Sin[i x], {i, n}];
coeff = Flatten[{cs0, Table[{c[i], s[i]}, {i, n}]}];
nlm = NonlinearModelFit[data, model, coeff, x];
Show[ListPlot[data], Plot[nlm[x], {x, 0, 2 π}]]

Better fit

Improve this answer
$\endgroup$
1
$\begingroup$

You can do Fourier analysis and decide, which Sin and Cos frequences over a choosen threshold you want to take into account and see the effect on fit. (Thanks to the ideas from other contributors!)

fou = Fourier[data[[All, 2]]];
tabc = Table[Cos[i x], {i, 1, Round[Length[data]/2]}];
tabs = Table[Sin[i x], {i, 1, Round[Length[data]/2]}];    
Manipulate[
 lp1 = {refou = Re[Rest@fou] /. (x_ /; Abs[x] < th) -> {}, 
 imfou = Im[Rest@fou] /. (x_ /; Abs[x] < th) -> {}}; 
 rc := refou /. (x_?NumericQ -> (x/x));
 rs := imfou /. (x_?NumericQ -> (x/x));
 cfun := DeleteCases[tabc rc[[1 ;; Round[Length[data]/2]]], {}];
 sfun := DeleteCases[tabs rs[[1 ;; Round[Length[data]/2]]], {}];
 basisFun = Join[{1}, cfun, sfun];
 fit = LinearModelFit[data, basisFun, x];
 {ListPlot[lp1, ImageSize -> 300, Filling -> Axis, 
    PlotRange -> {{0, Round[Length[data]/2]}, All}, 
    PlotStyle -> {Blue, Red}], 
  Plot[fit[x], {x, 0, 2 \[Pi]}, Epilog -> {Red, Point@data}, 
   ImageSize -> 300, 
   PlotRange -> {0, 1.1 Max[data[[All, 2]]]}]}, {{th, 10000}, 1000, 
  20000, Appearance -> "Labeled"}]

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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