[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