Re: [R] Need help to locate my mistake

From: Ted Harding <Ted.Harding_at_manchester.ac.uk>
Date: Sun, 02 Mar 2008 21:05:54 +0000 (GMT)


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 ------------------------------

______________________________________________
R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun 02 Mar 2008 - 21:13:57 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Mon 03 Mar 2008 - 00:30:19 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive