Example with data

Needs["ErrorBarPlots`"]

Reading the data

These are the known errors of the data:

σx=0.5;σy=0.8;

data=Import["/Users/ihn/Documents/Teaching/VP-Leitung/DataAnalysis/2014/9. Vorlesung/errorxydata.csv","CSV"];

Grid[Prepend[data,{Style["x",Bold],Style["y",Bold]}],Background->{None,{Gray,{LightGray,White}}},Dividers->{Black,{Black,Black}},Frame->True]

nn=Length[data];

Plotting the data

datae=Table[{data[[n]],ErrorBar[σx,σy]},{n,1,nn}];
g1=ErrorListPlot[datae,PlotRange->{{-5,2*nn+7-10},{5,2*nn+7}},AspectRatio->1,AxesLabel->{Style["x",24],Style["y",24]}]

Statistics of the data

xav=Mean[data[[All,1]]]
yav=Mean[data[[All,2]]]
Varx=CentralMoment[data[[All,1]],2]
Vary=CentralMoment[data[[All,2]],2]
rho =CentralMoment[data,{1,1}]/(Sqrt[Varx Vary])

-log-posterior for a and R

LL[a_,R_,σx_,σy_,ρ_,n_]:=Log[R a]+(n/2)Log[((σx)^2)(R^2) a+((σy)^2)(R^2)/a+((σx)^2)((σy)^2)]+(n/2)(((R^2)(a^2)-2 (R^2)ρ a+(((σx)^2)+((σy)^2))a+(R^2))/(((σx)^2)(R^2)(a^2)+((σx)^2)((σy)^2)a+(R^2)((σy)^2)))

Minimizing the posterior for a and R

res1=FindMinimum[LL[aa,R,σx/(Sqrt[Varx]),σy/(Sqrt[Vary]),rho,nn],{{aa,1},{R,1}}]

LLmin=res1[[1]];
aest=aa/.res1[[2]]
Resti=R/.res1[[2]]

Plotting the -log-posterior for a and R

We determine the uncertainties of a graphically from this plot:

Show[ContourPlot[LL[aa,R,σx/(Sqrt[Varx]),σy/(Sqrt[Vary]),rho,nn]-LLmin,{aa,0.935,1.07},{R,0.68,1.52},Contours->{0.25,1,2,3},FrameLabel->{Style["α'",24],Style["R'",24]},FrameTicksStyle->24,ContourLabels->(Text[Style[#3,18],{#1,#2},Background->White]&),ContourStyle->{Black,{Red,AbsoluteThickness[2]},Black,Black}],Graphics[{{Red,Disk[{aest,Resti},{0.01(1.07-0.935),0.01(1.52-0.68)}]},{Red,Dashed,Line[{{aest+0.0375,0.68},{aest+0.0375,1.52}}],Line[{{aest-0.0355,0.68},{aest-0.0355,1.52}}]}}]]

Δap = 0.0375;
Δam = 0.0355;

Estimate of the slope and the intercept parameters α and β

Slope parameter:

α=aest (Sqrt[Vary])/(Sqrt[Varx])

Δαp=Δap (Sqrt[Vary])/(Sqrt[Varx])
Δαm=Δam (Sqrt[Vary])/(Sqrt[Varx])

Intercept parameter:

β=yav-xav α

Δβ=(Sqrt[((((σy)^2)+((σx)^2)α)/nn)+(xav^2)Δαp Δαm])

Plotting the result of the fit

Show[g1,Plot[{α (x-xav) + yav,(α+Δαp) (x-xav) + yav,(α-Δαm) (x-xav) + yav},{x,-5,2*nn+7-10},PlotStyle->{{Red,AbsoluteThickness[2]},{Red,Dashed},{Red,Dashed}}]]