[R] Alternating between "for loops"

Rui Barradas ruipbarradas at sapo.pt
Wed Aug 1 21:52:11 CEST 2012


Hello,

Em 01-08-2012 20:02, Mercier Eloi escreveu:
> Your code almost gave me a headache. :-/

Agree. There's a good way of avoiding headaches, package formatR, 
function tidy.source.

My simplification is different in many places.
A common one is to treat 'observable' as a logical variable, not as a 
numeric one.
I've simplified some tests, the syntax in some places (namely, the 
condition 'if' is a function).
And eliminated some calls to runif(), whenever they could be done only once.

I think my code is the equivalent of Claudia's and I've worked it all. 
here it goes but without comments.



# # # Robust design
# # # with Markovian
# # # emigration and
# # # dummy time
# # # periods
J <- 24
N <- 10
S <- 0.9
PsiAD <- 0.06
PsiAd <- 0.04
PsiAA <- 0.4
PsiaA <- 0.3
p <- 0.5
c <- p
y <- matrix(0, N, J)
y[, 1] <- "A"
dtp <- rep(c(TRUE, TRUE, FALSE, FALSE), 6)[seq_len(J)]
for (j in which(dtp)) {
     for (q in 1:N) {
         (observable <- TRUE)
         if (j %% 2) {
             survive <- runif(1, 0, 1)
             decide <- runif(1, 0, 1)
             if (survive <= S) {
                 if (observable) {
                   observable <- (decide <= PsiAA)
                 } else {
                   observable <- (decide <= PsiaA)
                 }
                 if (observable) {
                   if (runif(1, 0, 1) <= p) y[q, j] <- "A"
                 }
             } else {
                 y[q, j] <- if (decide <= PsiAd) "d" else "D"
                 break
             }
         } else {
             if (observable) {
                 if (runif(1, 0, 1) <= c) y[q, j] <- "A"
             }
         }
     }
}

for (j in which(!dtp)) {
     for (q in 1:N) {
         if (j %% 2) {
             decide <- runif(1, 0, 1)
             if (observable) {
                 observable <- decide <= PsiAA
             } else {
                 observable <- decide <= PsiaA
             }
             if (observable) {
                 if (runif(1, 0, 1) <= p) y[q, j] <- "A"
             }
         } else {
             if (observable) {
                 if (runif(1, 0, 1) <= c) y[q, j] <- "A"
             }
         }
     }
}

# # # Robust design with markovian
# # # emigration
J <- 24
N <- 1000
S <- 0.9
PsiOO <- 0.4
PsiUO <- 0.3
p <- 0.5
c <- p
y <- matrix(0, N, J)
y[, 1] <- 1
for (q in 1:N) {
     (observable <- TRUE)
     for (j in 2:J) {
         if (j %% 2) {
             surviv <- runif(1, 0, 1)
             if (surviv <= S) {
                 decide <- runif(1, 0, 1)
                 if (observable) {
                   observable <- (decide <= PsiOO)
                 } else {
                   observable <- (decide <= PsiUO)
                 }
                 if (observable) {
                   y[q, j] <- if (runif(1, 0, 1) <= p) 1 else 0
                 }
             } else {
                 break
             }
         } else {
             if (observable) {
                 y[q, j] <- if (runif(1, 0, 1) <= c) 1 else 0
             }
         }
     }
}


Hope this helps,

Rui Barradas

> There are a lot of unnecessary tests and conditions. I tried to break
> down the code and write the tests that have been done when assigning a
> variable. I simplified your the first part but cannot guarantee that all
> criteria are met.
>
> #####COMMENTS#####
> for(j in which(dtp)){
>       for (q in 1:N){
>           (observable=1) #prefer TRUE/FALSE for boolean operators
>           if(j %% 2)
>           {
>               survive=runif(1,0,1)
>               if(survive<=S)
>               {
>                   stay=runif(1,0,1)
>               ####Is observable ?###
>                   if (observable==1)  #if (observable); but can
> observable!=1 here ???? It is defined as observable=1 above
>                   {
>                       if(stay<=PsiAA)
>                       {
>                           observable=1 #stay<=PsiAA && observable &&
> survive <= S && j%%2
>                           #not needed; observable is already equals to 1
>                       }else{
>                           observable=0  #stay>PsiAA && observable &&
> survive <= S && j%%2
>                           #ie. observable = !observable; ie. was 1, now is 0
>                       }
>                   }else{
>                       return=runif(1,0,1)
>                       if(return<=PsiaA)
>                       {
>                           observable=1 #return<=PsiAA && !observable &&
> survive <= S && j%%2
>                           #ie. observable = !observable; ie. was 0, now is 1
>                       }else{
>                           observable=0 #return>PsiAA && !observable &&
> survive <= S && j%%2
>                           #not needed; observable is already equals to 0
>                       }
>                   }
>               #######
>                   if(observable==1) #you could move this block above,
> where you assign the value to observable
>                   {
>                       capture=runif(1,0,1)
>                       if(capture<=p)
>                           {y[q,j]="A"} #capture<=p && observable &&
> survive <= S && j%%2
>                   }
>               #######
>               }else{
>                   deado=runif(1,0,1)
>                   if(deado<=PsiAd)
>                   {
>                       y[q,j]="d" #deado<=PsiAd && survive > S && j%%2
>                   }else{
>                       (y[q,j]="D") #deado<PsiAd && survive > S && j%%2
>                   } #use ifelse instead
>               #######
>                   if(y[q,j]%in%c("D","d")) # test not needed; y[q,j] will
> always be either D or d according the assignment above
>                   {
>                       break  #survive > S && j%%2
>                       #the use of break is not recommended, especially in
> this case with so many loops - it's hard to tell where the break will exit
>                   }
>               #######
>               }
>           }else{
>               if(observable==1)
>               {
>                   recap=runif(1,0,1)
>                   if(recap<=c)
>                   {
>                       y[q,j]="A"#recap<=c && observable && !(j%%2)
>                   }
>               }
>           }
>       }
> }
>
>
> ####SUGGESTED CODE######
>
> dtp <- rep(c(TRUE, TRUE, FALSE, FALSE), 6)[seq_len(J)]
> for(j in which(dtp)){
>       for (q in 1:N){
>           observable <-TRUE
>           if(j %% 2)
>           {
>               survive=runif(1,0,1)
>               if(survive<=S)
>               {
>                   stay=runif(1,0,1)
>                   return=runif(1,0,1)
>               ####Is observable ?###
>                   if (return<=PsiAA || stay>PsiAA)
>                   {
>                       observable = !observable
>                   }        #otherwise, observable stays identical
>               #######
>                   if(observable)
>                   {
>                       capture=runif(1,0,1)
>                       if(capture<=p)
>                           {y[q,j]="A"} #capture<=p && observable &&
> survive <= S && j%%2
>                   }
>               #######
>               }else{
>                   deado=runif(1,0,1)
>                   y[q,j]=ifelse(deado<=PsiAd, "d", "D") #deado<=PsiAd &&
> survive > S && j%%2
>                   break
>               #######
>               }
>           }else{
>               recap=runif(1,0,1)
>               if(recap<=c && observable)
>               {y[q,j]="A"}#recap<=c && observable && !(j%%2)
>           }
>       }
> }
>
> Regards,
>
> Eloi
>
> On 12-08-01 09:47 AM, Claudia Penaloza wrote:
>> I will accept any help you can give me, especially to free myself of
>>> SAS-brain inefficiencies...
>> My supervisor knows SAS not R, which may explain my code.
>> What I'm actually trying to do is simulate mark-recapture data with a
>> peculiar structure.
>> It is a multistate robust design model with dummy time periods... it
>> will basically be a matrix with a succession of letters (for different
>> states/age classes) and zeros that are generated following a certain
>> set of conditions (probability statements; drawing from a random
>> uniform distribution if an event happens).
>> Normally, survival and transition to another state occur between
>> primary sampling periods (in my R example, every two columns.. between
>> [,2] and [,3]) but in the "dummy time period" model survival occurs
>> first and then transition to another state, which is why I need my
>> "conditions" to alternate. Additionally, the robust design has
>> secondary sampling sessions that are within the same year, i.e.,
>> survival is assumed = 1, which is why my columns are paired. Each
>> state can also be in an unobservable state at any given time... all of
>> these details complicate the coding.
>> Below I've pasted what I have thus far... I have not debugged it yet
>> (the second loop isn't working yet and the first loop still has some
>> glitches).
>> Further below is properly working code for a robust design without
>> dummy time periods (it also lacks the dead states the dummy time
>> period model has, which happen to be part of the glitchy-ness).
>> Is there a better (more R-ish) way of doing this?
>> Thank you for all your help,
>> Claudia
>> ################################################################
>> # Robust design with Markovian emigration and dummy time periods
>> ################################################################
>> J = 24
>> N = 10
>> S = .9
>> PsiAD = 0.06
>> PsiAd = 0.04
>> PsiAA = .4
>> PsiaA = .3
>> p = .5
>> c = p
>> y <- matrix(0,N,J)
>> y[,1]="A"
>> dtp <- rep(c(TRUE, TRUE, FALSE, FALSE), 6)[seq_len(J)]
>> for(j in which(dtp)){
>>    for (q in 1:N){
>>      (observable=1)
>>      if(j %% 2){survive=runif(1,0,1)
>>                 if(survive<=S){
>>                   stay=runif(1,0,1)
>>                   if (observable==1){
>>                     if(stay<=PsiAA){
>>                       observable=1
>>                     }else{observable=0}
>>                     }else{
>>                       return=runif(1,0,1)
>>                       if(return<=PsiaA){
>>                         observable=1
>>                       }else{observable=0}}
>>                   if(observable==1){
>>                     capture=runif(1,0,1)
>>                     if(capture<=p){
>>                       y[q,j]="A"}
>>                   }}else{
>>                     deado=runif(1,0,1)
>>                     if(deado<=PsiAd){
>>                       y[q,j]="d"
>>                     }else(y[q,j]="D")
>>                       if(y[q,j]%in%c("D","d")){
>>                         break
>>                        }
>>                       }
>>                     }else{
>>                       if(observable==1){
>>                         recap=runif(1,0,1)
>>                         if(recap<=c){
>>                           y[q,j]="A"
>>                         }
>>                       }
>>                     }
>>                   }
>>                 }
>>
>> for(j in which(!dtp)){
>>    for (q in 1:N){
>>      if(j %% 2){
>>        stay=runif(1,0,1)
>>        if(observable==1){
>>          if(stay<=PsiAA){
>>            observable=1
>>          }else{observable=0}
>>        }else{
>>          return=runif(1,0,1)
>>          if(return<=PsiaA){
>>            observable=1
>>          }else{observable=0}
>>        }
>>        if(observable==1){
>>          capture=runif(1,0,1)
>>          if(capture<=p){
>>            y[q,j]="A"}
>>          }}else {
>>            if(observable==1){
>>              recap=runif(1,0,1)
>>              if(recap<=c){
>>                y[q,j]="A"
>>              }
>>            }
>>          }
>>        }
>>      }
>> ###########################################
>> ### Robust design with markovian emigration
>> ###########################################
>> J = 24
>> N = 1000
>> S = .9
>> PsiOO = .4
>> PsiUO = .3
>> p = .5
>> c = p
>> y <- matrix(0,N,J)
>> y[,1]=1
>> for (q in 1:N){
>>    (observable=1)
>>    for(j in 2:J){
>>      if(j %% 2){surviv=runif(1,0,1)
>>                 if(surviv<=S){
>>                   stay=runif(1,0,1)
>>                   if(observable==1){
>>                     if(stay<=PsiOO){
>>                     observable=1
>>                   }else{observable=0}
>>                     }else{
>>                       return=runif(1,0,1)
>>                       if(return<=PsiUO){
>>                         observable=1
>>                       }else{observable=0}}
>>                   if(observable==1){
>>                     cap=runif(1,0,1)
>>                     if(cap<=p){
>>                       y[q,j]=1}
>>                   }else y[q,j]=0
>>                 }else{break}
>>                 }else{
>>                   if (observable==1){
>>                     recap=runif(1,0,1)
>>                     if (recap<=c){
>>                       y[q,j]=1}
>>                     else{y[q,j]=0}
>>                     }
>>                   }
>>      }
>> }
>> On Tue, Jul 31, 2012 at 7:16 PM, David Winsemius
>> <dwinsemius at comcast.net <mailto:dwinsemius at comcast.net>> wrote:
>>
>>      Claudia;
>>
>>      The loop constructs will keep you mired in SAS-brain
>>      inefficiencies forever. Please, listen more carefully to Mercier.
>>
>>      --
>>      David
>>
>>
>>
>>      On Jul 31, 2012, at 3:54 PM, Mercier Eloi wrote:
>>
>>          On 12-07-31 02:38 PM, Claudia Penaloza wrote:
>>
>>              Thank you very much Rui (and Eloi for your input)... this
>>              is (the very
>>              simplified version of) what I will end up using:
>>
>>          Could we get the extended version ? Because right now, I don't
>>          know why
>>          you need such complicated code that can be done in 1 line.
>>
>>              J <- 10
>>              N <- 10
>>              y <- matrix(0,N,J)
>>              cols <- rep(c(TRUE, TRUE, FALSE, FALSE), 3)[seq_len(J)]
>>              for(j in which(cols)){
>>               for (q in 1:N){
>>                 if(j %% 2){
>>                   y[q,j]=1
>>                   }else{y[q,j]=2}
>>                 }
>>               }
>>              for(j in which(!cols)){
>>               for (q in 1:N){
>>                 if(j %% 2){
>>                   y[q,j]="A"
>>                   }else{y[q,j]="B"}
>>                 }
>>               }
>>
>>          There is no need for a double loop (on N) :
>>
>>          J <- 10
>>          N <- 10
>>          y <- matrix(0,N,J)
>>          cols <- rep(c(TRUE, TRUE, FALSE, FALSE), 3)[seq_len(J)]
>>          for(j in which(cols)){
>>             if(j %% 2){
>>                y[,j]=1
>>                }else{y[,j]=2}
>>            }
>>          for(j in which(!cols)){
>>              if(j %% 2){
>>                y[,j]="A"
>>                }else{y[,j]="B"}
>>            }
>>
>>          if you really wants to use this code.
>>
>>          Cheers,
>>
>>          Eloi
>>
>>              Cheers,
>>              Claudia
>>              On Mon, Jul 30, 2012 at 5:38 PM, Mercier Eloi
>>              <emercier at chibi.ubc.ca <mailto:emercier at chibi.ubc.ca>
>>              <mailto:emercier at chibi.ubc.ca
>>              <mailto:emercier at chibi.ubc.ca>>> wrote:
>>
>>                 Or, assuming you only have 4 different elements :
>>
>>                 mat<- matrix(rep(c(1,2,"A", "B"),each=10),10,10, byrow=F)
>>                 mat2 <- as.data.frame(mat)
>>
>>                 mat
>>
>>                       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
>>                  [1,] "1"  "2"  "A"  "B"  "1"  "2"  "A"  "B"  "1"  "2"
>>                  [2,] "1"  "2"  "A"  "B"  "1"  "2"  "A"  "B"  "1"  "2"
>>                  [3,] "1"  "2"  "A"  "B"  "1"  "2"  "A"  "B"  "1"  "2"
>>                  [4,] "1"  "2"  "A"  "B"  "1"  "2"  "A"  "B"  "1"  "2"
>>                  [5,] "1"  "2"  "A"  "B"  "1"  "2"  "A"  "B"  "1"  "2"
>>                  [6,] "1"  "2"  "A"  "B"  "1"  "2"  "A"  "B"  "1"  "2"
>>                  [7,] "1"  "2"  "A"  "B"  "1"  "2"  "A"  "B"  "1"  "2"
>>                  [8,] "1"  "2"  "A"  "B"  "1"  "2"  "A"  "B"  "1"  "2"
>>                  [9,] "1"  "2"  "A"  "B"  "1"  "2"  "A"  "B"  "1"  "2"
>>                 [10,] "1"  "2"  "A"  "B"  "1"  "2"  "A"  "B"  "1"  "2"
>>
>>                 mat2
>>                    V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
>>
>>                 1   1  2  A  B  1  2  A  B  1   2
>>                 2   1  2  A  B  1  2  A  B  1   2
>>                 3   1  2  A  B  1  2  A  B  1   2
>>                 4   1  2  A  B  1  2  A  B  1   2
>>                 5   1  2  A  B  1  2  A  B  1   2
>>                 6   1  2  A  B  1  2  A  B  1   2
>>                 7   1  2  A  B  1  2  A  B  1   2
>>                 8   1  2  A  B  1  2  A  B  1   2
>>                 9   1  2  A  B  1  2  A  B  1   2
>>                 10  1  2  A  B  1  2  A  B  1   2
>>
>>                 Cheers,
>>
>>                 Eloi
>>
>>
>>                 On 12-07-30 04:28 PM, Rui Barradas wrote:
>>
>>                     Hello,
>>
>>                     Maybe something along the lines of
>>
>>                     J <- 10
>>                     cols <- rep(c(TRUE, TRUE, FALSE, FALSE), 3)[seq_len(J)]
>>                     for(i in which(cols)) { do something }
>>                     for(i in which(!cols)) { do something else }
>>
>>                     Hope this helps,
>>
>>                     Rui Barradas
>>
>>                     Em 31-07-2012 00 <tel:31-07-2012%2000>
>>              <tel:31-07-2012%2000>:18, Claudia Penaloza
>>
>>                     escreveu:
>>
>>                         Dear All,
>>                         I would like to apply two different "for loops"
>>              to each
>>                         set of four columns
>>                         of a matrix (the loops here are simplifications
>>              of the
>>                         actual loops I will
>>                         be running which involve multiple if/else
>>              statements).
>>                         I don't know how to "alternate" between the loops
>>                         depending on which column
>>                         is "running through the loop" at the time.
>>                         ## Set up matrix
>>                         J <- 10
>>                         N <- 10
>>                         y <- matrix(0,N,J)
>>                         ## Apply this loop to the first two of every
>>              four columns
>>                         ([,1:2],
>>                         [,5:6],[,9:10], etc.)
>>                         for (q in 1:N){
>>                         for(j in 1:J){
>>                         if(j %% 2){
>>                         y[q,j]=1
>>                         }else{y[q,j]=2}
>>                         }
>>                         }
>>                         ## Apply this loop to the next two of every
>>              four columns
>>                         ([,3:4],[,7:8],[,11:12], etc.)
>>                         for (q in 1:N){
>>                         for(j in 1:J){
>>                         if(j %% 2){
>>                         y[q,j]="A"
>>                         }else{y[q,j]="B"}
>>                         }
>>                         }
>>                         I want to get something like this (preferably
>>              without the
>>                         quotes):
>>
>>                             y
>>
>>                         [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
>>                           [1,] "1"  "2"  "A"  "B"  "1"  "2"  "A"  "B"
>>               "1"  "2"
>>                           [2,] "1"  "2"  "A"  "B"  "1"  "2"  "A"  "B"
>>               "1"  "2"
>>                           [3,] "1"  "2"  "A"  "B"  "1"  "2"  "A"  "B"
>>               "1"  "2"
>>                           [4,] "1"  "2"  "A"  "B"  "1"  "2"  "A"  "B"
>>               "1"  "2"
>>                           [5,] "1"  "2"  "A"  "B"  "1"  "2"  "A"  "B"
>>               "1"  "2"
>>                           [6,] "1"  "2"  "A"  "B"  "1"  "2"  "A"  "B"
>>               "1"  "2"
>>                           [7,] "1"  "2"  "A"  "B"  "1"  "2"  "A"  "B"
>>               "1"  "2"
>>                           [8,] "1"  "2"  "A"  "B"  "1"  "2"  "A"  "B"
>>               "1"  "2"
>>                           [9,] "1"  "2"  "A"  "B"  "1"  "2"  "A"  "B"
>>               "1"  "2"
>>                         [10,] "1"  "2"  "A"  "B"  "1"  "2"  "A"  "B"
>>               "1"  "2"
>>
>>                         Any help greatly appreciated!
>>
>>                         Claudia
>>
>>
>>                 --
>>                 Eloi Mercier
>>                 Bioinformatics PhD Student, UBC
>>                 Paul Pavlidis Lab
>>                 2185 East Mall
>>                 University of British Columbia
>>                 Vancouver BC V6T1Z4
>>
>>
>>
>>
>>      David Winsemius, MD
>>      Alameda, CA, USA
>>
>>
>



More information about the R-help mailing list