I have a discrete-time stochastic model where individuals enter in state 0, can step up/down in their state value, and then exit. The numbers of individuals and time-steps are very large so I only save the times of entry, state-transitions (up, down), and exit. From these logs, I want to reconstruct the mean state value among active individuals (after entry & before exit) for all time.
Problem / Question: The solution below is not very fast (especially for large ni
and tf
). Can we obtain val.t / act.t
more efficiently? Ideally avoiding the loop, but I’m not sure if that’s possible due to the ragged nature of to.i
and tx.i
(each individual experiences a random number of events). I know we can speed-up with parallel::
but I’m already using all CPUs with model instances in parallel.
MWE
# config
set.seed(1)
ni = 1000 # num individuals (i)
ne = 100 # mean num state change events per i
tf = 365*100 # total timesteps
# generate event times (dummy model)
t0.i = runif(ni,1,tf/2) # model entry time per i
tf.i = t0.i + tf/2 # model exit time per i
ne.i = rpois(ni,ne) + 1 # num events per i
to.i = lapply(1:ni,function(i){ round(runif(ne.i[i],t0.i[i],tf.i[i])) }) # state up times
tx.i = lapply(1:ni,function(i){ to.i[[i]] + rpois(ne.i[i],365) }) # state down times
# mean value among active
act.t = numeric(tf) # will hold total num active
val.t = numeric(tf) # will hold total state value among active
for (i in 1:ni){ # for each i
ti = t0.i[i]:tf.i[i] # timesteps active
val.i = cumsum(tabulate(to.i[[i]],tf)-tabulate(tx.i[[i]],tf)) # state value
act.t[ti] = act.t[ti] + 1 # add to active
val.t[ti] = val.t[ti] + val.i[ti] # add to state value
}
# plot
plot(val.t/act.t,type='l') # mean value among active
lines(act.t/ni,col='red') # proportion active
Plot (FYI)