# Knowing the functional form, how to find best-fit parameters and their uncertainties for a set of raw data in Mathematica 8?

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]]


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

• Data is not a {x,y} data. Use NonlinearModelFit[Data[[All,1]], model, {a, b, c, d, f}, x] Commented Jul 14, 2016 at 20:55
• 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). Commented Jul 14, 2016 at 21:47
• @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? Commented Jul 14, 2016 at 21:59
• @MarcoB, thanks for the feedback. You are right. Now, I have a better model. I edited it. Commented Jul 14, 2016 at 22:04
• Implement @MarcoB 's suggestion of better starting values by replacing {a,b,c} with {{a,10},b,c}.
– JimB
Commented Jul 14, 2016 at 23:08

## 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"}
]


## 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.

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"]


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}
]


• Alternatively: weights = Normalize[weights, First]. Commented Jul 15, 2016 at 1:36
• 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. >>" Commented Jul 15, 2016 at 22:37
• 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??? Commented Jul 16, 2016 at 0:08
• 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]] Commented Jul 16, 2016 at 0:28
• There is a measuredData constructed from the data that you are missing. measuredData = data[[All, 1]] Commented Jul 16, 2016 at 0:42