# How to fit the data?

Here is my attempt.

data = {{0., 2.61}, {0.1, 2.62}, {0.2, 2.62}, {0.3, 2.62}, {0.4,
2.63}, {0.5, 2.63}, {0.6, 2.74}, {0.7, 2.98}, {0.8, 3.66}, {0.9,
5.04}, {1., 7.52}, {1.1, 10.74}, {1.2, 12.62}, {1.3, 10.17}, {1.4,
5}, {1.5, 2.64}, {1.6, 11.5}, {1.65, 35.4}};
NonlinearModelFit[data,a*Cosh[b*x^c*Sin[d*x^e]], {{a, 3}, {b, 0.2}, {c, 1}, {d, 1},
{e, 3}}, x, Method -> "LevenbergMarquardt"]


NonlinearModelFit::nrjnum: The Jacobian is not a matrix of real numbers at {a,b,c,d,e} = {3.,0.2,1.,1.,3.}.

The same issues appear with other initial points. As far as I know it, a good fit with the residuals of quantity $$0.01$$ exists and can be found with Mathcad (not MATLAB).

• The Jacobian contains Log[x] and your data contains x == 0 Commented Jan 28, 2019 at 17:09
• @Coolwater: There are other methods of fitting. Commented Jan 28, 2019 at 17:43
• It's all about providing good starting values when there are 7 parameters (a, b, c, d, e, and the error variance) to estimate with only 18 (or 17 legitimate) limited precision data points. Following @Coolwater 's approach gets you good starting values for use in NonlinearModelFit.
– JimB
Commented Jan 28, 2019 at 18:31
• @JimB: I tried you advice, obtaining "NonlinearModelFit::nrjnum: The Jacobian is not a matrix of real numbers at {a,b,c,d,e} = {2.75699,1.66039,2.92883,4.74005,0.70619}.". Commented Jan 28, 2019 at 18:37
• Drop the observation with x=0 and it will work just fine.
– JimB
Commented Jan 28, 2019 at 18:38

You can minimize the cost function by using NMinimize

  ClearAll["Global*"]

data = {{0., 2.61}, {0.1, 2.62}, {0.2, 2.62}, {0.3, 2.62}, {0.4,
2.63}, {0.5, 2.63}, {0.6, 2.74}, {0.7, 2.98}, {0.8, 3.66}, {0.9,
5.04}, {1., 7.52}, {1.1, 10.74}, {1.2, 12.62}, {1.3, 10.17}, {1.4,
5}, {1.5, 2.64}, {1.6, 11.5}, {1.65, 35.4}};

model[x_] := a Cosh[b x^c Sin[d x^e]]
cost = Total[(#2 - model[#1])^2 & @@@ data];
fit = NMinimize[{cost, 1 < a < 3, 1 < b < 3, 1 < c < 3, 1 < d < 3, 1 < e < 3},
{a, b, c, d, e}, Method -> "DifferentialEvolution"]


{0.00112965, {a -> 2.61748, b -> 1.71949, c -> 2.30924, d -> 1.50333, e -> 1.84597}}

Thread[{a, b, c, d, e} = {a, b, c, d, e} /. Last@fit];

Show[ListPlot[data], Plot[model[k], {k, 0, 1.7}]]


You can also use NonlinearModelFit with appropriate constrain on parameters and method choice.

ClearAll["Global*"]

data = {{0., 2.61}, {0.1, 2.62}, {0.2, 2.62}, {0.3, 2.62}, {0.4,
2.63}, {0.5, 2.63}, {0.6, 2.74}, {0.7, 2.98}, {0.8, 3.66}, {0.9,
5.04}, {1., 7.52}, {1.1, 10.74}, {1.2, 12.62}, {1.3, 10.17}, {1.4,
5}, {1.5, 2.64}, {1.6, 11.5}, {1.65, 35.4}};

model[x_] := a Cosh[b x^c Sin[d x^e]]
nlm = NonlinearModelFit[
data, {model[x], 1 < a < 3, 1 < b < 3, 1 < c < 3, 1 < d < 3,
1 < e < 3}, {a, b, c, d, e}, x,
Method -> {"NMinimize", {Method -> "DifferentialEvolution"}}, MaxIterations -> 1000]

nlm // Normal


2.61749 Cosh[1.71949 x^2.30926 Sin[1.50334 x^1.84596]]

nlm["BestFitParameters"]


{a -> 2.61749, b -> 1.71949, c -> 2.30926, d -> 1.50334, e -> 1.84596}

Same picture

• Why was my comment " ClearAll kills the data so the data should be inputted after ClearAll" deleted by an unknown to me person? This is unfair. Commented Jan 29, 2019 at 12:29
• Maybe somebody deleted because I modified the code and now it included the data in the code after ClearAll. Commented Jan 29, 2019 at 12:34

Below err is defined as the residual sum of squares. Using ?NumericQ on one of the arguments prevents FindMinimum from exact differentiation. In that way no problem arises even though the point with x == 0 is included.

err[a_, b_, c_, d_, e_?NumericQ] = Total[(a*Cosh[b*#^c*Sin[d*#^e]] - #2)^2 & @@@ data];

res[x_] = a*Cosh[b*x^c*Sin[d*x^e]] /. Last[FindMinimum[err[a, b, c, d, e], {
{a, 3, 1/100, 25}, {b, 2, 1/100, 25}, {c, 1, 1/100, 25},
{d, 23/10, 1/100, 25}, {e, 2, 1/100, 25}}, Method -> "InteriorPoint"]]

Show[Plot[res[x], {x, 0, 1.65}, PlotRange -> All], ListPlot[data]]


Alternative start values result in a smaller sum of squares:

FindMinimum[err[a, b, c, d, e],
{{a, 2.6174, 1/100, 25}, {b, 1.7195, 1/100, 25}, {c, 2.3092, 1/100, 25},
{d, 1.5033, 1/100, 25}, {e, 1.845, 1/100, 25}}, Method -> "InteriorPoint"]


{0.0011296543, {a -> 2.6174825, b -> 1.7194932, c -> 2.3092448, d -> 1.5033314, e -> 1.845972}}

These are highly likely optimal: When x == 0 the regression formula equals just a and because of the first data point a should presumably be close to 2.61.

If the last factor and Cosh are cancelled from the data values (using the right branch of ArcCosh) we get something we 2 peaks at high x-values:

data2 = Thread[{data[[All, 1]], PadLeft[{-1, -1}, 18, 1] ArcCosh[data[[All, 2]]/2.61]}]
ListPlot[data2]


We can solve for b and c such that the peaks are interpolated leaving only d and e for estimation:

sol = First[Solve[b #1^c Sin[d #1^e] == #2 & @@@ data2[[{-1, -6}]], {b, c}] /. C[1] -> 0 // Chop]


For different values of d and e let's look at the residual sum of squares excluding x == 0, because sol can't be evaluated in that instance

err0[a_, b_, c_, d_, e_?NumericQ] = Total[(a*Cosh[b*#^c*Sin[d*#^e]] - #2)^2 & @@@ Rest[data]];
search = Table[{d, e, Log[If[Im[#] == Im[#2] == 0,
Quiet[Min[err0[2.61, #, #2, d, e], 10^5.]], 10^5.]]} & @@
({b, c} /. sol), {d, 1/50, 7, 1/50}, {e, 1/50, 7, 1/50}] // Catenate;
ListPointPlot3D[search, PlotRange -> All, AxesLabel -> {x, y, z}]


the 2 smallest local minima of which coincide with the previous fits.

• Thank you. Likely a somewhat better fit than your $2.75699 \cosh \left(1.66039 x^{2.92883} \sin \left(4.74005 x^{0.70619}\right)\right)$ is $[a= 2.61748239217892,b= 1.71949328471199,c= 2.30924398627606,d= 1.50333104215966,e= 1.84597270613482]$. Can you kindly compare the results? Commented Jan 28, 2019 at 17:42
• Colleagues, let us call things by their proper names. The answer of @Coolwater is not bad, but is not optimal. Commented Jan 28, 2019 at 18:40
• That's interesting, but for $[a= 2.61748239217892,b= 1.71949328471199,c= 2.30924398627606,d= 1.50333104215966,e= 1.84597270613482]$ the sum of squared residuals equals $0.0000627585700895793$. Hope you feel the difference. Commented Jan 28, 2019 at 19:01
• You can see it at dropbox.com/s/o08c8x0cm0dutcl/deep%20fit.pdf?dl=0 . Commented Jan 28, 2019 at 19:50
• From where "Alternative start values {a,}, {b, 1.7195,}, {c, 2.3092,}, {d, 1.5033}, {e, 1.845}" which lead to {a -> 2.6174825, b -> 1.7194932, c -> 2.3092448, d -> 1.5033314, e -> 1.845972}} are taken? Commented Jan 28, 2019 at 20:17

Just an extended comment that for this particular data/model combination it is difficult to find the values of parameters that minimize the residual sum of squares. I think there are two main issues:

1. The surface is extremely bumpy which can get any iterative search algorithms to get lost. (This is not necessarily a Mathematic, Maple, R, SAS, or *MathCAD" issue.)

2. Some of the parameter estimators are highly correlated with each other. That also causes problems for finding the optimal parameter values.

Below is a contour plot of the sum of squares evaluated at @user64494 's optimal values of $$a$$, $$b$$, and $$c$$ in the neighborhoods of $$d$$ and $$e$$. One can see the extreme bumpiness and the highly correlated nature of the estimators of $$d$$ and $$e$$. The red dot indicates the location of the optimal values of $$d$$ and $$e$$.

• +1. Thank you for your valuable info. I think the increase of the data size to 30 changes nothing here. Commented Jan 28, 2019 at 21:00

NonlinearModelFit works without additional assumptions(Thanks @JimB 's comment) with Method->"NMinimize"

mod = NonlinearModelFit[Rest[data],a*Cosh[b*x^c*Sin[d*x^e]], { a , b , c, d,  e }, x,Method -> "NMinimize" ]
Show[{ListPlot[data], Plot[Normal[mod], {x, 0, 1.65}]}]


• 1. You changed the data by Rest. 2. Your fit $2.76855 \cosh \left(1.65609 x^{2.93246} \sin \left(4.74049 x^{0.705923}\right)\right)$ is far away from optimal ones (see dropbox.com/s/o08c8x0cm0dutcl/deep%20fit.pdf?dl=0 ). Commented Jan 28, 2019 at 20:10
• Yes I omitted the first data point(as @JimB recommended). My plot seems to fit the data quite well , not more. Commented Jan 28, 2019 at 20:16
• One of the aims of my question is to pay attention of Mathematica developers to the fact that NonlinearModelFit faces problems with some model functions. Commented Jan 28, 2019 at 20:21