fit_beamabs.nb

Least squares fit to determine the absorption coefficient

Version of 12 Apr 2010 (Sasa)

1. Data

pjvoda={11.1,10.1,9.2,8.3,7.5,6.8,5.9,5.1}

pjnig={9.2,6.8,5.1,3.7,2.7,2.1,1.5,1.1}

yj=pjnig/pjvoda

xj=0.01Range@Length@yj

data=Transpose@{xj,yj}

n=Length[data]

2. Model

params={α}

modelfn[x_,{α_}]:=Exp[-α x]

s2fn[data_,params_]:=Module[{m,d,x=First/@data,y=Last/@data},
m=modelfn[#,params]&/@x;
d=y-m;
d.d
]

m=Length[params]

3. Fit

pfit=FindFit[data,modelfn[x,params],params,x]

D[s2fn[data,params],#]&/@params/.pfit//N

Plot[modelfn[x,params/.pfit],{x,Min[First/@data],Max[First/@data]},PlotStyle->Thick,Frame->True,Epilog->{Red,AbsolutePointSize[5],Point/@Transpose@{First/@data,Last/@data}}]

4. Parameter errors

s2min=s2fn[data,params/.pfit]

rms=(Sqrt[(s2min/(n-m))])

atamat=(1/2)Outer[D[((n-m)/s2min) s2fn[data,params],#1,#2]&,params,params]/.pfit//N

covmat=Inverse[atamat]

pstderr=Table[covmat[[i,i]],{i,m}]//Sqrt

5. Summary

MapThread[Print[StringJoin[#1," = ",#2," +/- ",#3]]&,Map[ToString,{params,params/.pfit,pstderr},{2}]];