Alternating between "for loops"
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