3

I have been through multiple R-while-if questions and I didn't find a similar one. Please guide me to the answer if it has been asked.

So I have a vector V of N=200 entries of number 0. Every entry signifies a "bin". I skip the first 3 bins, for separate reasons. Now for i= 3:N number of bins, I use a while loop.

  • For each ith number of bin, I generate a random number between 0 and 1.
  • If the random number is less than a certain numerical value (here it's 0.66) then, I simply replace the 0 at index i in Vector V with 1 and at the same time, I want while loop to skip next two values.

  • if the random number is greater than the numerical value, I simply consider the next i.

    Here's the code I am running: 

    N = 200
Time <- rep(0, N)
Time <- replace(Time, 3, 1)
i = 5
while (i <= N){
p <- runif(1)
 if(p < 0.66){
   i= i+ 2; replace(Time, i, 1)
  } else {
    i <- i+1
   }
 }

I am not very good at R, and this is the combination I have been trying to use, to get R to do what I want, but it's clearly not working. I might be missing something obvious.

halfer
  • 19,824
  • 17
  • 99
  • 186
14thTimeLord
  • 363
  • 1
  • 14

1 Answers1

4

Your problem is that most functions in R do not modify objects they are applied to. Instead, they return a modified version of the original object which you have to assign, either to a new object, or to the original one (update it).

So instead of just running replace(Time, i, 1) (which does not update the Time vector), you need to assign the modified version back to Time like this:

Time <- replace(Time, i, 1)

or, better yet, use square bracket notation to update directly:

Time[i] <- 1

Moreover, you need to update vector Time, Time <- replace(Time, i, 1), before you update your index i, i = i + 2. If you want to apply your loop only to i > 3, then you need to set i = 4. Finally, in case runif(1) < 0.66, if you want to skip the next two entries in Time, you should update i as i <- i + 3.

To summarise, you should modify your loop as follows:

N = 200
Time <- rep(0, N)
Time[3] <- 1
i <- 4

while (i <= N){
    p <- runif(1)
    if(p < 0.66){
        Time[i] <- 1
        i <- i + 3
    } else {
        i <- i+1
    }
}
stas g
  • 1,503
  • 2
  • 10
  • 20
  • 2
    Good answer! For a fully complete answer, you should strengthen your recommendation to use square bracket notation over functions like `replace` or `subset`: https://stackoverflow.com/questions/9860090/why-is-better-than-subset – divibisan Aug 06 '18 at 16:49
  • 2
    @divibisan `replace` is `[<-`, so I'm not sure that criticism applies. I guess a real improvement for the OP would come from using vectorized code instead of a loop... – Frank Aug 06 '18 at 17:07
  • 1
    @Frank yes, `subset` is not a relevant function in this question. but i don't think vectorising this would be trivial, seeing that each step depends on the previous one? – stas g Aug 06 '18 at 17:11
  • 1
    @stasg You'd get a different outcome for a given seed, but could achieve the same distribution over outcomes by assigning 0-1 based on `w = runif(200) < 0.6` and then afterwards applying the rules: `w[1:3] = c(FALSE, FALSE, TRUE); w[c(which(w)+1L, which(w)+2L)] = FALSE; as.integer(w)` or something. (Not sure about this code... also, there may be more direct ways.) – Frank Aug 06 '18 at 17:14
  • 2
    @Frank: It's hard to say the best way to vectorize this program since it's not really clear (at least to me) what OP wants to do and so it's hard to say how to do it better. If you have an idea, though, you should write your own answer! – divibisan Aug 06 '18 at 17:17
  • 1
    @divibisan Yeah, I agree. Hard to know if their real problem is necessarily iterative or not. This is a good answer for the points that we can be sure about. – Frank Aug 06 '18 at 17:19
  • 2
    @Frank Good point on `replace`, I didn't look at the source before. I still think it shouldn't be used for purposes of clarity, but you're right that there shouldn't be any risks to it. – divibisan Aug 06 '18 at 17:20
  • Thank you @stasg, this worked. Also thank you to everybody in the comments. I am simulating a Poisson Process with dead period. So the 1's in the vector V (V is a vector of time-bins, 200 time-bins for a 2 second spike train) are where the spike happens in the timeline. I needed a code to ensure that no spike happens for a certain time period after each spike. That's why we skip two time bins, every time there's a `1` entry. This was for one spike train. I will need to reiterate it multiple times to get more spike train. Please let me know if changing something would help the process. – 14thTimeLord Aug 07 '18 at 07:15