Tracking the intersection point of simulated time series with a specific value over many runs in OpenBUGS
I have an OpenBUGS model that uses observable data (y.values) over time (x.values) to simulate many runs (~ 100000) with new y-value estimates (y.est) for each run. The observed data show a significant decrease from the maximum value.
I want to track how long it takes for each run to drop the maximum amount (T.max) to 10% of the maximum amount (T.10%). Since the maximum abundance value varies from run to run, 10% of this maximum will also vary from run to run, and therefore T.10% will vary from run to run.
Setting the parameter to store T.max is straightforward enough that it does not depend on the run to run, because the maximum value is large enough than any other value.
I can't figure out how to keep the intersection of y-est and T.10% values.
My first attempt was to determine if each y-est value is above or below T.10% using the function step()
:
above.below[i] <- step(T.10% - y.est[i])
This generates a string of them and zeros for each value of y.est (e.g. 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, etc.). If each run simply deviates continuously from maximum to minimum, I could use a function rank()
to determine how many values above.below[i]
occur above T. 10%:
decline.length <- rank(above.below[1:N], 0)
In this example, it decline.length
will equal the number '0 in the above line, which is 9. Unfortunately, the y-est values sometimes represent periods of growth after they have declined below T.10%. So a vector of values above.below
might look like this: 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 , 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, etc. Thus, it decline.length
will equal 14, not 9, given the subsequent 0s in the vector.
What I want to do is figure out how to store only the number "0" in above.below [1:10] and not above.below[1:N]
. Unfortunately, this is not always the 10th step where the first "1" occurs, so I have to do the maximum range value will above.below
vary from run to run during simulation.
I'm trying to accomplish this in OpenBUGS as it is a non-procedural language, but I think it can be done, I just don't know how. I hope that someone more familiar with the features step()
and rank()
can provide some expert guidance.
Any guidance is greatly appreciated!
source to share
Two solutions I suggest:
1) Calculate the total amount up to each time step:
for (i in 1:N){
above.below[i] <- step(T.10% - y.est[i])
cum.above.below[i] <- sum(above.below[1:i])
}
decline.length <- rank(cum.above.below[1:N], 0)
2) Calculate whether each year is above or below the threshold directly, without 1s and 0s:
for(i in 1:N){
above.below[i] <- step(T.10% - y.est[i])
dummy[i] <- above.below[i] * i + (1 - above.below[i]) * (N+1)
}
decline.length <- ranked(dummy[], 1)
So dummy
is i
when above.below
equal to 1, and dummy
equal to N + 1 when above.below
equal to 0.
source to share