Your code almost gave me a headache. :-/
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:
> 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
> <dwinsem...@comcast.net <mailto:dwinsem...@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
> <emerc...@chibi.ubc.ca <mailto:emerc...@chibi.ubc.ca>
> <mailto:emerc...@chibi.ubc.ca
> <mailto:emerc...@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@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.
