[R] Alternating between "for loops"

Rui Barradas ruipbarradas at sapo.pt
Sat Aug 11 20:15:55 CEST 2012


Hello,

I'm not sure it works, but try the following.

for(j in which(dtp)){
   for (q in 1:N){
     if(y[q, j] %in% c("d", "D")) break
     [...etc...]

and in the other loop the same,

for (j in which(!dtp)) {
   for (q in 1:N) {
     if(y[q, j] %in% c("d", "D")) break
     [...etc...]


Em 10-08-2012 20:42, Claudia Penaloza escreveu:
> This is what my code looks like now. However, there is one thing I
> can't/don't know how to fix.
> I can't get it to be "once dead always dead", i.e., in any given row, after
> a "D" or a "d" there should be only zeros.
> I've tried applying a flag to break the loop if dead but I can't get it to
> work.
> Could you please help me with this?
>
> Thank you for your time,
> Claudia
>
> J = 24
> N = 10
> S = .9
> PsiADd = 0.4
> 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 <- TRUE)
>      if(j %% 2){
>        if (runif(1,0,1)<=S) {
>          if (observable) {
>            observable <- (runif(1,0,1)<=PsiAA)
>          } else {
>            observable <- (runif(1,0,1)<=PsiaA)
>          }
>          if (observable) {
>            if (runif(1,0,1) <= p) y[q,j] <- "A"
>          }
>        } else {
>          y[q,j] <- ifelse(runif(1,0,1) <= PsiADd,"d","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) {
>        if (observable) {
>          observable <- runif(1, 0, 1) <= PsiAA
>        } else {
>          observable <- runif(1, 0, 1) <= 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"
>        }
>      }
>    }
> }
> On Wed, Aug 8, 2012 at 2:04 PM, Claudia Penaloza
> <claudiapenaloza at gmail.com>wrote:
>
>> Answers inserted in BLUE below
>>
>> On Thu, Aug 2, 2012 at 5:34 PM, Claudia Penaloza <
>> claudiapenaloza at gmail.com> wrote:
>>
>>> Thank you very much for all your suggestions. I am very sorry my code is
>>> so crude (it gives me a headache too!), I'm very new to R programing. I
>>> will have to look at your suggestions/questions very carefully (I'm barely
>>> holding on at this point) and get back to you (Dr. Winsemius) with some
>>> answers.
>>>
>>> Thank you!
>>> Claudia
>>>
>>> On Wed, Aug 1, 2012 at 2:11 PM, David Winsemius <dwinsemius at comcast.net>wrote:
>>>
>>>> On Aug 1, 2012, at 12:02 PM, Mercier Eloi wrote:
>>>>
>>>>   Your code almost gave me a headache. :-/
>>>> I had a similar reaction. However, my approach might have been to
>>>> request a more complete verbal description of the data structures being
>>>> operated on and the methods and assumptions being used. Generally it is
>>>> going to be easier to go from a procedural description to good R code than
>>>> it is to go from a SAS Data Proc ... even if it were tested and debugged...
>>>> and yours was not even debugged.
>>>>
>>>>
>>>>   ##############################**##############################**####
>>>>> # 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)
>>>>>
>>>> # So is 'y' a matrix of states where the N rows are levels and the J
>>>> columns are discrete times?
>>>> # Actually I decided that context suggested you were using letters to
>>>> represent states.
>>>>
>>> Yes, 'y' is a matrix of state, N rows are levels (independent of each
>> other) and J columns are discrete times.
>> Yes, I am using letters to represent states.
>>
>>>>   y[,1]="A"
>>>> # So we start with N <something>'s in state "A"?
>>>> # It seems as though it might be the case that every row is independent
>>>> of the others.
>>>> # .. and you are performing ( replicate()-ing in R terminology) this
>>>> test N times
>>>> # states:
>>>> #("A" "D")
>>>> #("a" "d")
>>>> #transitions:
>>>>     matrix( c( runif(1, 0, 1) <= 0.4, # A -> A
>>>>            runif(1,0,1) <= 0.3,  # a -> A
>>>>            runif(1,0,1) <= 0.06. # A -> D
>>>>            runif(1,0,1) <= 0.04  # A -> d) )
>>>>
>>>> # What is state "a"?
>>>> # How do you get from A to a?
>>>> # can you get from D or d to A or a?
>>>>
>>> Yes, we start with N "individuals" in state "A", each individual is
>> independent of each other.
>> State "a" is an unobserved (!observable) version of state "A"
>> (biologically, an individual that has temporarily left the sampling area
>> and is not available for capture)
>>   An individual in state "A" transitions to state "a" if it is observable
>> and doesn't stay (stay >= AA)
>> "D" and "d" are dead (absorbing) states. Once in either "D" or "d", the
>> individual can no longer transition and is no longer "captured" (the row
>> should only have zeros after a "D" or a "d")
>>
>>
>>>
>>>
>>>> # Not sure I have gotten the model assumptions down.
>>>> # Is D" or "d" an absorbing state such as "dead or "Dead"?
>>>>
>>> Yes.
>>> dtp <- rep(c(TRUE, TRUE, FALSE, FALSE), 6)[seq_len(J)]
>>>>
>>>> #Presumably the "dummy time periods"
>>>>
>>> Yes.
>>
>>>> for(j in which(dtp)){   # j will be c(1,2, 5,6, 9,10 , ..., 21,22)
>>>>
>>>>
>>>>
>>>>   for (q in 1:N){  # This can almost surely become vectorized
>>>>
>>>>     (observable=1)
>>>>     if(j %% 2){survive=runif(1,0,1)
>>>>                if(survive<=S){
>>>>                  stay=runif(1,0,1)
>>>>                  if (observable==1){
>>>>
>>>>
>>>> # It would help if the concept of "observable" were explained
>>>> # Is this something to do with the capture-recapture observables?
>>>>
>> Yes. All individuals start out "observable" (we need to capture them a
>> first time). Individuals can then stay in their current state or transition
>> to another one, if they are observable they can continue to be observable,
>> or they can become unobservable and viceversa. Transition depends on what
>> state you are in at any given time (A->A != A->a != a->A != a->a).
>>
>>>>                    if{
>>>>                      observable=1
>>>>                    }else{observable=0}
>>>>                    }else{
>>>>
>>>>
>>>> # After allowing for Mercier's astute observation that observable will
>>>> always be 1, I'm wondering if it can be replaced with just this code (and
>>>> not need the loop surrounding it.)
>>>>
>>>>
>>>>
>>> The individual will always "start" as observable... but as the loop
>> progresses through the columns in a row, it can go between being observable
>> and unobservable (at least that is what I wanted to code).
>>
>>>>                      observable <- (stay<=PsiAA)
>>>>
>>>> #--------------
>>>>
>>>>                      return=runif(1,0,1)
>>>>
>>>>
>>>> # better NOT to use the name "return" since it is a crucial R function
>>>> name.
>>>>
>>>>
>>>>
>>> Will change. Thank you.
>>>>
>>>>
>>>>                      if(return<=PsiaA){
>>>>                        observable=1
>>>>                      }else{observable=0}}
>>>>
>>>>
>>>> #------- perhaps:
>>>>
>>>>               randret <-  return=runif(N,0,1)
>>>>               observable <- randret <= PsiaA
>>>> # -------------------
>>>>
>>>>
>>>>                  if(observable==1){
>>>>                    capture=runif(1,0,1)
>>>>                    if(capture<=p){
>>>>
>>>>
>>>>      #---------- perhaps:
>>>>            randcap <- runif(N, 0,1)
>>>>            y[ randcap <= p, j] <- "A"
>>>> #That would replace with "A" (within the j-column) ...
>>>> #     only the rows in 'y' for which randcap were less than the randcap
>>>> threshold.
>>>>
>>>>      #------------
>>>>
>>>>
>>>>                      y[q,j]="A"}
>>>>                  }}else{
>>>>
>>>>
>>>>
>>>>
>>>>                    deado=runif(1,0,1)
>>>>                    if(deado<=PsiAd){
>>>>                      y[q,j]="d"
>>>>                    }else(y[q,j]="D")
>>>>
>>>>
>>>> # -------Perhaps
>>>>       deado <- runif(N, 0,1)
>>>>       y[ , j] <- ifelse( deado<=PsiAd, "d", "D")
>>>> #------------------------
>>>>
>>>>
>>>>
>>>>                      if(y[q,j]%in%c("D","d")){
>>>>                        break
>>>>
>>>>
>>>> #    ----------Really?  I thought that condition was assured at this
>>>> point?
>>>>
>>>>
>>>>
>>> I thought so too but non-zero elements appeared toward the right in rows
>> that had already "died".
>>
>>>>                       }
>>>>                      }
>>>>                    }else{
>>>>                      if(observable==1){
>>>>                        recap=runif(1,0,1)
>>>>                        if(recap<=c){
>>>>                          y[q,j]="A"
>>>>                        }
>>>>                      }
>>>>                    }
>>>>                  }
>>>>                }
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> 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.
>>>>
>>>>
>>>>
>>> I will run the three modified versions and see how things go. I hope my
>> answers are helpful.
>>
>>
>>



More information about the R-help mailing list