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]]
______________________________________________
R-help@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.
