selection of outputs from the function

Hi,

Another option is to vectorize your function and have the theta()
function return the data frame.  On my system, this approach was
roughly 10x faster than using your original function + for loop and
gave identical results.  I also made a few other minor changes to your
code to make it easier for me to understand.  The only real changes
are:

ei <- as.vector(Y - X %*% B.cap)
my.pi <- diag(P)
data.frame(ti = ti, ti.star = ti.star, pi = my.pi,
    pi.star = pi.star, LDi = LDi, CVRi = CVRi)

## instead of

ei <- Y[i, ] - X[i, ] %*% B.cap
pi <- P[i, i]
list(ti = ti, ti.star = ti.star, pi = my.pi,
    pi.star = pi.star, LDi = LDi, CVRi = CVRi)

Hope that helps,

Josh

#### Updated theta function ####
theta2 <- function(X, Y) {
  B.cap <- solve(crossprod(X)) %*% crossprod(X, Y)
  P <- X %*% solve(crossprod(X)) %*% t(X)
  Y.cap <- P %*% Y
  e <- Y - Y.cap
  dX <- nrow(X) - ncol(X)
  var.cap <- crossprod(e) / (dX - 1)
  ei <- as.vector(Y - X %*% B.cap)
  my.pi <- diag(P)
  var.cap.i <- (((dX - 1) * var.cap)/(dX - 2)) -
    (ei^2/((dX - 2) * (1 - my.pi)))

  ti <- ei / sqrt(var.cap * (1 - my.pi))
  ti.star <- ei / sqrt(var.cap.i * (1 - my.pi))
  pi.star <- my.pi + ei^2 / crossprod(e)

  LDi <- nrow(X) *
    log((nrow(X)/(nrow(X) - 1)) * ((dX - 2)/(ti.star^2 + dX - 2))) +
    ((ti.star^2 * (nrow(X) - 1)) / ((1 - my.pi) * (dX - 2))) - 1

  CVRi <- (((dX - 1 - ti^2)/(dX - 2))^nrow(X)) / (1 - my.pi)

  output <- data.frame(ti = ti, ti.star = ti.star, pi = my.pi,
    pi.star = pi.star, LDi = LDi, CVRi = CVRi)

  return(output)
}

selection of outputs from the function

Thread (6 messages)