4
$\begingroup$

I have the following data.

HSI = {{1.`, 4.502201319666759`}, {2.`, 4.498737206541754`}, {3.`, 
   4.4964025063938955`}, {4.`, 4.491168609823489`}, {5.`, 
   4.480225462058295`}, {6.`, 4.472178038829604`}, {7.`, 
   4.4743675246518775`}, {8.`, 4.4582896838826604`}, {9.`, 
   4.45235342111989`}, {10.`, 4.47731636067876`}, {11.`, 
   4.457492363168564`}, {12.`, 4.465951955374022`}, {13.`, 
   4.464606374521092`}, {14.`, 4.45898605964085`}, {15.`, 
   4.457726468220152`}, {16.`, 4.457199643172643`}, {17.`, 
   4.444311998210321`}, {18.`, 4.450761654864603`}, {19.`, 
   4.436634076088381`}, {20.`, 4.433101795058949`}, {21.`, 
   4.449683071507449`}, {22.`, 4.44430904025694`}, {23.`, 
   4.433387852516324`}, {24.`, 4.425938152412656`}, {25.`, 
   4.421420618879485`}, {26.`, 4.415008417671637`}, {27.`, 
   4.409255990949457`}, {28.`, 4.408889403725628`}, {29.`, 
   4.4067336007306634`}, {30.`, 4.389233972225948`}}

I fit it with the following equation.

nlm = NonlinearModelFit[HSI, {A + B (t)^z + (c (t)^z)*(Cos[\[Omega]*Log[t] + \[Phi]]), -1 < c < 1, 0.1 <= z <= 0.9, 4.8 <= \[Omega] <= 13, 0 <= \[Phi] <= 2 \[Pi]}, {A, B, z, c, \[Omega], \[Phi]}, t, MaxIterations -> 10000];

fitfun /. nlm["BestFitParameters"];Show[Plot[nlm[t], {t, 1, 30}], ListPlot[HSI]]

I get the following fit.

enter image description here

How do I get the root mean squared error for this fit?

$\endgroup$
  • 1
    $\begingroup$ Do you mean something like (I am using my own data) $${}$$ data = {{1, 2.9}, {2, 5.0}, {3, 7.1}, {4, 8.9}, {5, 10.9}, {6, 13.1}, {7, 14.9}, {8, 16.8}, {9, 18.5}, {10, 22.2}}; f[x_] = Fit[data, {1, x}, x];RootMeanSquare[data[[All, 2]] - Map[f, data[[All, 1]]]]? $\endgroup$ – Moo Jan 5 at 17:56
  • 3
    $\begingroup$ nlm["EstimatedVariance"]^0.5 is the more appropriate estimate for the root mean square for a regression as the denominator is n - p` where n is the sample size and p is the number of estimated parameters. RootMeanSquare just divides by n. $\endgroup$ – JimB Jan 5 at 22:54
5
$\begingroup$
HSI = {{1., 4.502201319666759}, {2., 4.498737206541754}, {3., 
    4.4964025063938955}, {4., 4.491168609823489}, {5., 
    4.480225462058295}, {6., 4.472178038829604}, {7., 
    4.4743675246518775}, {8., 4.4582896838826604}, {9., 
    4.45235342111989}, {10., 4.47731636067876}, {11., 
    4.457492363168564}, {12., 4.465951955374022}, {13., 
    4.464606374521092}, {14., 4.45898605964085}, {15., 
    4.457726468220152}, {16., 4.457199643172643}, {17., 
    4.444311998210321}, {18., 4.450761654864603}, {19., 
    4.436634076088381}, {20., 4.433101795058949}, {21., 
    4.449683071507449}, {22., 4.44430904025694}, {23., 
    4.433387852516324}, {24., 4.425938152412656}, {25., 
    4.421420618879485}, {26., 4.415008417671637}, {27., 
    4.409255990949457}, {28., 4.408889403725628}, {29., 
    4.4067336007306634}, {30., 4.389233972225948}};

nlm1 = NonlinearModelFit[
   HSI, {A + B t^z + (c t^z)*(Cos[ω*Log[t] + ϕ]), -1 < c < 1, 
    0.1 <= z <= 0.9, 4.8 <= ω <= 13, 0 <= ϕ <= 2 π}, {A, B, z,
     c, ω, ϕ}, t, MaxIterations -> 10000];

Rationalize input and use arbitrary-precision by setting WorkingPrecision

nlm2 = NonlinearModelFit[
   Rationalize[HSI, 
    0], {A + B t^z + (c t^z)*(Cos[ω*Log[t] + ϕ]), -1 < c < 1, 
    1/10 <= z <= 9/10, 24/5 <= ω <= 13, 0 <= ϕ <= 2 π}, {A, B,
     z, c, ω, ϕ}, t, MaxIterations -> 10000, 
   WorkingPrecision -> $MachinePrecision];

EDIT: As suggested by @JimB

nlm3 = NonlinearModelFit[
   Rationalize[HSI, 
    0], {A + B t^z + (c t^z)*(Cos[ω*Log[t] + ϕ]), -1 < c < 1, 
    1/10 <= z <= 9/10, 24/5 <= ω <= 13, 
    0 <= ϕ <= 2 π}, {{A, 45/10}, {B, -5/100}, {z, 3/10}, {c, 
     1}, {ω, 85/10}, {ϕ, 1}}, t, MaxIterations -> 10000, 
    WorkingPrecision -> $MachinePrecision];

Visually comparing the models,

Plot[{nlm1[t], nlm2[t], nlm3[t]}, {t, 1, 30},
 Epilog -> {Red, AbsolutePointSize[4], Point[HSI]},
 PlotStyle -> {Medium, Medium, Thick},
 PlotLegends -> Placed["Expressions", {0.8, 0.6}]]

enter image description here

Comment by @Moo showed how to calculate the RootMeanSquare. Comparing the results,

RootMeanSquare[HSI[[All, 2]] - # /@ HSI[[All, 1]]] & /@ {nlm1, nlm2, 
  nlm3}

(* {0.0120219, 0.0107395, 0.00731085} *)

indicates that nlm3 is the better fit (i.e., smaller RootMeanSquare).

$\endgroup$
  • $\begingroup$ Several of the parameter estimators are highly correlated with each other which makes convergence difficult. If one uses starting values {{a, 45/10}, {b, -5/100}, {z, 3/10}, {c, 1}, {\[Omega], 85/10}, {\[Phi], 1}}, then the root mean square is about 30% smaller than what one gets with nlm2. $\endgroup$ – JimB Jan 6 at 0:19
  • $\begingroup$ Uh, sorry. I changed A to a and B to b. Changing them back to A and B should make the code work. $\endgroup$ – JimB Jan 6 at 1:08
  • $\begingroup$ @JimB - added your input. Thanks. $\endgroup$ – Bob Hanlon Jan 6 at 1:28
  • $\begingroup$ Thank you very much! You helped me beyond what I asked for. $\endgroup$ – Linn Htet Aung Jan 6 at 2:12
7
$\begingroup$

This dataset and model have "issues" because of a few of the parameter estimators being highly correlated with each other. There are 7 parameters to be estimated: A, B, c, z, ω, ϕ, and σ. Given z, ω, ϕ, σ, we can obtain symbolic estimates of A, B, and c as the problem becomes just a simple linear regression. Then we can find the estimates of z, ω, ϕ, and σ.

(* Determine log likelihood *)
data = Rationalize[HSI, 0];
logL = LogLikelihood[NormalDistribution[0, σ], (#[[2]] - (a + b #[[1]]^z + (c #[[1]]^z)*
  Cos[ω*Log[#[[1]]] + ϕ]) &) /@ data];

(* Symbolic solution for a, b, and c given σ, ω, ϕ, and z *)
dlogL = D[logL, {{a, b, c}}];
sol = Solve[dlogL == 0, {a, b, c}];

(* Numerically solve for maximum likelihood estimates of σ, ω, ϕ, and z *)
logLs = logL /. sol;
mleωϕzσ = FindMaximum[logLs, {{σ, 7/1000}, {ω, 10}, {ϕ, 0}, {z, 8/10}}]
(* {105.054, {\[Sigma] -> 0.00729368, \[Omega] -> 10.3986, \[Phi] -> -0.580638, z -> 0.924522}} *)

(* Determine maximum likelihood estimates for a, b, and c *)
mleabc = sol /. mleωϕzσ[[2]]
(* {{a -> 4.50173, b -> -0.00402987, c -> 0.000382872}} *)

(* Now get estimate of parameter covariance and correlation matrix *)    
hessian = D[logL, {{a, b, c, z, ω, ϕ, σ}, 2}];
cov = -Inverse[hessian /. mleωϕzσ[[2]] /.  mleabc[[1]]]
cor = Table[cov[[i, j]]/(cov[[i, i]] cov[[j, j]])^0.5, {i, 1, 7}, {j, 1, 7}]

Now check of the root mean square:

RootMeanSquare[(#[[2]] - (a + b #[[1]]^z + (c #[[1]]^z)*Cos[ω*Log[#[[1]]] + ϕ]) /. 
  mleωϕzσ[[2]] /. mleabc[[1]] &) /@ Rationalize[HSI]]
(* *)

This matches with the maximum likelihood estimate of σ and is even a bit smaller than in other answers (0.00729368).

$\endgroup$
  • 1
    $\begingroup$ Why the commas in cov[[i, i]], cov[[j, j]] in the definition of cor? Is that a typo? $\endgroup$ – chris Jan 6 at 9:51
  • $\begingroup$ @chris Yep, a typo. Sorry about that. Thanks for catching that. I'll edit it. $\endgroup$ – JimB Jan 6 at 15:09
6
$\begingroup$

You can also use RootMeanSquare on the property "FitResiduals":

RootMeanSquare[#["FitResiduals"]] & /@ {nlm1, nlm2, nlm3}

{0.0120112, 0.034652979361999, 0.007328510886305}

where nlm1, nlm2 and nlm3 are from Bob Hanlon's 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.