3
$\begingroup$

I am trying to plot the best-fit curve (with known functional form but free parameters) passing through the following set of discrete data points which would take into account the size of error bars using Least Squares Fitting Procedure. I would like to find the best parameters along with their uncertainties as follows:

Needs["ErrorBarLogPlots`"]
Needs["ErrorBarPlots`"]
h=0.7;
Data = {{{9.6 + Log10[h], Log[10]*10^(-2.292)}, ErrorBar[Log[10]*10^(-2.292)*{10^(-0.063) - 1, 10^0.072 - 1}]}, {{9.7 + Log10[h], Log[10]*10^(-2.347)}, ErrorBar[Log[10]*10^(-2.347)*{10^(-0.029) - 1, 10^0.031 - 1}]}, {{9.8 + Log10[h], Log[10]*10^(-2.289)}, ErrorBar[Log[10]*10^(-2.289)*{10^(-0.028) - 1, 10^0.030 - 1}]}, {{9.9 + Log10[h], Log[10]*10^(-2.308)}, ErrorBar[Log[10]*10^(-2.308)*{10^(-0.036) - 1, 10^0.040 - 1}]}, {{10 + Log10[h], Log[10]*10^(-2.325)}, ErrorBar[Log[10]*10^(-2.325)*{10^(-0.027) - 1, 10^0.028 - 1}]}, {{10.1 + Log10[h], Log[10]*10^(-2.253)}, ErrorBar[Log[10]*10^(-2.253)*{10^(-0.073) - 1, 10^0.087 - 1}]}, {{10.2 + Log10[h], Log[10]*10^(-2.342)}, ErrorBar[Log[10]*10^(-2.342)*{10^(-0.028) - 1, 10^0.030 - 1}]}, {{10.3 + Log10[h], Log[10]*10^(-2.372)}, ErrorBar[Log[10]*10^(-2.372)*{10^(-0.025) - 1, 10^0.027 - 1}]}, {{10.4 + Log10[h], Log[10]*10^(\[Minus]2.327)}, ErrorBar[Log[10]*10^(-2.327)*{10^(-0.033) - 1, 10^0.036 - 1}]}, {{10.5 + Log10[h], Log[10]*10^(-2.332)}, ErrorBar[Log[10]*10^(-2.332)*{10^(-0.028) - 1, 10^0.030 - 1}]}, {{10.6 + Log10[h], Log[10]*10^(-2.384)}, ErrorBar[Log[10]*10^(-2.384)*{10^(-0.026) - 1, 10^0.028 - 1}]}, {{10.7 + Log10[h], Log[10]*10^(-2.360)}, ErrorBar[Log[10]*10^(-2.360)*{10^(-0.031) - 1, 10^0.033 - 1}]}, {{10.8 + Log10[h], Log[10]*10^(-2.493)}, ErrorBar[Log[10]*10^(-2.493)*{10^(-0.028) - 1, 10^0.029 - 1}]}, {{10.9 + Log10[h], Log[10]*10^(-2.644)}, ErrorBar[Log[10]*10^(-2.644)*{10^(-0.036) - 1, 10^0.039 - 1}]}, {{11 + Log10[h], Log[10]*10^(-2.734)}, ErrorBar[Log[10]*10^(-2.734)*{10^(-0.036) - 1, 10^0.039 - 1}]}, {{11.1 + Log10[h], Log[10]*10^(-2.978)}, ErrorBar[Log[10]*10^(-2.978)*{10^(-0.047) - 1, 10^0.052 - 1}]}, {{11.2 + Log10[h], Log[10]*10^(-3.114)}, ErrorBar[Log[10]*10^(-3.114)*{10^(-0.057) - 1, 10^0.066 - 1}]}, {{11.3 + Log10[h], Log[10]*10^(-3.46)}, ErrorBar[Log[10]*10^(-3.46)*{10^(-0.083) - 1, 10^0.10 - 1}]}, {{11.4 + Log10[h], Log[10]*10^(-3.67)}, ErrorBar[Log[10]*10^(-3.67)*{10^(-0.10) - 1, 10^0.10 - 1}]}, {{11.5 + Log10[h], Log[10]*10^(-4.12)}, ErrorBar[Log[10]*10^(-4.12)*{10^(-0.20) - 1, 10^0.30 - 1}]}, {{11.6 + Log10[h], Log[10]*10^(-4.35)}, ErrorBar[Log[10]*10^(-4.35)*{10^(-0.20) - 1, 10^0.40 - 1}]}, {{11.7 + Log10[h], Log[10]*10^(-5.09)}, ErrorBar[Log[10]*10^(-5.09)*{10^(-0.40) - 1, 10^1.00 - 1}]}, {{11.8 + Log10[h], Log[10]*10^(-5.05)}, ErrorBar[Log[10]*10^(-5.05)*{10^(-0.40) - 1, 10^1.00 - 1}]}};
model = 
  Log[10]*Exp[-10^(x - Log10[h] - a)]*10^b*(10^(x - Log10[h] - a))^(c + 1);
curve = 
 NonlinearModelFit[Data, model, {a,b,c}, x]

And, then to plot both raw data and best-fit on the same plot using Show command as follows:

Show[ErrorListLogPlot[Data, 
  PlotRange -> {{9.25, 11.95}, {5*10^-6, 0.2}}, 
  PlotStyle -> {Red, Thick}, Joined -> False, Frame -> True, 
  FrameLabel -> {Style["X", FontSize -> 24], Style["Y", FontSize -> 24]}, 
  FrameTicksStyle -> Directive[FontSize -> 24]], 
  LogPlot[curve[x], {x, 9.25, 11.95}, PlotRange -> {5*10^-6, 0.2}, 
  PlotStyle -> Black]]

enter image description here However, the second step (plotting) will be possible only once the free parameters and their uncertainties are found. Your help is greatly appreciated,

$\endgroup$
  • 1
    $\begingroup$ Data is not a {x,y} data. Use NonlinearModelFit[Data[[All,1]], model, {a, b, c, d, f}, x] $\endgroup$ – Sumit Jul 14 '16 at 20:55
  • 1
    $\begingroup$ Your model does not seem to easily reproduce the trend in your experimental data. You will need to provide much better starting values for your parameters. Also, you can use 1/errors as weights in your fit (see the Weights option), but not in the format in which you have them for plotting. I'd suggest that you ignore the weigthing issue first, and get some decent starting values (if a, b, c etc have physical meaning, then you should be able to provide some better estimate; if they don't, then simplify your model). $\endgroup$ – MarcoB Jul 14 '16 at 21:47
  • $\begingroup$ @Sumit, thank for the clarification. But I am not getting the uncertainties of the central values. Is there any other approach so I can get those? $\endgroup$ – Benjamin Jul 14 '16 at 21:59
  • $\begingroup$ @MarcoB, thanks for the feedback. You are right. Now, I have a better model. I edited it. $\endgroup$ – Benjamin Jul 14 '16 at 22:04
  • 2
    $\begingroup$ Implement @MarcoB 's suggestion of better starting values by replacing {a,b,c} with {{a,10},b,c}. $\endgroup$ – JimB Jul 14 '16 at 23:08
4
$\begingroup$

Step 1 - Find reasonable starting values

Typically one uses Manipulate in order to get values that are in the ball park.

Using your data and model create measuredData to be the {x,y} pairs extracted from the data.

h = 0.7;

model = Log[10]*Exp[-10^(x - Log10[h] - a)]*10^b*(10^(x - Log10[h] - a))^(c + 1)

(* 10^b (10^(0.154902 - a + x))^(1 + c) E^-10^(0.154902 - a + x) Log[10] *)

measuredData = data[[All, 1]]

Create a Manipulate and see if it is possible to locate reasonable starting values.

Manipulate[
 Show[
  LogPlot[
   Log[10]*Exp[-10^(x - Log10[h] - a)]*10^
     b*(10^(x - Log10[h] - a))^(c + 1), {x, 9.4, 11.7}, 
   PlotStyle -> Black],
  ListLogPlot[measuredData, PlotStyle -> Red],
  PlotRange -> {Automatic, All}
  ],

 {{a, 11}, 5, 15, Appearance -> "Open"},
 {{b, -2}, -3, 3, Appearance -> "Open"},
 {{c, -1}, -3, 3, Appearance -> "Open"}
 ]

Mathematica graphics

Step 2 - NonlinearModelFit with weights

Create the weights assuming that the uncertainty for each data point equals the difference between the positive and negative error divided by the measured value, this quantity squared.

Mathematica graphics

weights = Map[(measuredData[[#,2]]/(data[[#,2,1,2]] - data[[#,2,1,1]]))^2 &,
               Range[Length@data]]

Normalize the weights:

weights = Normalize[weights]

(* {0.0531704, 0.27532, 0.294666, 0.170646, 0.328493, \
0.0372979, 0.294666, 0.366696, 0.207596, 0.294666, 0.340004, \
0.241928, 0.305815, 0.175643, 0.175643, 0.100157, 0.0641387, \
0.0282171, 0.0244968, 0.00284079, 0.00149458, 0.0000573521, \
0.0000573521} *)

Now run NonlinearModelFit.

weightedCurve =  NonlinearModelFit[measuredData, model,
     {{a, 11}, {b, -2}, {c, -1}}, x, Weights -> weights]

weightedCurve["BestFitParameters"]

(* {a -> 10.8517, b -> -2.10571, c -> -0.830111} *)

weightedCurve["ParameterTable"]

Mathematica graphics

I am unable to locate the ErrorBarLogPlots package. Below is a plot without the error bars.

Show[
 LogPlot[weightedCurve[x], {x, 9.4, 11.7}, PlotStyle -> Black],
 ListLogPlot[measuredData, PlotStyle -> Red],
 PlotRange -> {Automatic, All}
 ]

Mathematica graphics

$\endgroup$
  • $\begingroup$ Alternatively: weights = Normalize[weights, First]. $\endgroup$ – J. M. is away Jul 15 '16 at 1:36
  • $\begingroup$ Dear Jack, thanks for this. I am working to see if this will work for me. I am receiving the following error message for Step 1 without any fit drawn: "LogPlot::exclul: {Im[10^(0.154902 -par1+x)]-0,Im[10^(0.154902 -par1+x)]-0} must be a list of equalities or real-valued functions. >>" $\endgroup$ – Benjamin Jul 15 '16 at 22:37
  • $\begingroup$ For the second method, I am receiving again nonsense answer for the weights command as follows: {72773.9 ErrorBar[{-0.00158728, 0.00211964}]^2, 485450. ErrorBar[{-0.000668975, 0.000766276}]^2...etc and subsequently the next lines are producing lots of error messages. What am I doing incorrectly??? $\endgroup$ – Benjamin Jul 16 '16 at 0:08
  • $\begingroup$ For the second step, I am entering this line first which doesn't seem to act accordingly: weights = Map[(Data[[#,2]]/(Data[[#, 2, 1, 2]] - Data[[#, 2, 1, 1]]))^2 &, Range[Length@Data]] $\endgroup$ – Benjamin Jul 16 '16 at 0:28
  • 1
    $\begingroup$ There is a measuredData constructed from the data that you are missing. measuredData = data[[All, 1]] $\endgroup$ – Jack LaVigne Jul 16 '16 at 0:42

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.