[R] Alternating between "for loops"
David Winsemius
dwinsemius at comcast.net
Wed Aug 1 22:11:16 CEST 2012
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.
> 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?
# Not sure I have gotten the model assumptions down.
# Is D" or "d" an absorbing state such as "dead or "Dead"?
> dtp <- rep(c(TRUE, TRUE, FALSE, FALSE), 6)[seq_len(J)]
#Presumably the "dummy time periods"
> 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?
> 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.)
observable <- (stay<=PsiAA)
#--------------
> return=runif(1,0,1)
# better NOT to use the name "return" since it is a crucial R function
name.
> 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?
> }
> }
> }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.
>
> 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
>>
>>
>
>
> --
> Eloi Mercier
> Bioinformatics PhD Student, UBC
> Paul Pavlidis Lab
> 2185 East Mall
> University of British Columbia
> Vancouver BC V6T1Z4
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at 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.
David Winsemius, MD
Alameda, CA, USA
More information about the R-help
mailing list