[R] Need help to locate my mistake

[(Ted Harding)]

On 02-Mar-08 20:18:29, Louise Hoffman wrote:
> Dear readers
> 
> I would like to make General Linear Model (GLM) for the
> following data set http://louise.hoffman.googlepages.com/fuel.csv
> 
> The code I have written is
> 
> fuelData<-read.table('fuel.csv',header=TRUE, sep=',')
> n<-dim(fuelData)[1]
> xOnes<- matrix(1,nrow=n,ncol=1)
> x<-cbind(xOnes,fuelData[,3])
> y<-fuelData[,4]
> theta<-((t(x)%*%x)^(-1))%*%t(x)%*%y
> 
> which gives
> 
>> theta
>             [,1]
> [1,] 215.8374077
> [2,]   0.1083491
> 
> When I do it in Matlab I get theta to be
> 79.69
> 0.18
> 
> which is correct. ~79 is the crossing of the y-axis.
> 
> Have I perhaps written theta wrong? The formula for theta is
> (alpha,beta)^T = (x^T * x)^(-1) * x^T * Y
> 
> where ^T means transposed.

Unfortunately, x^(-1) is not the inverse of x:

x<-matrix(c(2,4,4,5),nrow=2)
x
#      [,1] [,2]
# [1,]    2    4
# [2,]    4    5

x^(-1)
#      [,1] [,2]
# [1,] 0.50 0.25
# [2,] 0.25 0.20

i.e. it is the matrix which you get by applying the operation
(...)^(-1) to each element of x.

In R, the inverse of a non-singular matrix is (somewhat obscurely)
denoted by solve(x):

solve(x)
#            [,1]       [,2]
# [1,] -0.8333333  0.6666667
# [2,]  0.6666667 -0.3333333

solve(x)%*%x
#      [,1]         [,2]
# [1,]    1 1.110223e-16
# [2,]    0 1.000000e+00
(Note the slight rounding error); whereas

(x^(-1))%*%x
#      [,1] [,2]
# [1,]  2.0 3.25
# [2,]  1.3 2.00

Best wishes,
Ted.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at manchester.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 02-Mar-08                                       Time: 21:05:50
------------------------------ XFMail ------------------------------