Inspired by Martin Morgan 's decision to a closely related issue, here is a really vectorized way to generate P_IDusing functions cummax. As soon as you notice that P_IDis closely linked to cumsumfrom !(r < 0.5):
set.seed(1)
N <- 10
r <- c(0, runif(N-1))
H_ID <- cumsum(r < .5)
r_ <- r >= .5 # flip the coins that generated H_ID.
z <- cumsum(r_) # this is almost P_ID; just need to subtract the right amount...
# ... and the right amount to subtract is obtained via cummax
P_ID <- 1 + z - cummax( z * (!r_) )
> cbind(H_ID, P_ID)
H_ID P_ID
[1,] 1 1
[2,] 1 2
[3,] 2 1
[4,] 3 1
[5,] 3 2
[6,] 3 3
[7,] 3 4
[8,] 4 1
[9,] 5 1
[10,] 5 2
I have not done detailed time tests, but it is probably evil, because these are all internal, vectorized functions
source
share