Linear Model Fit

Simple simulation to check properties of the linear model

Preparations for the matrix method

A=Table[{1,0.025+i*0.05},{i,0,19}]

xvec=#2&@@@A

σ=0.2*(1.0-2.0 #+2.0*((#)^2))&/@xvec

Create covariance matrix (the data) and weight matrix

V=DiagonalMatrix[(σ^2)];
W=Inverse[V];

Calculate the covariance matrix of the parameters and their uncertainties

Note: The calculation only needs the errors in the data and the x values, not the actual data

Va=Inverse[Transpose[A].W.A]
(Sqrt[Diagonal[Va]])
TableForm[{%},TableHeadings->{None,{"(Subscript[σ, 0])","(Subscript[σ, 1])"}}]

Simulation of the data and determination of the fit parameters and the (χ^2)

(Subscript[y, i])=(Subscript[a, 0])+(Subscript[a, 1])(Subscript[x, i])+(Subscript[ϵ, i])
(Subscript[a, 0])=2
(Subscript[a, 1])=-3
(Subscript[ϵ, i]) will be taken from a normal distribution  𝒩(0,(Subscript[σ, i]))

yvec=Transpose[{RandomReal[NormalDistribution[2-3 #1,#2]]&@@@Transpose[{xvec,σ}]}]
Flatten[{a=Va.Transpose[A].W.yvec,
Transpose[yvec].W.yvec-Transpose[a].Transpose[A].W.yvec}]

The whole thing is done 100,000 times

n=100000;data=Table[yvec=Transpose[{RandomReal[NormalDistribution[2-3 #1,#2]]&@@@Transpose[{xvec,σ}]}];Flatten[{a=Va.Transpose[A].W.yvec,
Transpose[yvec].W.yvec-Transpose[a].Transpose[A].W.yvec}],{n}];

Split the fiter results into separate tables

data0=#1&@@@data;
data1=#2&@@@data;
data2=#3&@@@data;

Parameter (Subscript[a, 0]) (y-intercept)

Histogram[data0]
Mean[data0]
StandardDeviation[data0]

Good agreement with the value determined above: (Subscript[σ, 0])=0.0654339

Parameter (Subscript[a, 1]) (slope)

Histogram[data1]
Mean[data1]
StandardDeviation[data1]

Good agreement with the value determined above: (Subscript[σ, 1])=0.118387

SumOfSquares

Histogram[data2,{0.5},Epilog->First@Plot[0.5n PDF[ChiSquareDistribution[18],x],{x,0,60},PlotStyle->Red]]
Mean[data2]
Variance[data2]

Good agreement with a  (χ^2) distribution with 18 d.o.f

Median

median=18(((1-(2/(9 18))))^3)//N

Split into 2 lists:small and large (χ^2)

dataL=Select[data,#[[3]]<median&];
dataH=Select[data,#[[3]]>median&];

Histogram[{dataL[[All,1]],dataH[[All,1]]}]

Note: The (χ^2) does not say anything about how close the parameters are to the “true” value.

Covariance ellipse

{u,s,v}=SingularValueDecomposition[Va]
ellipse=u.(Sqrt[Diagonal[s]] {Cos[t],Sin[t]})

Show[ListPlot[data[[1;;10000,1;;2]],PlotStyle->{Blue,PointSize[0.002]}],
ParametricPlot[ellipse+{2,-3},{t,0,2π},PlotStyle->Red]]

Data points within the ellipse

dataE=Select[data,Total[((1/Sqrt[Diagonal[s]]) (Transpose[u].(#[[1;;2]]-{2,-3})))^2]<1&];
Dimensions[dataE]
CDF[ChiSquareDistribution[2],1]//N

It is expected that 39.3% of the points are within the ellipse

Show[ListPlot[data[[1;;10000,1;;2]],PlotStyle->{Blue,PointSize[0.002]}],ListPlot[dataE[[1;;3935,1;;2]],PlotStyle->{Yellow,PointSize[0.002]}],
ParametricPlot[ellipse+{2,-3},{t,0,2π},PlotStyle->Red]]