(* 
 This is a simple example of plotting different combinations of variables where each are stored as a vector David Stuart, January 2013
*)

(* 
 First we need the error bar plotting package 
*)
Needs["ErrorBarPlots`"];

(* 
 Next we define the data for each variable as vectors, which is a bunch of points of (a,b) with an uncertainty for b of db.
*)

a={1,2,3,4,5};
b={0.1,1.95,4.1,6.0,9.1};
db={0.06,0.05,0.04,0.3,1.8};

(* 
 First let's make a simple plot of y vs x by making a table of the data and plotting it 
*)
xy=Table[{a[[i]],b[[i]]},{i,1,Length[a]}];

ListPlot[xy]

(* 
 Now, fit the data to a simple line and then overlay the fitted function and the original data using the Show command.
*)
linearFit = Fit[xy,{1,x}, x];

Show[ListPlot[xy,PlotStyle->Red],Plot[linearFit,{x,0,5}]]

(* 
 But that fit and plot did not include the error bars. We do that now by making a Table including the error bar specification 
*)
xyerr=Table[{{a[[i]],b[[i]]},ErrorBar[db[[i]]]},{i,1,Length[a]}];

(* And plot it using ErrorListPlot *)
ErrorListPlot[xyerr]

(* 
 Do a non-linear fit including the uncertainties. The uncertainties are handled as weights which are 1/error^2. 
*)
fitResult = NonlinearModelFit[xy,intercept+slope*x,{{intercept,1.},{slope,2.}},x, Weights->1/db^2]

fitResult["ParameterTable"]

Show[
Plot[fitResult[x],{x,0,6}, PlotLabel->"Linear fit to b(a)",AxesLabel->{"a","b"},PlotRange->{-2.,10.}],
ErrorListPlot[xyerr]]

(* 
 Print the 68% confidence level for each of the parameters 
*)
fitResult["ParameterConfidenceIntervals",ConfidenceLevel->0.68]

(* 
Now calculate the residuals which are how far each point is from the fitted value.
*)

(* 
First get the simple residuals directly from the fit 
*)

resid = fitResult["FitResiduals"]

(* 
Then get the residuals with uncertainties for easier plotting 
*)
residWithErrors = Table[{{xy[[i,1]],xy[[i,2]] - fitResult[xy[[i,1]]]},ErrorBar[db[[i]]]},{i,1,Length[a]}];

ErrorListPlot[residWithErrors, PlotLabel->"Residuals",PlotRange->{-1,1.}, AxesOrigin->{0,0}, AxesLabel->{"a","b-Fit"}]

(* Print a table of the output and confidence intervals *)
fitResult["SinglePredictionConfidenceIntervalTable", ConfidenceLevel->0.68]