Consider the following problem. There is a gym with a turnstile that marks every time someone enters or exits the gym. At the beginning and end of the day the gym is empty: the number of entries and exits on each day is equal. The problem is to infer the distribution of gym-session lengths. This would be straightforward if we could match up entrances and exits, as we would then know how long each individual gym session lasted, but we do not have this information.
To make this more concrete, we have the following as data:
with the constraint that matching_exit is a permutation of 1..N. The duration of session i is
t_exit[matching_exit[i]] - t_enter[i]
If matching_exit were data we could create a Stan model, perhaps using some chosen parametric form and moderately informative prior, to infer the distribution of session times. But it’s a latent variable of a sort that Stan doesn’t support, and the usual tricks for marginalizing out discrete variables don’t seem to work here.
Is there some modeling trick I’ve missed that would allow Stan to handle a problem like this? If not, what tool would you recommend for problems like this where part of the problem is matching up data items?
My intuition is that with some assumptions about the arrival process and session time this can be reduced to something like an M/M/\infty or M/M/c queue problem.
I think you should be able to write the likelihood for the number of gym occupants at a given time N(t) in terms of arrival process/session length parameters based on the transient solution for the particular queue model.
The general idea should hold for more complex assumptions about the arrival process and session time distribution, the transient solution will just become more complicated. The c versus \infty assumption is based on if gym goers are always able to start their workout right away or if they sometimes have to wait.
That would require assuming an exponential distribution for gym-session duration, correct? It crucially depends on the exit rate being a function only of the number of gym occupants at any given time, and not of how long they have been there.
My queue theory is a bit rusty – I think the assumptions are typically made about about the arrival process and the service time distributions (possibly also the wait for service) and the departure process is derived as logical result of these two.
For the nice closed form N(t|\theta) in the M/M/\infty case, I think the session time has to be exponentially distributed. I believe there is also a general case where the session time is i.i.d. but not necessarily exponential where there are approximation strategies for N(t|\theta). For a more general case of non-iid session times where, say, the session time depends on the current number of occupants, I’m not sure.
What do you care about modeling here? If I was running a business like a gym or a restaurant, I’d want to know how many people to expect on any given day and the distribution with which they’d arrive so I can model things like variance of how crowded it is. We can’t get either of those from a model of gym-session lengths. I think we’d need to move to something like a Prophet-style generalized additive model for crowdedness by day so you could take into account seasonal effects, day-of-week, holiday effects and trends.
The total duration of stays and the average can be estimated directly.
average_stay = (sum(t_exit) - sum(t_enter)) / N
How much do you care about variance? Without assumptions on variance, we can’t calculate alignments.
For each entrance time, it’s easy to find all the potential exit times that could match it, but there’s not an easy way to select one and it’s even harder to keep this all consistent. Given that it’s a permutation problem, it may just be NP-hard to do the matchups in general even once you have the distribution of stays.
You can also lower-bound the durations by taking sort(t_exit) - sort(t_enter), e.g., [3, 6, 9] - [1, 2, 5] = [2, 4, 4]. The average duration length isn’t going to change no matter how they’re aligned, only the variance will. For example, I could take the order [9, 3, 6] - [1, 2, 5] = [8, 1, 1], which is higher variance, but the same average time.
The gym-session problem itself is not of particular importance to me; rather, it’s a simple representative of a class of problems that I have run into on multiple occasions, that involve this sort of matching of different items of data.
Identity stitching is another example of this sort of problem, in web analytics. You can set cookies in a user’s browser to use as an ID, and this works well for identifying different page hits that come from the same user session, but people change browsers or clear their cache and so on, which means their ID can frequently change, resulting in multiple apparently distinct users that are actually the same person. There are various sorts of clues (geolocation of IP address, behavior, etc.) that may suggest that two apparently distinct users are the same one, and you would like to use this information to probabilistically infer which distinct user IDs correspond to the same person.
The common thread here is that the obvious model involves a discrete latent variable that is a large vector assigning a label to each datum, where the number of distinct labels is very large and grows linearly with the size of the data set. In such problems you can’t marginalize out the discrete latent variable. You could perhaps write a custom MCMC sampler that includes updates that swap labels or something, but (1) that’s a lot of work, (2) getting it to sample efficiently looks difficult. (Still, it would be nice if Stan had a way to interleave such discrete updates of discrete parameters with the highly efficient HMC on the continuous parameters, so one could at least quickly try out such hybrid inference.)
Doing some more research, it looks like one possibility for the gym problem in particular might be particle MCMC. Not sure if that approach would work for something like the identity stitching problem, though.
I have a few ideas that might help you in developing a model for these problems. For the device tracking problem, when I worked at Comscore my colleagues developed an IP-colocation method coupled with a community detection algorithm to accurately identify devices in a household without persistent cookies or even persistent IP addresses. See https://dl.acm.org/doi/epdf/10.1145/3219819.3219852.
The graph is built by discretizing time where the edge strengths are determined by the inverse of the number of distinct device IDs on IP address at time t. This is described in detail in an earlier paper, https://pages.cs.wisc.edu/~pb/kdd17_final.pdf.
Once a majority of cookies are gone, due to browser policies, the graph collapses though at what point that happens was never researched.
If you have a subsample, such as a panel of users, you can use this information to estimate assignment in a probabilistic manner. First, estimate the panel subsample distribution of gym times by estimating or assigning one of k distributions to their gym times. Since many of the gym goers assignment is unknown you estimate a simplex over which the k distributions sum to the total gym time as given by the assigned cluster distribution for user i. You do assume some parametric distributions for the k clusters as well as the number of clusters a prior. This is written up in a Google research paper Virtual People: Actionable Reach Modeling.