Question: Parsing mistake in MAPLE

In equation(19), the value for F = 14.616.. is wrong. It was produced by a bult execution of the whole program. However, executing only equation(19) in one step results in the correct answer F=24.665. The values for SSE,SSR, which are 1 x 1 matrices, are correct. Indices(1x1) provide access to the elements as constants so that arithmetic operations can be performed.

The correct result follows an execution at eqn(19) as a stand alone. An incorrect result follows execution of the entire program.

  restart:               

interface(rtablesize = 50); with(LinearAlgebra); with(plots); with(Statistics); with(ArrayTools);

Y := Matrix(12, 1, [2, 3, 2, 7, 6, 8, 10, 7, 8, 12, 11, 14]); n := ArrayNumElems(Y);

This is the data, with n rows and one column.

                                [2 ]

                                [  ]

                                [3 ]

                                [  ]

                                [2 ]

                                [  ]

                                [7 ]

                                [  ]

                                [6 ]

                                [  ]

                                [8 ]

                           Y := [  ]

                                [10]

                                [  ]

                                [7 ]

                                [  ]

                                [8 ]

                                [  ]

                                [12]

                                [  ]

                                [11]

                                [  ]

                                [14]

                            n := 12

type(Y, Matrix);

                              true

X := Matrix(12, 3, [1, 0, 2, 1, 2, 6, 1, 2, 7, 1, 2, 5, 1, 4, 9, 1, 4, 8, 1, 4, 7, 1, 6, 10, 1, 6, 11, 1, 6, 9, 1, 8, 15, 1, 8, 13]); k := NumElems(X)/n;

                             [1  0  2 ]

                             [        ]

                             [1  2  6 ]

                             [        ]

                             [1  2  7 ]

                             [        ]

                             [1  2  5 ]

                             [        ]

                             [1  4  9 ]

                             [        ]

                             [1  4  8 ]

                        X := [        ]

                             [1  4  7 ]

                             [        ]

                             [1  6  10]

                             [        ]

                             [1  6  11]

                             [        ]

                             [1  6  9 ]

                             [        ]

                             [1  8  15]

                             [        ]

                             [1  8  13]

                             k := 3

Xprimed := Transpose(X);

                                                                                                                                                                                                     Calculate Multiple Linear Regression in x1, and x2

                 [1  1  1  1  1  1  1  1   1   1  1   1 ]

                 [                                      ]

      Xprimed := [0  2  2  2  4  4  4  6   6   6  8   8 ]

                 [                                      ]

                 [2  6  7  5  9  8  7  10  11  9  15  13]

XprimedX := Xprimed . X;

                              [12   52   102 ]

                              [              ]

                  XprimedX := [52   296  536 ]

                              [              ]

                              [102  536  1004]

XprimedXinverse := MatrixInverse(XprimedX);

                      [   309          77              145 ]

                      [   ---          ---      uminus0--- ]

                      [   317          317             634 ]

                      [                                    ]

                      [    77          411             141 ]

   XprimedXinverse := [   ---         ----      uminus0----]

                      [   317         2536             1268]

                      [                                    ]

                      [       145         141       53     ]

                      [uminus0---  uminus0----      ---    ]

                      [       634         1268      634    ]

XprimedY := Xprimed . Y;

                                   [90 ]

                                   [   ]

                       XprimedY := [482]

                                   [   ]

                                   [872]

Betahat := XprimedXinverse . XprimedY;

                               [   1704   ]

                               [   ----   ]

                               [   317    ]

                               [          ]

                               [   3819   ]

                    Betahat := [   ----   ]

                               [   1268   ]

                               [          ]

                               [       815]

                               [uminus0---]

                               [       634]

Betahat := convert(Betahat, float);

                          [    5.375394322     ]

                          [                    ]

               Betahat := [    3.011829653     ]

                          [                    ]

                          [&uminus0;1.285488959]

 We choose the Null Hypothesis that Y= β0+Β1*x1+Β2*x2+ϵ . 

Performing multiple linear regression leads to the best fit linear equation:  

  Y= 5.38+3.01*x1-1.29x2

fit := plot3d(5.38+3.01*X1-1.29*X2, X1 = 0 .. 8, X2 = 2 .. 15);

data := pointplot3d({[0, 2, 2], [2, 5, 7], [2, 6, 3], [2, 7, 2], [4, 7, 10], [4, 8, 8], [4, 9, 6], [6, 9, 12], [6, 10, 7], [6, 11, 8], [8, 13, 14], [8, 15, 11]}, axes = normal, color = red, symbol = soliddiamond, symbolsize = 40, scaling = unconstrained, title = `\` Data vs Best Fit ' `);

display({data, fit});

                                                                                                               Comparison of Data and Best Fit

Yprimed := Transpose(Y); YprimedY := Yprimed . Y; BetahatPrime := Transpose(Betahat);

      Yprimed := [2  3  2  7  6  8  10  7  8  12  11  14]

                       YprimedY := [840]

BetahatPrime := [5.375394322  3.011829653  &uminus0;1.285488959]

SSE := YprimedY-BetahatPrime . XprimedY; SampleVariance := (YprimedY-BetahatPrime . XprimedY)/(n-k-1);

                   SSE := [25.4589905220000]

              SampleVariance := [3.18237381525000]

SSE (below) is a 1x1 vector which has been changed into a constant                                                  Covariance of β                                                                                                                                                     

`covβ` := Typesetting[delayDotProduct](SSE(1, 1), XprimedXinverse, true);

covβ := [

  [24.8164923384795, 6.18404501638486, &uminus0;5.82263978815458], 

  [6.18404501638486, 4.12604302229574, &uminus0;2.83100762113723], 

  [&uminus0;5.82263978815458, &uminus0;2.83100762113723, 

  2.12827523291167]]

SSR := Transpose(Betahat) . Transpose(X) . Y-n*MeanY(1, 1)^2;

                   SSR := [139.541009478000]

SST := SSR+SSE; F := SSR(1, 1)*(n-k-1)/(k*SSE);

F is found in textbook eqn (8.5)

                         SST := [165.]

                    F := [14.6160295824160]

RsquaredNonCentered := Determinant(SSR)/Determinant(SST); R := sqrt(%);

            RsquaredNonCentered := 0.845703087745454

                     R := 0.919621165342259

MeanY := Mean(Y);

                  MeanY := [7.50000000000000]

RsquaredCentered := `~`[`.`](Transpose(Betahat) . Transpose(X) . Y-n*MeanY(1, 1)^2, 1/(Yprimed . Y-n*MeanY(1, 1)^2));

            RsquaredCentered := [0.845703087745454]

                                                                                                         We note the same value for R -data not centered, and that for R-data which is centered.

Transpose(Betahat) . Transpose(X) . Y; MeanY(1, 1)^2; Yprimed . Y;

# Compare SSR, Mean of Y (data), and SST to results in the book. 

                       [814.541009478000]

                        56.2500000000000

                             [840]

                                                     Table 8.2 on Page 189 of textbook

                                    ANOVA for Overall Regression Test for Data in Table 7.1

________________________________________________________________________________________________________

 Source           df                                                   SS                                                        MS                                                 F

 _________________________________________________________________________________________________________

 Include β1    2                                        139.5410                                    69.7705                                24.665

Error           9                                          25.4590                                      2.8288

         Total                                             165.0000

_________________________________________________________________________________________

and from calculations using the defining equations in step(12) through step(13) and summarized in Step(19):

F := SSR(1, 1)*(n-k-1)/(k*SSE(1, 1)); SSE, SSR, SSR+SSE; Total := Yprimed . Y-n*MeanY(1, 1)^2;

                     F := 14.6160295824160

         [25.4589905220000], [139.541009478000], [165.]

                        Total := [165.]