Python and R

Note that using solve can be numerically unstable for certain problems.

On Fri, Feb 20, 2009 at 6:50 AM, Kenn Konstabel <lebatsnok at gmail.com> wrote:

Decyphering formulas seems to be the most time consuming part of lm:

mylm1 <- function(formula, data) {
    # not perfect but works
    F <- model.frame(formula,data)
    y <- model.response(F)
    mt <- attr(F, "terms")
    x <- model.matrix(mt,F)
    coefs <- solve(crossprod(x), crossprod(x,y))
    coefs
    }

mylm2 <- function(x, y, intercept=TRUE) {
    if(!is.matrix(x)) x <- as.matrix(x)
    if(intercept) x <- cbind(1,x)
    if(!is.matrix(y)) y <- as.matrix(y)
    solve(crossprod(x), crossprod(x,y))
    }

system.time(for(i in 1:1000) mylm2(EuStockMarkets[,-1],
EuStockMarkets[,"DAX"]))

   user  system elapsed
   6.43    0.00    6.53

system.time(for(i in 1:1000) mylm1(DAX~., EuStockMarkets))

   user  system elapsed
  16.19    0.00   16.23

system.time(for(i in 1:1000) lm(DAX~., EuStockMarkets))

   user  system elapsed
  21.43    0.00   21.44

So if you need to save time, I'd suggest something close to mylm2 rather
than mylm1.

Kenn



On Thu, Feb 19, 2009 at 8:04 PM, Gabor Grothendieck
<ggrothendieck at gmail.com> wrote:

On Thu, Feb 19, 2009 at 8:30 AM, Esmail Bonakdarian <esmail.js at gmail.com>
wrote:

Thanks for the suggestions, I'll have to see if I can figure out how to
convert the relatively simple call to lm with an equation and the data
file
to the functions you mention (or if that's even feasible).

X <- model.matrix(formula, data)

will calculate the X matrix for you.

Not an expert in statistics myself, I am mostly concentrating on the
programming aspects of R. Problem is that I suspect my colleagues who
are providing some guidance with the stats end are not quite experts
themselves, and certainly new to R.

Cheers,

Esmail

Kenn Konstabel wrote:

lm does lots of computations, some of which you may never need. If
speed
really matters, you might want to compute only those things you will
really
use. If you only need coefficients, then using %*%, solve and crossprod
will
be remarkably faster than lm

# repeating someone else's example
# lm(DAX~., EuStockMarkets)

 y <- EuStockMarkets[,"DAX"]
 x <- EuStockMarkets
 x[,1]<-1
colnames(x)[1] <- "Intercept"

lm(y ~ x-1)
solve(crossprod(x), t(x))%*%y    # probably this can be done more
efficiently

# and a naive timing

 > system.time( for(i in 1:1000) lm(y ~ x-1))

  user  system elapsed
 14.64    0.33   32.69

 > system.time(for(i in 1:1000) solve(crossprod(x), crossprod(x,y)) )

  user  system elapsed
  0.36    0.00    0.36


Also lsfit() is a bit quicker than lm or lm.fit.

Regards,
Kenn