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) <[EMAIL PROTECTED]> Fax-to-email: +44 (0)870 094 0861 Date: 02-Mar-08 Time: 21:05:50 ------------------------------ XFMail ------------------------------ ______________________________________________ R-help@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.
