Hello
I'm a population ecologist and use R for all my stats and modelling.
Recently I have been using R to numerically solve integral projection
models. This involves constructing several large matrices. The current code
by Easterling (Size-specific sensitivity: Applying a new structured
population model. Ecology, 2000, 81, 694-708) uses nested loops to construct
the matrices. To speed up the code I decided to replace these with outer
functions. One of these works just as expected and the other almost works!
Having spent 2 days trying to work out what's going wrong I have decided to
ask for help. The R-code is given below, the first part sets various
parameters, then there are several functions. The main functions are called
bigmatrix and new.bigmatrix which construct the matrices, these should give
the same answers, as all I've done is replace the nested loops with outer
functions. The Bmatrix calculation is correct but the Pmatrix calculation
gets some of the answers wrong. The final lines of code test the functions.
Any advice on what's going wrong would be greatly appreciated, also any tips
for speeding up the code would fantastic.
Many thanks
Mark
Dr Mark Rees
Imperial College, Silwood Park,
Ascot, Berks, SL5 7PY
Direct line 0207 5942488
Fax 0207 594 2339
#Parameters
minsize<-log(0.1); maxsize<-log(32);
p.vec.names<-rep(NA,11)
p.vec<-rep(0,11);
p.vec[1]<- 0.36 ; p.vec.names[1]<-"1st survival param";
p.vec[2]<- 0.17 ; p.vec.names[2]<-"2nd survival param";
p.vec[3]<- -18.27; p.vec.names[3]<-"1st flow param ";
p.vec[4]<- 6.91 ; p.vec.names[4]<-"2nd flow param ";
p.vec[5]<- 0.96 ; p.vec.names[5]<-"ag ";
p.vec[6]<- 0.59 ; p.vec.names[6]<-"bg ";
p.vec[7]<- 0.67^2; p.vec.names[7]<-"sigma2 growth ";
p.vec[8]<- 1.036; p.vec.names[8]<-"intercept seeds ";
p.vec[9]<- 2.22; p.vec.names[9]<-"slope seeds ";
p.vec[10]<- log(10)*(-0.304+0.27);
p.vec.names[10]<-"mean kids size ";
p.vec[11]<- (log(10)^2)*(0.172^2+0.28^2);
p.vec.names[11]<-"sigma2 kids size ";
p.est<-0.48*0.0205 #actual
# Part (II) ##############################################################
# Compute the kernel component functions from the fitted models
#survival function
sx<-function(x) {
u<-p.vec[2]*x+p.vec[1];
if(u>=0.7) u<-0.7;
if(u<0) u<-0;
return(u);
}
#flowering functiom
fx<-function(x) {
u<-exp(p.vec[3]+p.vec[4]*x);
return(u/(1+u));
}
#growth function
gxy<-function(x,y) {
sigmax2<-p.vec[7];
sigmax<-sqrt(sigmax2);
mux<-p.vec[5]+p.vec[6]*x;
fac1<-sqrt(2*pi)*sigmax;
fac2<-((y-mux)^2)/(2*sigmax2);
return(exp(-fac2)/fac1);
}
#survival-growth function
pxy<-function(x,y) { return(sx(x)*(1-fx(x))*gxy(x,y)) }
#fecundity function
fxy<-function(x,y) {
nkids<-p.est*exp(p.vec[8]+p.vec[9]*x);
kidsize.mean<- p.vec[10];
kidsize.var<- p.vec[11];
fac1<-sqrt(2*pi)*sqrt(kidsize.var);
fac2<-((y-kidsize.mean)^2)/(2*kidsize.var);
f<-fx(x)*nkids*exp(-fac2)/fac1;
return(f);
}
# Part (III)
################################################################
# THE KERNEL K(y,x)
# note the chopoff
Kyx<-function(y,x) {
xeval<-max(x,minsize); xeval<-min(xeval,maxsize);
yeval<-max(y,minsize); yeval<-min(yeval,maxsize);
#return(pxy(xeval,yeval)+fxy(xeval,yeval))
return(pxy(x,y)+fxy(x,y))
};
############## The 'big matrix' M of size n x n this uses outer and
doesn't
work, don't know why??!!
new.bigmatrix<-function(n) {
# upper and lower integration limits
L<-minsize; U<-maxsize;
# boundary points b and mesh points y
b<-L+c(0:n)*(U-L)/n;
y<-0.5*(b[1:n]+b[2:(n+1)]);
# construct the matrix
I<-diag(n);
M<-matrix(0,n,n);
P<-matrix(0,n,n);
B<-matrix(0,n,n);
P<-t(outer(y,y,pxy))
B<-t(outer(y,y,fxy))
M<-P+B
K<-M;
M<-(U-L)*M/n;
P<-(U-L)*P/n;
B<-(U-L)*B/n;
return(list(matrix=M,kernel=K,meshpts=y,Pmatrix=P,Bmatrix=B,Imatrix=I));
}
bigmatrix<-function(n) {
# upper and lower integration limits
L<-minsize; U<-maxsize;
# boundary points b and mesh points y
b<-L+c(0:n)*(U-L)/n;
y<-0.5*(b[1:n]+b[2:(n+1)]);
# loop to construct the matrix
I<-diag(n);
M<-matrix(0,n,n);
P<-matrix(0,n,n);
B<-matrix(0,n,n);
for (i in 1:n){
for(j in 1:n){
M[i,j]<-Kyx(y[i],y[j])
P[i,j]<-pxy(y[j],y[i])
B[i,j]<-fxy(y[j],y[i])
}
}
K<-M;
M<-(U-L)*M/n;
P<-(U-L)*P/n;
B<-(U-L)*B/n;
return(list(matrix=M,kernel=K,meshpts=y,Pmatrix=P,Bmatrix=B,Imatrix=I));
}
#Test code
M<-bigmatrix(5)
N<-new.bigmatrix(5)
M$Bmatrix==N$Bmatrix
M$Pmatrix==N$Pmatrix
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._