I. Example of Weighted Linear Fit, with check of uncertainties in slope and intercept via Monte Carlo calculation.
II. Use of fitting in calibration: Uncertainties in results.
III. Naive method of calibration

<<"ErrorBarPlots`"

I. Example of Weighted Linear Fit, with check of uncertainties via Monte Carlo calculation.

Function & Data Set

Here is a definition of a linear function.  I define it in this seemingly clumsy way so that the techniques can 
be generalized to any function of a given set of parameters --- it could be defined in a more straightforward way.

a = {1,2};
f[x_,a_]:=a[[1]]+a[[2]]*x

data={{0,0.8214,0.1},{1,2.8471,0.3},{2,4.8520,0.5},{3,7.5347,0.7},{4,10.2464,0.9},{5,10.2707,1.1},{6,12.8011,1.3},{7,13.7108,1.5},{8,17.8501,1.7},{9,15.3667,1.9},{10,19.3933,2.1}};

MatrixForm[data]

Separate the list "data" into a list that has x-y pairs only ("xydata"), and a list that has the weights for each point ("w").  This is convenient for the Mathematica syntax in the function LinearModelFit below.

xydata = Table[{data[[i,1]],data[[i,2]]},{i,1,Length[data]}];
w = Table[1/data[[i,3]]^2,{i,1,Length[data]}];

dataplot=ErrorListPlot[{Table[{{data[[i,1]],data[[i,2]]},ErrorBar[data[[i,3]]]},{i,1,Length[data]}]},PlotStyle->PointSize[0.018],PlotRange->All,AxesLabel->{x,y},BaseStyle->{FontFamily->"Times",FontSize->16}]

Define weighted χ2 and unweighted χ2 functions

Define functions that calculate χ2.  The function χ2 includes weights.

χ2[data_,p_] :=Sum[(f[data[[i,1]],p]-data[[i,2]])^2/data[[i,3]]^2,{i,1,Length[data]}]

Minimize weighted χ2 using Mathematica LinearModelFit Function

Clear[x];
lm = LinearModelFit[xydata,{1,x},x,Weights->w]

You can get other statistical information about the fit.  Note that the CovarianceMatrix of Mathematica is the curvature matrix, α,  TIMES  the reduced χ2 of the original data set.. (See my notes and error_test.nb)

lm["BestFitParameters"]

lm["ParameterErrors"]

lm["ParameterTable"]

MatrixForm[lm["CovarianceMatrix"]]

Note that gnuplot returns the CorrelationMatrix .  This is related to CovarianceMatrix  as follows:
              (CorrelationMatrix)_{ij} = (CovarianceMatrix}_{ij}/(Sqrt[(CovarianceMatrix)_{ii}] Sqrt[(CovarianceMatrix)_{jj}])
See Bevington Eq. (7.14)

To use information returned by Mathematica  I give appropriate names to the parameters:

{b,m} = lm["BestFitParameters"]
χ2min =χ2[data,{b,m}]
rχ2 = χ2min/(Length[data]-2)

Mathematica will also return the minimum value of the reduced chi squared :

lm["EstimatedVariance"]

Alternatively, Minimize weighted χ2 by explicit minimization of χ2

a = {bb,mm};      (* Reset slope and intercept as variables to be adjusted.          *)
χ2[data,a];       (* Remove semicolon to see explicit form for χ2                   *)

Find minimum of χ2.  Note that because χ2 is a quadratic in the paramaters b and m it is possible to get explicit formulas for b and m.  See, for example, Taylor Chapter 8.  (Taylor treats case in which uncertainties are all equal.)  In this case you can use the Solve function; when χ2 isn't such a simple function of the parameters NSolve may be necessary.

best =Solve[{D[χ2[data,a],bb]==0,D[χ2[data,a],mm]==0},{bb,mm}]

b = bb/.best[[1]]

m = mm/. best[[1]]

Graphs and residuals

funplot=Plot[lm[x],{x,0,10}];

Show[dataplot,funplot]

Plot residuals:

ErrorListPlot[{Table[{{data[[i,1]],f[data[[i,1]],{b,m}]-data[[i,2]]},ErrorBar[data[[i,3]]]},{i,1,Length[data]}]},PlotRange->{{-1,11},All},PlotStyle->PointSize[0.018],AxesLabel->{x,"Residuals"},BaseStyle->{FontFamily->"Times",FontSize->16}]

Plot Normalized Residuals :

ErrorListPlot[{Table[{{data[[i,1]],(f[data[[i,1]],{b,m}]-data[[i,2]])/data[[i,3]]},ErrorBar[1.]},{i,1,Length[data]}]},PlotRange->{{-1,11},All},PlotStyle->PointSize[0.018],AxesLabel->{x,"Normalized Residuals"},BaseStyle->{FontFamily->"Times",FontSize->16}]

Generate hypothetical data sets to find variation in best fit slope and intercept

Define function that will be used to sample y - values from a normal distribution.

pickY[x_,σ_]:= Random[NormalDistribution[x,σ]]

Assume fitted model is correct and assume errors are stated correctly.  Synthesize a hypothetical data set consistent with these assumptions that is statistically equivalent to the original data set:

hdata = Table[{x = data[[i,1]],pickY[f[x,{b,m}],data[[i,3]]],data[[i,3]]},{i,1,Length[data]}];
MatrixForm[hdata]

npt=1000;
params=Table[{0,0},{i,1,npt}];
For[i=1,i≤npt,++i,
{hdata = Table[{x = data[[i,1]],pickY[f[x,{b,m}],data[[i,3]]],data[[i,3]]},{i,1,Length[data]}],
hxydata=Table[{hdata[[i,1]],hdata[[i,2]]},{i,1,Length[data]}],hw=Table[1/hdata[[i,3]]^2,{i,1,Length[data]}],
Clear[x],
hlm=LinearModelFit[hxydata,{1,x},x,Weights->hw],
params[[i]]= hlm["BestFitParameters"]}]

paramplot=ListPlot[params]

Find average of b and m from all of the HYPOTHETICAL data sets.  The average should be the m and b given by the FIT to data, not b=1, m=2.  (Remeber, hdata generated from original fits; we're only learning about reasonable spread here.)

The actual spreads here can be different from predicted Standard Errors.  Predicted SEs include a factor of Sqrt[χ2/ndf], which reflects "quality" of the data.  If reduced χ2 =1, agreement should be very good.

bbar = Mean[Table[params[[i,1]],{i,1,Length[params]}]]
bdev = StandardDeviation[Table[params[[i,1]],{i,1,Length[params]}]]
mbar = Mean[Table[params[[i,2]],{i,1,Length[params]}]]
mdev = StandardDeviation[Table[params[[i,2]],{i,1,Length[params]}]]

Histogram of slopes from hypothetical data sets:

Histogram[Table[params[[i,2]],{i,1,Length[params]}]]

Examine χ2 with near minimum - Contour Plots

Plot contours of Δχ2 = χ2 - χ2min.   (NOTE: Not contours of reduced Δχ2.)

?ContourPlot

chicontour=ContourPlot[χ2[data,{a1,a2}]-χ2[data,{b,m}],{a1,b-0.25,b+0.25},{a2,m-0.2,m+0.2},Contours->{0,1,2,3,4,5,6},ContourShading->False,FrameLabel->{"b","m"}]

Show[chicontour,paramplot]

Note that another way to find error in slope is to look at range of a2 values defined by the projection of the Δχ2 =1 contour onto the a2 axis.  Similarly, the error in the intercept is the range of a1 values defined by the projection of the contour on the a1 axis.

The joint probability for (a1,a2) can be investigated by looking at the number of points entirely contained within a Δχ2 contour.    About 67% of the pairs (a1,a2) should fall inside the χ2 = 2.3 contour - see p. 555 of Numerical Recipes for discussion.                        
Check this:

nmb = 0;
For[i=1,i≤Length[params],++i,If[χ2[data,{params[[i,1]],params[[i,2]]}]-χ2[data,{b,m}]<2.30,++nmb]]
nmb/Length[params]//N

II. Use of fitting in calibration: Uncertainties in results.

Method A: Generation of multiple statisticallly equivalent data sets

Redo weighted fitting:

best= LinearModelFit[xydata,{1,x},x,Weights->w]

{b,m} = best["BestFitParameters"]
σ2 = best["CovarianceMatrix"]

Assume that you measure some value of y,  and call this value y0, with uncertainty Δy0.

y0= 20;
Δy0 = 2;

Simulate experiment many times:

npt= 1000;
params = Table[{0,0,0},{i,1,npt}];
For[i=1,i≤npt,++i,{hdata = Table[{x =data[[i,1]],pickY[f[x,{b,m}],data[[i,3]]],data[[i,3]]},{i,1,Length[data]}],
hxydata = Table[{hdata[[i,1]],hdata[[i,2]]},{i,1,Length[data]}],
hY = Random[NormalDistribution[y0,Δy0]];
hw = Table[1/hdata[[i,3]]^2,{i,1,Length[data]}],
Clear[x],
hbest= LinearModelFit[hxydata,{1,x},x,Weights->hw] ,
params[[i,1]]= hbest["BestFitParameters"][[1]],
params[[i,2]]=hbest["BestFitParameters"][[2]],
params[[i,3]] = (hY-hbest["BestFitParameters"][[1]])/hbest["BestFitParameters"][[2]]}]

bbar = Mean[Table[params[[i,1]],{i,1,Length[params]}]]
bdev = StandardDeviation[Table[params[[i,1]],{i,1,Length[params]}]]
mbar = Mean[Table[params[[i,2]],{i,1,Length[params]}]]
mdev = StandardDeviation[Table[params[[i,2]],{i,1,Length[params]}]]

x0bar = Mean[Table[params[[i,3]],{i,1,Length[params]}]]      (* Should agree with following line *)
x0best = (y0-b)/m
x0dev = StandardDeviation[Table[params[[i,3]],{i,1,Length[params]}]]

FINAL ANSWER: For measurement of y0 = 20 ±  2, we infer that x0 = 9.7 ± 1.1

If you repeat the above procedures for y0 = 3 ± 0.3, you will find that x0 =1.09 ± 0.16
Notice that the uncertainties are very different!

Method B: Using the covariance matrix.

The covariance matrix can be used to determine uncertainties (see Bevington, Numerical Recipes, ML notes.)

This function is inverse of linear relationship; it is the function you use to find X given a value of Y

ff[y_,b_,m_]:= (y-b)/m

This following calculates  the curvature matrix α; the covariance matrix is the inverse of the curvature matrix.

α = Sum[w[[i]]{{1,data[[i,1]]},{data[[i,1]],data[[i,1]]^2}},{i,1,Length[data]}];
MatrixForm[Inverse[α]]

Mathematica calculates the Covariance Matrix, but note that the CovarianceMatrix of Mathematica is the curvature matrix, α, TIMES the reduced χ2 of the original data set . (See my notes and error_test.nb)

(* Output of following two lines should match. *)
MatrixForm[(χ2min/(Length[data]-2))Inverse[α]]   
MatrixForm[lm["CovarianceMatrix"]]

β is the "vector" of first derivatives of the function ff evaluated at the point of interest. The derivatives are east to evaluate “by hand,” but just for fun I’m letting Mathematica do it.

β={D[ff[y,bb,mm],{bb,1}]/.{y->y0,bb->b,mm->m},D[ff[y,bb,mm],{mm,1}]/.{y->y0,bb->b,mm->m}}

The uncertainty (due to the CORRELATED uncertainties in m and b is given by (see Bevington, p. 59):

Sqrt[β.Inverse[α].β]  (* Uncertainty if you believe your error bars *)
Sqrt[β.σ2.β]    (* Uncertainty given by Mathematica's scaled error bars *)
Sqrt[β.Inverse[α].β]Sqrt[rχ2]

NOTE: Above gives uncertainties due to correlated fit parameters (slope & intercept).  There is also uncertainty due to ΔY; this is independent of other uncertainties, so it can simply be added in quadrature.  You can check relative contribution to uncertainties here.  You could also increase dimension of σ2 by 1 by adding single diagonal element, and adding element to β too.

Total Uncertainty:

Sqrt[β.Inverse[α].β +(ΔY D[ff[y,bb,mm],{y,1}]/.{mm->m,bb->b,y->Y})^2]

Same total uncertainty as found above via Monte Carlo methods: 1.1.

III. Naive Method of Calibration

Naive Method of Calibration

This is an example of what happens if you do the simplest WRONG analysis of the calibration data.  You eventually are going to measure a value Y, and at first you might think to fit directly to get X as a function of Y (e.g.., wavelength as a function of pixel number).  There are real problems with this, including 
   1) The uncertainty is in the y values (e.g., pixel number).  This means that you can't use the standard packages to minimize χ2.
   2) You don't really have an obvious way to find an uncertainty in your final answer.
But let's proceed anyway.  Let's do an unweigthed fit and see what happens.

yxdata = Table[{xydata[[i,2]],xydata[[i,1]]},{i,1,Length[xydata]}];

y0 = 20;

lm =LinearModelFit[yxdata,{1,x},x,Weights->w]

lm[y0]

This gives a value of 9.565 with unknown uncertainty.  The value we got above from a better analysis was X = 9.7 ± 1.1