Fitting counting data

The model

Assume our counting data follows the function

f[x_,A_,B_]:=A Exp[-(x^2)]+B

Plot[f[x,5,1],{x,-5,5},AxesLabel->{"x","f(x)"}, AxesOrigin->{0,0}]

Our measurement outcome bins the data into bins of width 1/2. Bin centers are at half integer numbers.

The data

In a counting experiment, the counts in the bins are measured. The result is

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

xi=data[[All,1]];
ci=data[[All,2]];

Here we show the distribution:

g1=Show[Graphics[{Blue,EdgeForm[Black],Table[Rectangle[{xi[[n]]-0.25,0},{xi[[n]]+0.25,N[ci[[n]]]}],{n,1,Length[xi]}]}],Axes->False,Frame->True,AspectRatio->1, FrameTicksStyle->24,FrameLabel->{Style["(E-(Subscript[E, 0]))/Δ",24],Style["Counts",24]}]

The maximum likelihood estimate

We assume a constant prior for A and B.

The -log of the posterior distribution is then, up to a constant

L[A_,B_]:=Evaluate[-(∑_(n=1)^(Length[xi])) Log[PDF[PoissonDistribution[f[xi[[n]],A,B]/2],ci[[n]]]]];

The Contour plot reveals a minimum around A=100, B=20:

fig1=ContourPlot[L[A,B],{A,5,200},{B,1,40},Contours->20,FrameLabel->{Style["A",24],Style["B",24]},FrameTicksStyle->24]

We find the minimum numerically:

res=FindMinimum[{L[A,B],A>0,B>0},{A,100},{B,20}]

Lmin=res[[1]]

A0=A/.res[[2]]

B0=B/.res[[2]]

Here is a contour plot of the posterior probability. The contour line L(A,B)=1 gives the standard errors for A and B, if projected on the coordinate axes. We determine the errors graphically.

Show[ContourPlot[2 (L[A,B]-Lmin),{A,75,110},{B,18,22.2},Contours->{0.25,1,2,3},FrameLabel->{Style["A",24],Style["B",24]},FrameTicksStyle->24,ContourLabels->(Text[Style[#3,24],{#1,#2},Background->White]&),ContourStyle->{Black,{Red,AbsoluteThickness[2]},Black,Black}],Graphics[{Red,Disk[{A0,B0},{0.4,.045}],Dashed,Line[{{A0,18},{A0,22.2}}],Line[{{75,B0},{110,B0}}],Line[{{83.2,18},{83.2,20.4}}],Line[{{100.5,18},{100.5,19.6}}],Line[{{75,18.97},{94,18.97}}],Line[{{75,21.12},{90,21.12}}]}]
]

We determine the uncertainties from the graph:

A0-83.2
100.5-A0

B0-18.97
21.12-B0

Here is a 3D-Plot of the posterior probability:

Plot3D[Exp[-(L[A,B]-Lmin)],{A,65,120},{B,15.6,24.6},AxesLabel->{Style["A",24],Style["B",24],Style["prob(A,B)",24]},TicksStyle->24,ColorFunction->"SandyTerrain",MeshFunctions->{#3&},Mesh->{{Exp[-1/2],Exp[-1],Exp[-3/2]}},MeshStyle->{Red,AbsoluteThickness[2]},PlotPoints->50,PlotRange->All]

Now we wish to generate a plot that shows the result of the fit with the data. First a plot of the result of the fit:

g2=Plot[{f[x,A0,B0]/2,B0/2},{x,-10,10},PlotRange->All,PlotStyle->{{Red,AbsoluteThickness[2]},{Red,Dashed, AbsoluteThickness[2]}}];
Show[g1,g2,AxesLabel->{"x","counts"},AxesOrigin->{-11,0}]

Here we calculate the Hesse-Matrix:

H={{Evaluate[(∂_A) (∂_A) (L[A,B])]/.res[[2]],Evaluate[(∂_A) (∂_B) (L[A,B])]/.res[[2]]},{Evaluate[(∂_A) (∂_B) (L[A,B])]/.res[[2]],Evaluate[(∂_B) (∂_B) (L[A,B])]/.res[[2]]}}

H//MatrixForm

Its inverse is the covariance matrix:

iH=Inverse[H]

iH//MatrixForm

From the elements of the covariance matrix we find the uncertainties in the parameters and the correlation coefficient:

σA=(Sqrt[iH[[1,1]]])

σB=(Sqrt[iH[[2,2]]])

ρ = iH[[1,2]]/σA/σB