I don't really have another place to write up, so here is a computational problem I keep coming back to: getting individual probs of survival from an array of times and an array of probs
Let me sharpen the problem, I have 2 arrays for probability of survival and the corresponding time. Everything following will be written in Python
>>> pr = [0.9, 0.8, 0.7, 0.5, 0.2, 0.1]
>>> t = [1, 2, 3, 5, 6, 9]
So S(t=1)=0.9, S(t=2)=0.8, and so on.
Given an array of individual times
>>> ti = [0.9, 5, 7.5, 6, 9, 10]
we want to get the probability of survival at that time. So our output should look like
>>> [1.0, 0.5, 0.2, 0.2, 0.1, 0.1]
The first element is 1 because S(0)=1 by definition. So, how can we get this array?
My first thought is that we can write a loop. We look over each element in `ti`, then use that element to find the last valid index in `t` (for ti=7.5 that would be 6), then use that index to look up the element in `pr`.
It might be faster to convert `pr` and `t` to a dictionary (in Python) but this is still expensive as we need to loop. If there are lots of elements in `ti`, this can become far too slow. So the idea is to vectorize this procedure. The following is the best solution I’ve come up with at this point
So what does this do. Setup is done on lines 1-3. Lines 1-2 tack on S(0)=1 to the start of the `t` and `pr` arrays. Line 3 adds the maximum time we saw in `ti` then drops the first element (zero)
The clever part (in my opinion) is the remainder. Line 4 creates a matrix of indicators where the columns are the individuals and rows are whether their value in `ti` was >= the corresponding `t`. I do a similar thing in line 5 but now < the shifted times. When these are multiplied together in 6, we get a matrix where the only non-zero value in a column corresponds to the final time the person was observed
Then in the final step, we multiply that by `pr`. Therefore, we can sum over the rows, as the only non-zero row in this matrix will be the probability at their corresponding individual time. There is probably a better way to do this, but I think it’s a clever vectorization of the problem
The censoring events before the first event time in `t`, and censoring times greater than the last event time in `t` were the annoying parts to deal with. The first 3 lines were primarily to setup the input arrays to deal with those issues
So first, I am going to convert `t`, `ti`, `pr` to NumPy arrays so we can use NumPy to help speed things along. After that
>>> t = np.insert(t, 0, 0)
>>> pr = np.insert(pr, 0, 1)
>>> shift_t = np.append(t, np.max(ti)+1)[1:]
>>> upper = (ti >= t[:, None]).astype(int)
>>> lower = (ti < shift_t[:, None]).astype(int)
>>> t_ind = upper * lower
>>> np.sum(t_ind * pr[:, None], axis=0)