Fitting the inclination of a plane with the matrix formula AND Mathematica

Mathematical definition and simulation

a plane that goes through the origin is is given by z=height(x,y;a1,a2)=a1.x+a2.y

height[xy_,a_]:=a.xy

r={{1,1},{1,2},{2,1},{2,2}} (* {x,y} coordinates where the experimental values were taken *)

We are going to simulate the measurement, then we choose arbitrary true values for the parameters

a0={-1.2,1.8}

Plot3D[height[{x,y},a0],{x,0,3},{y,0,3},AxesLabel->Automatic]

We choose arbitrary standard deviations, but not all equal, so they enter effectively the calculation

σ=Range[0.1,0.25,0.05]

This is already in the format that Mathematica will work

hSim=Table[Flatten[{r[[i]],RandomVariate[ NormalDistribution[ height[r[[i]],a0],σ[[i]] ] ]}],{i,Length[r]}]

plane0=Plot3D[height[{x,y},a0],{x,0,3},{y,0,3},AxesLabel->Automatic]

planeExp=ListPointPlot3D[hSim,PlotRange->{{0,3},{0,3},{-4,4}}, AxesLabel->Automatic,PlotStyle->PointSize[Medium]]

Show[plane0,planeExp ] (* may be some points are located under the plane - turn to see *)

Fit using the Matrix formalism

Some thinking and perhaps scribbling shows that the design matrix is r.

r //MatrixForm (* just for illustration purposes *)

Now, the experimental values at these r[[i]] coordinates are

z=hSim[[All,3]]

The variance matrix is diagonal,

Notice the trick to define AND visualize the matrix - first we store the two-dimension list then we display it. This is not possible to do in a single instruction, we definitively need this compound expression.

v=DiagonalMatrix[σ^2];v//MatrixForm

Evaluate the estimated parameters using the matrix formula - this problem does not require neither care nor optimization

aHat=LinearSolve[Transpose[r].Inverse[v].r,Transpose[r].Inverse[v].z]

va=Inverse[Transpose[r].Inverse[v].r];MatrixForm[va]

hf=Table[ height[r[[i]],aHat],{i,Length[r]}]

z-hf

χ2=(z-hf).Inverse[v].(z-hf)

Round[ 100(1-CDF[ChiSquareDistribution[2],χ2]) ,0.1]

Show[Plot3D[{height[{x,y},aHat],height[{x,y},a0]},{x,0,3},{y,0,3},AxesLabel->Automatic,PlotStyle->{Red,Blue}],planeExp ]

The hat matrix and the residuals

hMatrix=r.Inverse[ Transpose[r].Inverse[v].r ] .Transpose[r].Inverse[v]

z

hMatrix.{0.5963320712532979`,2.,-1.0570294850258657`,1.059620322760166`}

hMatrix.z

Table[height[{r[[i,1]],r[[i,2]]},aHat],{i,Length[r]}]

varResidues=(IdentityMatrix[4]-hMatrix).v

Sqrt[Diagonal[varResidues]]

σ

normalizedResiduals=(z-hf)/Sqrt[Diagonal[varResidues]]

residuals = (z-hf)/σ

Fit using LinearModelFit

In LinearModelFit, it is possible to use as many independent variables as needed. Each point is a list, with values for all coordinates, listed in the 3rd parameter. The 2nd parameter gives the composition of the given coordinates that define the contribution of each parameter; they are the g-functions of the formalism we have seen last week.

TableForm[hSim,TableHeadings->{Range[1,Length[hSim] ],{"x","y", "z"}}]

fitMath=LinearModelFit[hSim,{x,y},{x,y},IncludeConstantBasis->False,Weights->1/σ^2,VarianceEstimatorFunction->(1&)]

fitMath["ParameterTable"]

aHat

Table[fitMath[r[[i,1]],r[[i,2]] ],{i,Length[r]}]

hf

fitMath["CovarianceMatrix"]//MatrixForm

va//MatrixForm

fitMath["FitResiduals"]

(z-hf)

fitMath["StandardizedResiduals"]

normalizedResiduals