3
$\begingroup$

I have two points stored in data:

data = {{0.0166667, 2.86927*10^-12}, {0.0333333, 1.12725*10^-11}};

data[[All,1]] contains the time t.

data[[All,2]] contains the first 2 points of the so called mean squared displacement.

To calculate the particle mass m the 2 points should be fitted with the following function:

k * T / m * t^2 (* fitting function *)

k = 1.3806488*10^-23;
T = 300;
(* m = particle mass in kg, is the fitting parameter *)

Manually I get a good fit if I take m = 4.1 * 10^-13kg.

dataPlot = 
 ListLinePlot[data, PlotStyle -> {Blue}, 
  Epilog -> {PointSize[Medium], Point[data]}];
manualFitPlot = 
  Plot[k*300/(4.1*10^-13)*t^2, {t, 1/60, 2/60}, PlotStyle -> {Green}];
Show[dataPlot, manualFitPlot, PlotRange -> All]

enter image description here

When I try to use NonlinearModelFit a wrong result for the mass is obtained:

nlm = NonlinearModelFit[data, k*300/m*t^2, {m}, t];
nlm["ParameterTable"]

$\small \begin{array}{l|llll} \text{} & \text{Estimate} & \text{Standard Error} & \text{t-Statistic} & \text{P-Value} \\ \hline m & 0.422029 & 4.36728\times 10^{11} & 9.66343 \times 10^{-13} & 1 \\ \end{array}$

I tried different initial values for m and also rescaled the small values by multiplying with (*10^10). It did not help.

How can I solve this?

$\endgroup$
  • $\begingroup$ Are you really attempting to fit two parameters (m and the variance about the line) with just two data points? $\endgroup$ – JimB Apr 20 '16 at 16:18
3
$\begingroup$

I wouldn't be so quick to rule out rescaling your data or working in reduced units, since this works fine:

data = {{0.0166667, 2.86927*10^-12}, {0.0333333, 1.12725*10^-11}};
k = 1.3806488*10^-23;
T = 300;
nlm = NonlinearModelFit[data, k*300/(m*10^-13)*t^2, {m}, t, 
   Method -> "NMinimize"];

(* Gives m = 4.07829 *)

I think the problem really lies with the numbers being so small that the default gradient-based fitting methods evaluate the gradient as zero in machine precision.

So NMinimize works, as does PrincipalAxis, since both involve derivative-free methods. You can look into it further in the documentation.


Note that you can also do it without scaling, but specifying a starting value (this time using the PrincipalAxis method, which takes two starting values).

nlm = NonlinearModelFit[data, k*300/m*t^2, {{m, 10^-14, 5*10^-13}}, t, 
         EvaluationMonitor :> Print["x = ", m], 
         Method -> "PrincipalAxis"];
(* Gives m = 4.09157*10^-13 *)
(* Also outputs the steps *)
$\endgroup$
  • $\begingroup$ I made a mistake, I rescaled only data[[All,2]] ... thanks a lot $\endgroup$ – mrz Apr 20 '16 at 15:48
  • $\begingroup$ @mrz you don't need to scale - using a gradient-free method will also work - see my edit, but remember you are using machine-precision numbers. $\endgroup$ – dr.blochwave Apr 20 '16 at 15:50
  • $\begingroup$ So, which fitting value is the better one? $\endgroup$ – mrz Apr 20 '16 at 15:54
  • $\begingroup$ That's up to you to decide. NMinimize probably gives the better results with the Nelder-Mead algorithm than PrincipalAxis but I don't know for sure - the crucial point really is derivative-free methods. $\endgroup$ – dr.blochwave Apr 20 '16 at 16:03
  • $\begingroup$ @mrz Probably stick with NMinimize. $\endgroup$ – dr.blochwave Apr 20 '16 at 16:09
1
$\begingroup$

Just a quick glance, you can try this though. I didn't have time to plot the residuals (or take a look at them for that matter...no pun intended) but it could be worth taking a look at. The function seems to predict your data points with relative accuracy. Function

$\endgroup$
1
$\begingroup$

This model is linear in 1/m, hence use LinearModelFit

lm = LinearModelFit[data, {k*300 t^2}, t, IncludeConstantBasis -> False]

Then just invert the result:

1/lm["BestFitParameters"][[1]]

that gives

4.07829*10^-13

If you give different weights to each point, you will find another result. The "best" result will be obtained with proper weighting.

$\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.