(* Various options *)
SetDirectory[NotebookDirectory[]];
Needs["ErrorBarPlots`"];

(* Some fake data with noise *)
SeedRandom[1264645];
(* A model *)
model[x_,a_,b_]:=a x^2+b x
(* Mock data with noise *)
data=Table[{x,model[x,3,-1]+RandomVariate[NormalDistribution[0,.1]],RandomReal[{.05,.1}]},{x,0,1,0.02}];
ndata=Length[data];
(* Plot of the data and model *)
Show[ErrorListPlot[data,PlotRange->All],Plot[model[x,3,-1],{x,0,1},PlotStyle->Red],Frame->True,FrameLabel->{"x","y"},Axes->False]

(* Define the chi^2 *)
chi2[a_,b_]:=Sum[((data[[i,2]]-model[data[[i,1]],a,b])/data[[i,3]])^2,{i,1,ndata}]

(* Find the best-fit by minimizing the (χ^2) *)
fmin=FindMinimum[chi2[a,b],{a,1},{b,1}]
fmin[[1]]/(ndata-2)

(* Plot the residues *)

residues=Table[((data[[i,2]]-model[data[[i,1]],a,b])/data[[i,3]])/.fmin[[2]],{i,1,ndata}];

Histogram[residues,Automatic,"PDF",Frame->True,Epilog->Plot[(1/(Sqrt[2Pi]))Exp[-x^2/2],{x,-3,3},PlotStyle->Red][[1]]]

Histogram[RandomVariate[NormalDistribution[0,1],10000],Automatic,"PDF",Frame->True,Epilog->First@Plot[PDF[NormalDistribution[0,1],x],{x,-4,4},PlotStyle->Red]]

(* Estimation of errors *)
(* The parameters of the model *)
p={a,b};
(* The Fisher matrix // see Numerical Recipes chap. 15 *)
Fij=Table[(1/2)D[chi2[a,b],p[[i]],p[[j]]]/.fmin[[2]],{i,1,2},{j,1,2}];
Print["(Subscript[F, ij])= ",Fij//MatrixForm]
(* Covariance matrix *)
Cij=Inverse[Fij];
Print["(Subscript[C, ij])= ",Cij//MatrixForm]
(* The 1σ parameter errors *)
Print["1σ errors = ",(Sqrt[Cij])//Diagonal]

(* Contours *)

(* The values of ((Δχ)^2) // See Sigmas.nb on the site *)
Table[{"((Δχ)^2)("<>ToString[xσ]<>"σ,M="<>ToString[M]<>") =",N[2 InverseGammaRegularized[(M/2),1-Erf[((xσ)/(Sqrt[2]))]],6]},{M,1,4},{xσ,1,4}]//TableForm

(* The values of ((Δχ)^2) for 1,2,3σ *)
dx2=Table[2 InverseGammaRegularized[(M/2),1-Erf[((xσ)/(Sqrt[2]))]]/.M->2//N,{xσ,1,3}]

?ContourPlot

ContourPlot[x^2+y^2==1,{x,-1.5,1.5},{y,-1.5,1.5},AspectRatio->1,Frame->True,Axes->False]

(* Find the contours by solving (χ^2)(a,b)=((Subscript[χ, min])^2)+((Δχ)^2)*)
(* 1sigma contour *)
imp1s=ContourPlot[chi2[a,b]==fmin[[1]]+dx2[[1]],{a,2.9,3.3},{b,-1.2,-0.95},AspectRatio->1,Frame->True,Axes->False]//AbsoluteTiming;
(* 2sigma contour *)
imp2s=ContourPlot[chi2[a,b]==fmin[[1]]+dx2[[2]],{a,2.5,3.8},{b,-1.5,-0.5},AspectRatio->1,Frame->True,Axes->False]//AbsoluteTiming;
(* 3sigma contour *)
imp3s=ContourPlot[chi2[a,b]==fmin[[1]]+dx2[[3]],{a,2.0,4.3},{b,-1.9,-0.0},AspectRatio->1,Frame->True,Axes->False]//AbsoluteTiming;
Show[imp1s[[2]],imp2s[[2]],imp3s[[2]],PlotRange->All]
(* Uncomment the following if you want to save the data of the contours. Very helpful if the calculation takes for ever :) . Later, use Import[] to read it back in *)
(*Export["contour_1s_model_1.m",imp1s]*)

(* Mean errors for a and b from contours *)
bd=FindRoot[{chi2[a,b]==fmin[[1]]+2.30,(D[chi2[a,b],a])==0},{a,3.2},{b,-1.1}][[All,2]][[2]];
bu=FindRoot[{chi2[a,b]==fmin[[1]]+2.30,(D[chi2[a,b],a])==0},{a,2.2},{b,-0.6}][[All,2]][[2]];
au=FindRoot[{chi2[a,b]==fmin[[1]]+2.30,(D[chi2[a,b],b])==0},{a,3.2},{b,-1.1}][[All,2]][[1]];
ad=FindRoot[{chi2[a,b]==fmin[[1]]+2.30,(D[chi2[a,b],b])==0},{a,2.2},{b,-0.6}][[All,2]][[1]];

(* Cramer-Rao bound works!!!!!! *)
{{au-fmin[[2,1,2]],fmin[[2,1,2]]-ad}//Mean,{bu-fmin[[2,2,2]],fmin[[2,2,2]]-bd}//Mean}
(Sqrt[Cij])//Diagonal

(* Add color and other stylistic options *)
cont=Show[Graphics[Style[{Hue[0.72,.4,.9],Polygon[imp3s[[2,1,1]]]},Antialiasing->True],Frame->True],Graphics[Style[{Hue[0.52,.4,.9],Polygon[imp2s[[2,1,1]]]},Antialiasing->True],Frame->True],Graphics[Style[{Hue[0.02,.4,.9],Polygon[imp1s[[2,1,1]]]},Antialiasing->True],Frame->True],AspectRatio->1,PlotRange->{{2.6,3.6},{-1.4,-0.7}},FrameLabel->{"a","b"},Epilog->{{Black,PointSize[Medium],Point[{fmin[[2,1,2]],fmin[[2,2,2]]}]},{Green,PointSize[Medium],Point[{3,-1}]},{DotDashed,Line[{{3,-1.4},{3,-1}}],{DotDashed,Line[{{0.8,-1},{3,-1}}]}}},TextStyle->{FontFamily->"Times",FontSize->22},PlotRangeClipping->True]

(* Export in .EPS or .PDF format *)
Export["fig_contour_model_1.eps",cont,ImageSize->600]
Export["fig_contour_model_1.pdf",cont,ImageSize->600]