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.
#####COMMENTS#####
for(j in which(dtp)){
for (q in 1:N){
(observable=1) #prefer TRUE/FALSE for boolean operators
if(j %% 2)
{
survive=runif(1,0,1)
if(survive<=S)
{
stay=runif(1,0,1)
####Is observable ?###
if (observable==1) #if (observable); but can
observable!=1 here ???? It is defined as observable=1 above
{
if(stay<=PsiAA)
{
observable=1 #stay<=PsiAA && observable &&
survive <= S && j%%2
#not needed; observable is already equals to 1
}else{
observable=0 #stay>PsiAA && observable &&
survive <= S && j%%2
#ie. observable = !observable; ie. was 1, now is 0
}
}else{
return=runif(1,0,1)
if(return<=PsiaA)
{
observable=1 #return<=PsiAA && !observable &&
survive <= S && j%%2
#ie. observable = !observable; ie. was 0, now is 1
}else{
observable=0 #return>PsiAA && !observable &&
survive <= S && j%%2
#not needed; observable is already equals to 0
}
}
#######
if(observable==1) #you could move this block above,
where you assign the value to observable
{
capture=runif(1,0,1)
if(capture<=p)
{y[q,j]="A"} #capture<=p && observable &&
survive <= S && j%%2
}
#######
}else{
deado=runif(1,0,1)
if(deado<=PsiAd)
{
y[q,j]="d" #deado<=PsiAd && survive > S && j%%2
}else{
(y[q,j]="D") #deado<PsiAd && survive > S && j%%2
} #use ifelse instead
#######
if(y[q,j]%in%c("D","d")) # test not needed; y[q,j] will
always be either D or d according the assignment above
{
break #survive > S && j%%2
#the use of break is not recommended, especially in
this case with so many loops - it's hard to tell where the break will exit
}
#######
}
}else{
if(observable==1)
{
recap=runif(1,0,1)
if(recap<=c)
{
y[q,j]="A"#recap<=c && observable && !(j%%2)
}
}
}
}
}
####SUGGESTED CODE######
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)
if(survive<=S)
{
stay=runif(1,0,1)
return=runif(1,0,1)
####Is observable ?###
if (return<=PsiAA || stay>PsiAA)
{
observable = !observable
} #otherwise, observable stays identical
#######
if(observable)
{
capture=runif(1,0,1)
if(capture<=p)
{y[q,j]="A"} #capture<=p && observable &&
survive <= S && j%%2
}
#######
}else{
deado=runif(1,0,1)
y[q,j]=ifelse(deado<=PsiAd, "d", "D") #deado<=PsiAd &&
survive > S && j%%2
break
#######
}
}else{
recap=runif(1,0,1)
if(recap<=c && observable)
{y[q,j]="A"}#recap<=c && observable && !(j%%2)
}
}
}
Regards,
Eloi
On 12-08-01 09:47 AM, Claudia Penaloza wrote:
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