#I am new with R programmin, I am trying to apply a monte carlo method by sampling every #time from a uniform distribution to get the final biomass that fall within the range that I put in #the condition #This takes more then 4 hours running , while when I use ready package for the function it only #take some minutes
## Prior r range
n=15000
rmin=0.1
rmax=0.9
r <- runif(n, rmin, rmax)
#Prior range for K
kmin=max(Catch)
kmax=2000000
k <- runif(n, kmin, kmax)
rk<-cbind(r,k)
#prior range for initial depletion and final depletion
dstart=c(0.8,1)
dinit<-seq(dstart[1],dstart[2],0.1)
dfin=c(0.2,0.4)
#process error
sigma<-0.1
vt<-rnorm(1,0,sigma)
2.2. Monte Carlo filtering with schaefer function
B<-matrix(NA,nrow=length(k),ncol=nyr, byrow=TRUE)
for (j in 1:length(k)){
for (i in dinit){
B[j,1]<-k[j]*i*exp(vt)
for (t in 1:(nyr-1)){
B[j,t+1] = B[j,t]+r[j]*B[j,t]*(1-(B[j,t]/k[j]))*exp(vt)-Catch[t]
matall<-cbind( rk,B)
lastyear<-ifelse((matall[,ncol(matall)]/k)>dfin[1] & (matall[,ncol(matall)]/k)<dfin[2],matall[,ncol(matall)],NA)
}
}
}
matall2<-cbind(matall[,1:ncol(matall)-1],lastyear)
Res<-na.omit(matall2)
Hi, loops (specially when you have three nested loops) can be very slow (depending on the speed of the processor, number of iterations, operations inside them). Your code runs roughly, as there is not enough data on your post, please have a look at
With 5 years of dummy information it runs in 2 minutes in my desktop, so it will also depends on how many years you have of Catches.
But anyway, your MC loop is wrong. You have the following line inside the third loop:
This line output is not associated to anything (nor j, t or i), neither catch by 'cibinding' to a vector/matrix so any time that you compute it, you lost it. The value that then you are 'cibinding' to matall2 is a single value (of the last iteration of the loop) and it would be the same for all rows. If this is not wrong, but a desired property of your MC, then you can speed up by calculating it only on the last iteration. In my home example, from two minutes if goes down to 30 secs. I just added a row with some conditions to execute lasteyar:
B<-matrix(NA,nrow=length(k),ncol=nyr, byrow=TRUE)
for (j in 1:length(k)){
for (i in dinit){
B[j,1]<-k[j]*i*exp(vt)
for (t in 1:(nyr-1)){
B[j,t+1] = B[j,t]+r[j]*B[j,t]*(1-(B[j,t]/k[j]))*vte-Catch[t]
matall<-cbind( rk,B)
if (j == length(k) & i == max(dinit) & t == (nyr-1)) # that's the line to be added
lastyear<-ifelse((matall[,ncol(matall)]/k)>dfin[1] & (matall[,ncol(matall)]/k)<dfin[2],matall[,ncol(matall)],NA)
}
}
}
The results of your loop and mine are identical (matall2 with your loop, matall3 with my modification)
identical(matall3, matall2)
[1] TRUE
cheers
Last minute edition, I forgot to mention that you include in your loop exp(vt) B[j,1]<-k[j]*i*exp(vt). As vt is a constant, I put it outside vte <- exp(vt) and in the loop (only in one place where it is computed). The correct way should be:
B<-matrix(NA,nrow=length(k),ncol=nyr, byrow=TRUE)
for (j in 1:length(k)){
for (i in dinit){
B[j,1]<-k[j]*i*vte
for (t in 1:(nyr-1)){
B[j,t+1] = B[j,t]+r[j]*B[j,t]*(1-(B[j,t]/k[j]))*vte-Catch[t]
matall<-cbind( rk,B)
if (j == length(k) & i == max(dinit) & t == (nyr-1)) # that's the line to be added
lastyear<-ifelse((matall[,ncol(matall)]/k)>dfin[1] &
(matall[,ncol(matall)]/k)<dfin[2],
matall[,ncol(matall)],NA)
}
}
}
And I got a similar increase of performance (now in my laptop, from 7 and a half minutes to 2 and a half minutes) with a toy example. So I can see that you can easily go for hours with real datasets. Performance on loops depends on the processor speed (and laptops are normally not very good at it, even if the processors are reported with high CPU clock speed).
If by using a package makes it much faster it is most likely because the package has compiled code (C, or c++ via Rcpp, or, as it seems you are working with catch data, maybe TMB, that is quite used in Fisheries/stock analysis, and really fast)
You should check carefully the loop. I have just another view and found another bottleneck, The following lines:
for (i in dinit){
B[j,1]<-k[j]*i*vte
Does nothing but calculate B[j,1] several times (for each value of dinit), overwritting it, so you have a bug there. As before, if it is not a bug, you can eliminate the loop as I have done below (and now it takes 45 seconds, vs the original 7.5 minutes, for achieving the same results)
s <- Sys.time()
B<-matrix(NA,nrow=length(k),ncol=nyr, byrow=TRUE)
for (j in 1:length(k)){
B[j,1]<-k[j]*dinit[length(dinit)]*vte
for (t in 1:(nyr-1)){
B[j,t+1] = B[j,t]+r[j]*B[j,t]*(1-(B[j,t]/k[j]))*vte-Catch[t]
matall<-cbind( rk,B)
if (j == length(k) & i == max(dinit) & t == (nyr-1)) # that's the line to be added
lastyear<-ifelse((matall[,ncol(matall)]/k)>dfin[1] &
(matall[,ncol(matall)]/k)<dfin[2],
matall[,ncol(matall)],NA)
}
}
And finally, you can go a bit faster taking the first operation out of the loop (37 seconds now):
B<-matrix(NA,nrow=length(k),ncol=nyr, byrow=TRUE)
B[,1]<-k*dinit[length(dinit)]*vte
for (j in 1:length(k)){
for (t in 1:(nyr-1)){
B[j,t+1] = B[j,t]+r[j]*B[j,t]*(1-(B[j,t]/k[j]))*vte-Catch[t]
matall<-cbind( rk,B)
if (j == length(k) & i == max(dinit) & t == (nyr-1)) # that's the line to be added
lastyear<-ifelse((matall[,ncol(matall)]/k)>dfin[1] &
(matall[,ncol(matall)]/k)<dfin[2],
matall[,ncol(matall)],NA)
}
}
But alas, you didn't provide some reproducible example to have a more real look, I may be wrong, but I think the problem is in your code.
Many people discourage other for using loops (I like them, and If you do bayesian approaches you won't survive without them). They are beautiful, but tricky...