[R] selection of outputs from the function

[Joshua Wiley]

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

On Fri, Dec 24, 2010 at 8:56 AM, ufuk beyaztas <ufukbeyaztas at gmail.com> wrote:
>
> Hi Dear All,
> This is a function which contains Covariance Ratio and Likelihood Distance
> values (CVRi, LDi). i want to compute the all row's values, that is run this
> function for nrow(X) times. The X and Y matrices are;
>
> X<-matrix(c(1125,920,835,1000,1150,990,840,650,640,583,570,570,510,555,460,275,510,165,244,79,232,268,271,237,192,202,184,200,180,165,151,171,243,147,286,198,196,210,327,334,7160,8804,8108,6370,6441,5154,5896,5336,5041,5012,4825,4391,4320,3709,3969,3558,4361,3301,2964,2777,85.9,86.5,85.2,83.8,82.1,79.2,81.2,80.6,78.4,79.3,78.7,78.0,72.3,74.9,74.4,72.5,57.7,71.8,72.5,71.9,8905,7388,5348,8056,6960,5690,6932,5400,3177,4461,3901,5002,4665,4642,4840,4479,4200,3410,3360,2599),nrow=20)
> Y<-matrix(c(1.5563,0.8976,0.7482,0.7160,0.3130,0.3617,0.1139,0.1139,-0.2218,-0.1549,0.0000,0.0000,-0.0969,-0.2218,-0.3979,-0.1549,-0.2218,-0.3979,-0.5229,-0.0458),nrow=20)
>
> theta <- function(X,Y) {
> B.cap <- solve(t(X)%*%X)%*%t(X)%*%Y
> P <- X%*%solve(t(X)%*%X)%*%t(X)
> Y.cap <- P%*%Y
> e <- Y-Y.cap
> var.cap<-(t(e)%*%e)/(nrow(X)-ncol(X)-1)
> ei <- Y[i,]-X[i,]%*%B.cap
> pi <- P[i,i]
> var.cap.i <-
> (((nrow(X)-ncol(X)-1)*var.cap)/(nrow(X)-ncol(X)-2))-(ei^2/((nrow(X)-ncol(X)-2)*(1-pi)))
> ti <- ei/(sqrt(var.cap)*sqrt(1-pi))
> ti.star <- ei/(sqrt(var.cap.i)*sqrt(1-pi))
>
> X.star <- cbind(X,Y)
> pi.star <- pi+(ei^2/(t(e)%*%e))
>
>
> LDi <-  nrow(X)*log(((nrow(X))/(nrow(X)-1))*(((nrow(X)-ncol(X)-2))/
> (ti.star^2+nrow(X)-ncol(X)-2)))+((ti.star^2*(nrow(X)-1))/((1-pi)*(nrow(X)-ncol(X)-2)))-1
> CVRi <- (((nrow(X)-ncol(X)-1-ti^2)/(nrow(X)-ncol(X)-2))^(nrow(X)))/(1-pi)
>
> list(ti=ti,ti.star=ti.star,pi=pi,pi.star=pi.star,LDi=LDi,CVRi=CVRi) }
>
> obj<-list()
> for(i in 1:nrow(X)){
> X<-X
> Y<-Y
> out<-theta(X,Y)
> obj<-c(obj,list(out))}
> obj
>
> Finally i get values...
> Is there any way to get the outputs as a list or data.frame like
>    pi CVRi
> 1   1    1
> 2   2    2
> 3   3    3
> 4   4    4
> 5   5    5
> 6   6    6
> 7   7    7
> 8   8    8
> 9   9    9
> 10 10   10
> 11 11   11
> 12 12   12
> 13 13   13
> 14 14   14
> 15 15   15
> 16 16   16
> 17 17   17
> 18 18   18
> 19 19   19
> 20 20   20
> for all values (pi,pi.star,ti,ti.star,CVRi,LDi)...
> Thanks so much for any idea !
> --
> View this message in context: http://r.789695.n4.nabble.com/selection-of-outputs-from-the-function-tp3163361p3163361.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> 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.
>



-- 
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.com/