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
<dwinsem...@comcast.net <mailto:dwinsem...@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
             <emerc...@chibi.ubc.ca <mailto:emerc...@chibi.ubc.ca>
             <mailto:emerc...@chibi.ubc.ca
             <mailto:emerc...@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




______________________________________________
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.

Reply via email to