Hello Everyone!
Here's what I'm trying to do. I'm working on generating occurrences of
precipitation based upon precipitation occurrence for a station during the
previous day and two stations that have already been generated by joint
probablities and 1st order Markov chains or by the same generation process.
This has to be done for each remaining stations for each month.
> genmat # 7 stations in this example, line_before is the climatology of the
last day of the previous month. Stations 4 and 6 have been generated already
in this example
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
line_before 1 1 1 0 1 1 1
NA NA NA 1 NA 0 NA
NA NA NA 0 NA 0 NA
NA NA NA 0 NA 0 NA
NA NA NA 0 NA 0 NA
NA NA NA 0 NA 0 NA
NA NA NA 0 NA 0 NA
NA NA NA 0 NA 0 NA
NA NA NA 1 NA 0 NA
NA NA NA 1 NA 1 NA
NA NA NA 1 NA 1 NA
NA NA NA 0 NA 0 NA
NA NA NA 0 NA 0 NA
NA NA NA 1 NA 1 NA
NA NA NA 0 NA 0 NA
NA NA NA 0 NA 0 NA
NA NA NA 0 NA 0 NA
NA NA NA 0 NA 0 NA
NA NA NA 0 NA 0 NA
NA NA NA 1 NA 1 NA
NA NA NA 0 NA 0 NA
NA NA NA 1 NA 1 NA
NA NA NA 1 NA 1 NA
NA NA NA 1 NA 1 NA
NA NA NA 0 NA 0 NA
NA NA NA 0 NA 1 NA
NA NA NA 0 NA 0 NA
NA NA NA 0 NA 0 NA
NA NA NA 1 NA 1 NA
NA NA NA 1 NA 1 NA
NA NA NA 1 NA 1 NA
NA NA NA 0 NA 0 NA
> num # station to generate
[1] 2
> use1 # 1st station to use in generation
[1] 6
> use2 # 2nd station to use in generation
[1] 4
> genmat = event.gen2(genmat,use1,use2,num,ortho_obs_used) # Generation
function shown below
> genmat # genmat - after it has gone through station 2
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
line_before 1 1 1 0 1 1 1
NA 0 NA 1 NA 0 NA
NA 0 NA 0 NA 0 NA
NA 0 NA 0 NA 0 NA
NA 0 NA 0 NA 0 NA
NA 0 NA 0 NA 0 NA
NA 0 NA 0 NA 0 NA
NA 0 NA 0 NA 0 NA
NA 0 NA 1 NA 0 NA
NA 1 NA 1 NA 1 NA
NA 1 NA 1 NA 1 NA
NA 1 NA 0 NA 0 NA
NA 0 NA 0 NA 0 NA
NA 1 NA 1 NA 1 NA
NA 0 NA 0 NA 0 NA
NA 0 NA 0 NA 0 NA
NA 0 NA 0 NA 0 NA
NA 0 NA 0 NA 0 NA
NA 0 NA 0 NA 0 NA
NA 1 NA 1 NA 1 NA
NA 0 NA 0 NA 0 NA
NA 1 NA 1 NA 1 NA
NA 0 NA 1 NA 1 NA
NA 1 NA 1 NA 1 NA
NA 0 NA 0 NA 0 NA
NA 1 NA 0 NA 1 NA
NA 0 NA 0 NA 0 NA
NA 0 NA 0 NA 0 NA
NA 1 NA 1 NA 1 NA
NA 1 NA 1 NA 1 NA
NA 1 NA 1 NA 1 NA
NA 0 NA 0 NA 0 NA
Where event.gen2 is this function:
event.gen2 = function(genmat,use1,use2,num,ortho_obs_used){
for(r in 2:nrow(genmat)){
ran = runif(1,0,1)
if(genmat[r,use1]==0 & genmat[r,use2]==0){
genmat[r,num]<-ifelse(genmat[r-1,num]==0,ifelse(ran<ortho_obs_used$Pr[1],1,0),ifelse(ran<ortho_obs_used$Pr[4],1,0))
}
if(genmat[r,use1]==0 & genmat[r,use2]==1){
genmat[r,num]<-ifelse(genmat[r-1,num]==0,ifelse(ran<ortho_obs_used$Pr[2],1,0),ifelse(ran<ortho_obs_used$Pr[5],1,0))
}
if(genmat[r,use1]==1 & genmat[r,use2]==0){
genmat[r,num]<-ifelse(genmat[r-1,num]==0,ifelse(ran<ortho_obs_used$Pr[3],1,0),ifelse(ran<ortho_obs_used$Pr[7],1,0))
}
if(genmat[r,use1]==1 & genmat[r,use2]==1){
genmat[r,num]<-ifelse(genmat[r-1,num]==0,ifelse(ran<ortho_obs_used$Pr[6],1,0),ifelse(ran<ortho_obs_used$Pr[8],1,0))
}
gc()
}
genmat
}
####
ortho_obs_used is a data frame that contains the probablity of precipitation
occurring on a given day for a specific set of condtions.
For instance ortho_obs_used$Pr[1] is the probablity of rain at station s for
day d, given that there was no rain at station s for day d-1 and no rain at
either of the other two stations for day d.
The event.gen2 function handles the generation, and it runs quickly for the
5 remaining stations and one month, but I have to run this for 317 stations
over 48 months or more, and it becomes a really bad bottleneck. So what I'd
like to know is if there is anyway that I can re-write this function to work
without a loop. I couldn't find anything from previous posts about getting
out of loops where the previous iteration is required to determine the next
calculation.
Sorry for the length of the post, but I thought it best to try to explain
what I was doing first, before diving into my question
Thanks in advance!
Adrienne Wootten
Graduate Research Assistant/Environmental Meteorologist
M.S. Atmospheric Science
NC State University
State Climate Office of North Carolina
Raleigh, NC 27695
[[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.
