From: Ajay Narottam Shah <ajayshah_at_mayin.org>

Date: Mon 06 Jun 2005 - 15:31:37 EST

Date: Mon 06 Jun 2005 - 15:31:37 EST

On my machine (ibook @ 1.2 GHz, OS X "panther", R 2.1), I get:

> measurement(ols.lf1)

[1] 7.45

> measurement(ols.lf1.inlogs)

[1] 6.9

Here is the self-contained bug-demonstration code:

ols.lf1 <- function(theta, y, X) {

beta <- theta[-1]

sigma2 <- theta[1]

if (sigma2 <= 0) return(NA)

n <- nrow(X)

e <- y - X%*%beta

logl <- ((-n/2)*log(2*pi)) - ((n/2)*log(sigma2)) - ((t(e)%*%e)/(2*sigma2))
return(-logl) # since optim() does minimisation by default.
}

# A variant: Shift to logs for sigma2 so I can do unconstrained maximisation -- ols.lf1.inlogs <- function(truetheta, y, X) { ols.lf1(c(exp(truetheta[1]), truetheta[-1]), y, X) }

# We embark on one numerical experiment

set.seed(101)

X <- cbind(1, runif(20000))

theta.true <- c(2,4,6) # error variance = 2, intercept = 4, slope = 6.
y <- X %*% theta.true[-1] + sqrt(theta.true[1]) * rnorm(20000)

measurement <- function(f) {

N <- 200

cost <- system.time(

for (i in 1:N) { j <- f(c(2,3,4), y, X) } ) return(cost[1]*1000/N) # cost per lf in milliseconds}

measurement(ols.lf1)

measurement(ols.lf1.inlogs)

-- Ajay Shah Consultant ajayshah@mayin.org Department of Economic Affairs http://www.mayin.org/ajayshah Ministry of Finance, New Delhi ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.htmlReceived on Mon Jun 06 19:09:52 2005

*
This archive was generated by hypermail 2.1.8
: Fri 03 Mar 2006 - 03:32:22 EST
*