MATH 376 -- Probability and Statistics 2

April 5, 2006

Example 1.

Normal equations for fitting the linear model

 Z = beta[0]+beta[1]*x+beta[2]*y+epsilon

to the data points [0,0,2], [1,0,1], [0,1,4], [1,1,5]

>    X:=matrix([[1.0,0,0],[1,1,0],[1,0,1],[1,1,1]]);

X := matrix([[1.0, 0, 0], [1, 1, 0], [1, 0, 1], [1, 1, 1]])

>    Z:=matrix([[2],[1],[4],[5]]);

Z := matrix([[2], [1], [4], [5]])

>    XtX:=multiply(transpose(X),X);

XtX := matrix([[4.00, 2., 2.], [2., 2, 1], [2., 1, 2]])

>    XtZ:=multiply(transpose(X),Z);

XtZ := matrix([[12.0], [6], [9]])

>    beta:=linsolve(XtX,XtZ);

beta := matrix([[1.500000000], [0.], [3.]])

>    PlP:=plot3d(beta[1,1]+beta[2,1]*x+beta[3,1]*y,x=-0.5..1.5,y=-0.5..1.5,axes=boxed):

>    SP1:=plot3d([0.05*cos(theta)*sin(phi),0.05*sin(theta)*sin(phi),2+0.05*cos(phi)],theta=0..2*Pi,phi=0..Pi,color=blue):

>    SP2:=plot3d([1+0.05*cos(theta)*sin(phi),0.05*sin(theta)*sin(phi),1+0.05*cos(phi)],theta=0..2*Pi,phi=0..Pi,color=blue):

>    SP3:=plot3d([0.05*cos(theta)*sin(phi),1+0.05*sin(theta)*sin(phi),4+0.05*cos(phi)],theta=0..2*Pi,phi=0..Pi,color=blue):

>    SP4:=plot3d([1+0.05*cos(theta)*sin(phi),1+0.05*sin(theta)*sin(phi),5+0.05*cos(phi)],theta=0..2*Pi,phi=0..Pi,color=blue):

>    display3d(SP1,SP2,SP3,SP4,PlP,scaling=constrained);

[Maple Plot]

We see the plane in R^3   that comes closest to containing the 4 data points.

Example 2.

Normal equations for fitting the linear model

  Y = beta[0]+beta[1]*x+beta[2]*x^2+epsilon

to the data points  [0,1], [1,4], [2,7], [3,8].

>    with(linalg):

We let X = matrix([[1, 0, 0], [1, 1, 1], [1, 2, 4], [1, 3, 9]])    and   Y = matrix([[1], [4], [7], [8]])

The normal equations are

X^t*X*beta = X^t*Y

>    X:=matrix([[1.0,0,0],[1,1,1],[1,2,4],[1,3,9]]);

X := matrix([[1.0, 0, 0], [1, 1, 1], [1, 2, 4], [1, 3, 9]])

>    Y:=matrix([[1],[4],[7],[8]]);

Y := matrix([[1], [4], [7], [8]])

>    XtX:=multiply(transpose(X),X);

XtX := matrix([[4.00, 6., 14.], [6., 14, 36], [14., 36, 98]])

>    XtY:=multiply(transpose(X),Y);

XtY := matrix([[20.0], [42], [104]])

>    beta:=linsolve(XtX,XtY);

beta := matrix([[.9000000000], [3.900000000], [-.5000000000]])

>    with(plots):

>    Pts:=[[0,1],[1,4],[2,7],[3,8]];

Pts := [[0, 1], [1, 4], [2, 7], [3, 8]]

>    PP:=plot(Pts,style=point,symbol=circle,color=blue):

>    RP:=plot(beta[1,1] + beta[2,1]*x + beta[3,1]*x^2,x=-1..4):

>    display(PP,RP);

[Maple Plot]

>