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,

Improve this question
$\endgroup$
7
  • 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, 2016 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, 2016 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, 2016 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, 2016 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, 2016 at 23:08

1 Answer 1

Reset to default
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

Improve this answer
$\endgroup$
5
  • $\begingroup$ Alternatively: weights = Normalize[weights, First]. $\endgroup$
    – J. M.'s lack of A.I. ♦
    Jul 15, 2016 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, 2016 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, 2016 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, 2016 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, 2016 at 0:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.