Message Boards Message Boards

GROUPS:

NonlinearModelFit - syntax and interpreting the fit

Posted 9 years ago
15512 Views
|
15 Replies
|
7 Total Likes
|

Hi,

I have a set of data I am trying to fit into a non-linear model fit.

The equation is

y = 1 - ( 1 / ( 1 + 2 J x ) ^3 )

I have values for y and x, and want to calculate the value of J.

I tried the following code with data as my x/y values, using an initial guess of 0.01 for b:

nlm = NonlinearModelFit[  data, [((1/(1 + (2 b x ))^3))], {b, 0.01}, x] 

But I get an error. Can you help/guide me?

Thanks in advance.

POSTED BY: syed islam
15 Replies
Posted 9 years ago

The goodness-of-fit is many times in the eye of the beholder. As a statistician I see that the first five values have positive residuals and the last three have negative residuals. That suggests that the model does not fit well (either the proposed curve form doesn't fit or the assumption of additive independent random errors is not supported or both). Also, R-squared is not so meaningful when describing the strength of the relationship of curves forced through the origin. In addition, I’d argue that the point {0,0} is likely NOT a data point as the curve form has to go through the origin.

A better measure of goodness is the root mean square error which is in the same units as the dependent variable. Given the undesirable residuals you might consider a slightly more complex model where the exponent "3" is replaced with a parameter.

data = {{1800, 0.258738807}, {3600, 0.362822881}, {5400, 
    0.473489839}, {7200, 0.54841429}, {9000, 0.614915391}, {10800, 
    0.638559079}, {14400, 0.68690743}, {18000, 0.730710491}};

(* Fit the model using the data *)
nlm = NonlinearModelFit[data, 1 - ((1/(1 + (b (x/1000)))^k)), {b, k}, 
  x]

(* Parameter estimates *)
nlm["BestFitParameters"]

(* Root mean square error *)
nlm["EstimatedVariance"]^0.5

(* Show the fit *)
Show[Plot[nlm[x], {x, 0, 18000}, PlotRange -> All], ListPlot[data]]

with output

{b -> 0.225117, k -> 0.818872}
0.0140814

Graph of data and model fit

I've also explicitly divided the independent variable by 1,000 as standardizing independent variables (either by subtracting the mean and dividing by the standard deviation or subtracting the minimum and dividing by the range) many times fixes numerical instabilities.

POSTED BY: Jim Baldwin
Posted 9 years ago

Thanks David Reiss and David Keith. I originally tried this using Excel solver and got a similar value, though my supervisor was not happy with this and recommended mathematica, which I am new too.

As the

nlm["AdjustedRSquared"]

is 0.992 I guess that means my data is a good fit for the value of b.

POSTED BY: syed islam
Posted 9 years ago

Hi David. Are you sure? I used b->b/10000, so in your Out[85] there is still an expression which contains b. And that is the expression defining the model, with parameter b.

POSTED BY: David Keith

Whoops! My bad! I need more coffee. Care to join me?

Yes your approach pushes the b into the region where whatever the Automatic Method is using will properly find it rather than send the search in the wrong direction.

POSTED BY: David Reiss
Posted 9 years ago

Coffee sounds great! But we'd need to meet somewhere in the middle, like maybe the Denver airport. ;-}

POSTED BY: David Keith

I think Oregon is the better choice...

POSTED BY: David Reiss
Posted 9 years ago

April through October. In the Columbia River Gorge, or the Mt Hood National Forest.

Winter is Mathematica's season. And Christmas parties!

POSTED BY: David Keith
Posted 9 years ago

In the problem as stated, the best fit for b appears to be around 0.00008. For some reason NonlinearModelFit does not find it. If we rewrite the model slightly, so a good fit for b is about 0.8, it does well.

nlm = NonlinearModelFit[data, 
   1 - ((1/(1 + (2 b x/4.0952977))^3)) /. b -> b/10000, {b}, x];

Show[ListPlot[data], Plot[nlm[x], {x, 0, 18000}]]

enter image description here

POSTED BY: David Keith

Yes, but David your example is not fitting the function since it is specifying the b inside and so the actual function it is trying to fit,

In[85]:= 1 - ((1/(1 + (2 b x/4.0952977))^3)) /. b -> b/10000

Out[85]= 1 - 1/(1 + 0.0000488365 b x)^3

has no parameters...

POSTED BY: David Reiss

Yes, it is a lousy fit ;-)

In[45]:= nlm["AdjustedRSquared"]

Out[45]= 0.156392

When doing NonlinearModelFit the default Method option is Method->Automatic and this tells Mathematica to choose a method that it thinks is the right one. Apparently it is failing miserably in this case. So the next step is to try a specific Method. The documentation tells us that there are the following choices:

"ConjugateGradient", "Gradient", "LevenbergMarquardt", "Newton", "NMinimize", "QuasiNewton"

So best to explore them with a Manipulate:

Manipulate[
 data = {{0, 0}, {1800, 0.258738807}, {3600, 0.362822881}, {5400, 
    0.473489839}, {7200, 0.54841429}, {9000, 0.614915391}, {10800, 
    0.638559079}, {14400, 0.68690743}, {18000, 0.730710491}};
 nlm = NonlinearModelFit[data, 1 - ((1/(1 + ( b x))^3)), {b}, x, 
   Method -> method];
 Column[{
   Show[Plot[nlm[x], {x, 0, 18000}, PlotRange -> All], ListPlot[data]],
   nlm["AdjustedRSquared"]}
  ],
 {method, {"ConjugateGradient", "Gradient", "LevenbergMarquardt", 
   "Newton", "NMinimize", "QuasiNewton"}}
 ]

Note that I changed your function to

1 - ((1/(1 + ( b x))^3))

as there are no reasons to include all the extra constants--they cal all be subsumed into your b.

POSTED BY: David Reiss
Posted 9 years ago

David,

The plot function works, but the graph is just showing y = 1, and not an nonlinear graph.

Plot[nlm[x], {x, 0, 18000}, PlotRange -> All]

any ideas?

POSTED BY: syed islam

The <<>> in the formatted output is a shorthand to indicate that there is an expression there that is too large to fit in the formatted version.

If you execute

InputForm[nlm]

you will see that everything you expect (and more) is there.

For plotting your expression is correct. Note that Mathematica chooses not to display all of the plot's y-range since the "interesting" part of the function would then note be visible. To see the entire y-range you can set PlotRange->All as in

Plot[nlm[x], {x, 0, 18000}, PlotRange -> All]
POSTED BY: David Reiss
Posted 9 years ago

Hi David,

Thanks. That worked. The actual data is:

data = {{0, 0}, {1800, 0.258738807}, {3600, 0.362822881}, {5400, 0.473489839}, {7200, 0.54841429}, {9000, 0.614915391}, {10800, 0.638559079}, {14400, 0.68690743}, {18000, 0.730710491}};

I wrote the equation wrong. The actual equation is:

nlm = NonlinearModelFit[data, 1 - ((1/(1 + (2 b x/4.0952977))^3)), {b}, x].

Which gives an output of:

1 - 1/(1 + <<20>> x)^3

What do the <<>> mean?

For a plot of the graph is the following code correct:

Plot[nlm[x], {x, 0, 18000}]

Thanks,

POSTED BY: syed islam

What error do you get? The form of the NonlinearModelFit that you wrote has two problems. First is that the square brackets are incorrect--they shouldn't be there. Second is that the parameter specification should be {b} rather than {b, 0.01}.

So, with artificial data:

In[7]:= data={{0,1},{1,0},{3,2},{5,4}, {6,4},{7,5}};
In[8]:= NonlinearModelFit[ data, ((1/(1 + (2 b x ))^3)), {b}, x]
Out[8]= FittedModel[1/(1-0.0612015 x)^3]
POSTED BY: David Reiss
Posted 4 years ago

I don't understand why {b, 0.01} is incorrect for an initial guess for b?

POSTED BY: James McBride
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract