Skip to content

Alternating between "for loops"

8 messages · Mercier Eloi, David Winsemius, Claudia Penaloza +1 more

Rui Barradas
# 
Hello,

Em 01-08-2012 20:02, Mercier Eloi escreveu:
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
David Winsemius
#
On Aug 1, 2012, at 12:02 PM, Mercier Eloi wrote:

            

      
      
      
    

        
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.

            

      
      
      
    

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

            

      
      
      
    

        
# 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?

# Not sure I have gotten the model assumptions down.
# Is D" or "d" an absorbing state such as "dead or "Dead"?

            

      
      
      
    

        
#Presumably the "dummy time periods"

            

      
      
      
    

        
# It would help if the concept of "observable" were explained
# Is this something to do with the capture-recapture observables?

            

      
      
      
    

        
# 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.)

                     observable <- (stay<=PsiAA)

#--------------

            

      
      
      
    

        
# better NOT to use the name "return" since it is a crucial R function  
name.

            

      
      
      
    

        
#------- perhaps:

              randret <-  return=runif(N,0,1)
              observable <- randret <= PsiaA
# -------------------

            

      
      
      
    

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

     #------------

            

      
      
      
    

        
# -------Perhaps
      deado <- runif(N, 0,1)
      y[ , j] <- ifelse( deado<=PsiAd, "d", "D")
#------------------------

            

      
      
      
    

        
#    ----------Really?  I thought that condition was assured at this  
point?

            

      
      
      
    

        

      
      
    

        
David Winsemius, MD
Alameda, CA, USA

  
  


  


      

    

    
1 day later

    

      
      

        
      

    

    
5 days later

    

      
      

        
      

    

    
1 day later

    

      
      

        
      

    

    

      
      

        


  

        

      Rui Barradas
    

    

      
      #
    

  


  

      
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:

            

      
      
      
    

  
  


  


      

    

    
2 days later