R help - Nov 2010 - On efficiency, Vectorize and loops

In my last e-mails, I have asked for help regarding
1. 'defining functions inside loops'
2. 'integrating functions / vector arithmetics'
3. 'vectors out of lists?'
4. 'numerical integration'

Since some of these topics seemed to be relevant (I'm guessing by the # of
replies I got), I'm posting a modified section of my code. Any thoughts on
improvements would be highly welcome!

In the end of this code I put a little test on comparing the different
functions I must integrate. I'm really interested in efficiency since
I'll
need to do that 'inner.product' operation 95*95*5 times on each of my B
bootstrap loops, where B should be 1000 but is more likely to be 500 or 200
if I'm willing to graduate before March!

# 'code0' is used for simulation and test purposes only;
# The model is the same as in Bathia, Yao, Ziegelmann p.13
# The code used for estimation is 'code1'

rm(list=ls(all=TRUE))

n<-100 # Number of observations
d<-4   # Dimension of X
a<-0   # Lower bound
b<-1   # Upper bound
p<-5   # [See Bathia, Yao, Ziegelmann p. 5)


#######################################
# Building the error function 'epsilon'
#######################################

# Building the vector process 'Z' (which is a iiN(0,I_d) vector process)
# Time varies columnwise
Z<-matrix(rnorm(10*n),nrow=10,ncol=n)

# Defining the deterministic base functions 'zeta'
zeta<-function(u,i){
  sqrt(2)+sin(pi*i*u)
}

zeta2<-function(i){force(i);
  function(u){
    sqrt(2)+sin(pi*outer(i,u))
  }
}

zeta3<-function(u,i){
    sqrt(2)+sin(pi*outer(i,u))
}

# Building the random functional process 'epsilon'
epsilon<-lapply(1:n, function(t)local({force(t);
  function(u){
  epsilon<-0
  for (j in 1:10){
    epsilon<-epsilon+2^(1-j)*Z[j,t]*zeta(u,j)
  }
  epsilon<-epsilon
  }
}))

epsilon2<-lapply(1:n, function(t)local({force(t);
  function(u){
    crossprod(2^(1-(1:10))*Z[,t],zeta2(1:10)(u))
  }
}))

epsilon3<-lapply(1:n, function(t)local({force(t);
  function(u){
    crossprod(2^(1-(1:10))*Z[,t],zeta3(u,1:10))
  }
}))

###########################
# Building the function 'X'
###########################

# Building the vector process 'xi'. Time varies columnwise.
# Each line is a realization of an AR(1) process.
xi<-matrix(0,nrow=d,ncol=n)
for (i in 1:d){
  theta<-(-1)^i * (.9 - .5*i/d)
  xi[i,]<-arima.sim(model=list(ar=theta), n, rand.gen = rnorm)
}
rm(i, theta)

#Defining the deterministic base functions 'phi'
phi<-function(u,i){
  sqrt(2)*cos(pi*i*u)
}

phi2<-function(i){force(i);
  function(u){
  sqrt(2)*cos(pi*outer(i,u))
  }
}

phi3<-function(u,i){
  sqrt(2)*cos(pi*outer(i,u))
}

# Building the random functional process 'X'

X<-lapply(1:n, function(t)local({force(t);
  function(u){
    X<-0
    for (i in 1:d){
      X<-X+xi[i,t]*phi(u,i)
    }
  X<-X
  }
}))

X2<-lapply(1:n, function(t)local({force(t);
  function(u){
    crossprod(xi[,t],phi2(1:d)(u))
  }
}))

X3<-lapply(1:n, function(t)local({force(t);
  function(u){
    crossprod(xi[,t],phi3(u,1:d))
  }
}))

#####################################
# Defining the functional process 'Y'
#####################################

Y<-lapply(1:n, function(t)local({force(t);
  function(u){
    X[[t]](u)+epsilon[[t]](u)
  }
}))

Y2<-lapply(1:n, function(t)local({force(t);
  function(u){
    X2[[t]](u)+epsilon2[[t]](u)
  }
}))

Y3<-lapply(1:n, function(t)local({force(t);
  function(u){
    X3[[t]](u)+epsilon3[[t]](u)
  }
}))

Y4<-lapply(1:n, function(t)local({force(t);
  function(u){
  crossprod(xi[,t],sqrt(2)*cos(pi*outer(1:d,u))) +
crossprod(2^(1-(1:10))*Z[,t],sqrt(2)+sin(pi*outer(1:10,u)))
  }
}))

######################
# remove later

Ybar<-function(u){
  Ybar<-0
  for (t in 1:n){
    Ybar<-Ybar+Y[[t]](u)
  }
  Ybar<-1/n*Ybar
}

Ybar2<-function(u){
  Ybar2<-0
  for (t in 1:n){
    Ybar2<-Ybar2+Y2[[t]](u)
  }
  Ybar2<-1/n*Ybar2
}

Ybar3<-function(u){
  Ybar3<-0
  for (t in 1:n){
    Ybar3<-Ybar3+Y3[[t]](u)
  }
  Ybar3<-1/n*Ybar3
}

Ybar4<-function(u){
  Ybar4<-0
  for (t in 1:n){
    Ybar4<-Ybar4+Y4[[t]](u)
  }
  Ybar4<-1/n*Ybar4
}

# Defining the inner product on L^2([a,b])
inner.product<-function(f,g){
  core<-function(u){
    core<-f(u)*g(u)
  }
  inner.product<-integrate(core,a,b)$value
}

###########################################################
# Defining the deviation functions (Y[[t]] - Ybar) used as an input
# by the function inner.product below
Y.dev<-lapply(1:n, function(t)local({force(t);
  function(u){
    Y[[t]](u)-Ybar(u)
  }
}))

Y.dev2<-lapply(1:n, function(t)local({force(t);
  function(u){
    Y2[[t]](u)-Ybar2(u)
  }
}))

Y.dev3<-lapply(1:n, function(t)local({force(t);
  function(u){
    Y3[[t]](u)-Ybar3(u)
  }
}))

Y.dev4<-lapply(1:n, function(t)local({force(t);
  function(u){
    Y4[[t]](u)-Ybar4(u)
  }
}))

print(inner.product(Y.dev[[n]], Y.dev[[n]]))
print(inner.product(Y.dev2[[n]], Y.dev2[[n]]))
print(inner.product(Y.dev3[[n]], Y.dev3[[n]]))
print(inner.product(Y.dev4[[n]], Y.dev4[[n]]))

system.time(for(i in 1:100)inner.product(Y.dev[[n]], Y.dev[[n]]))
system.time(for(i in 1:100)inner.product(Y.dev2[[n]], Y.dev2[[n]]))
system.time(for(i in 1:100)inner.product(Y.dev3[[n]], Y.dev3[[n]]))
system.time(for(i in 1:100)inner.product(Y.dev4[[n]], Y.dev4[[n]]))


# And this is the output I get:
[1] 13.64025
[1] 13.64025
[1] 13.64025
[1] 13.64025
   user  system elapsed
  13.09    0.00   13.19
   user  system elapsed
  11.86    0.00   12.42
   user  system elapsed
  11.17    0.00   11.44
   user  system elapsed
  10.59    0.00   10.79

Thanks folks!

	[[alternative HTML version deleted]]