Generalized priority-based model for smartphone screen touches

The distribution of intervals between human actions such as email posts or keyboard strokes demonstrates distinct properties at short versus long timescales. For instance, at long timescales, which are presumably controlled by complex process such as planning and decision making, it has been shown that those interevent intervals follow a scale-invariant (or power-law) distribution. In contrast, at shorter timescales—which are governed by different processes such as sensorimotor skill—they do not follow the same distribution and we know little about how they relate to the scale-invariant pattern. Here, we analyzed 9 million intervals between smartphone screen touches of 84 individuals which span several orders of magnitudes (from milliseconds to hours). To capture these intervals, we extend a priority-based generative model to smartphone touching events. At short timescale, the model is governed by refractory effects, while at longer timescales, the intertouch intervals are governed by the priority difference between smartphone tasks and other tasks. The ﬂexibility of the model allows us to capture interindividual variations at short and long timescales, while its tractability enables efﬁcient model ﬁtting. According to our model, each individual has a speciﬁc power-law exponent which is tightly related to the effective refractory time constant suggesting that motor processes which inﬂuence the fast actions are related to the higher cognitive processes governing the longer interevent intervals. DOI:


I. INTRODUCTION
Human actions such as mail correspondences, library loans, or website visits are not equally distributed in time but are typically structured in bursts followed by long periods of inactivity [1][2][3]. Several types of models have been proposed to capture the power-law structure of interevent time distribution (for a review see Ref. [3]). Priority-based queuing models [4][5][6][7][8][9] rely on one (or multiple) list(s) of tasks to be executed, where each task is associated with a priority level which directly influences the timing of its execution. This class of models have been pioneered by Barabási [4] and then generalized to multiple interacting queues [7,9], time-varying priorities [10], or priorities which depend on the position within the list of tasks [11]. Those models provide an interesting interpretation for the origin of the power-law scaling for long intervals (they come from prioritizing tasks) but are usually not designed to capture short interevent timings.
Poisson-based models belong to another class [12][13][14][15]. They rely on the assumption that the event rate is governed by a Poisson process whose rate can change over time. Because of this dynamic rate assumption, those Poisson-based models can easily accommodate a precise description of short interevent intervals. In the simplest case where the Poisson rates are piecewise constant and stochastically jump at each event time, the power-law exponent can be directly obtained from the distribution of Poisson rates [12]. This approach has been extended to continuously changing Poisson rates [12][13][14][15]. In particular, the framework proposed by Malmgren et al. [14] provides a circadian explanation for the origin of power-law distributions. Self-exciting point processes (also called Hawkes processes [16,17]) have also been used to provide a mechanistic interpretation of power law interevent intervals [18,19]. Those models are very flexible (they can accommodate both short and long interevent intervals), but lack the priority-related interpretation. Indeed, it is unclear how those Poisson-based models relate to priority-based models.
Here, we start from a priority-based framework, generalise it on different levels and apply it to smartphone touchscreen interaction data (see Fig. 1). First, unlike the priority-based model from Barabási, our priority-based model aims at predicting interevent distribution instead of response-time distribution. Second, our model is the continuous-time extension of the classical priority-based model. Under this limit, we can compute analytically the interevent distribution and show that our generalized priority-based model can be mapped to Poisson-based models. Third, our priority-based model does not only describe long interevent intervals but also includes a detailed description of short interevent intervals and thereby overcomes the need to define an arbitrary onset of the powerlaw distribution [20]. In particular, it assumes that the agent remains in a so-called refractory state during a short time after each event, where the probability of generating a new event is reduced.
Finally, because our model is based on arbitrary priority distribution and not on specific priority distribution imposed by the presence of lists (with discrete number of items), it can produce any power-law exponent. We found that for each subject, the intertouch interval (ITI) distribution is different and well captured by the model. We also found that from those fitted parameters, we can quantify the relative priority placed on smartphone actions.

A. Discrete-time model
In the first step, we propose a discrete-time generative model for smartphone touches. This model extends existing priority-based models by including refractoriness [4,7]. The output of the model is the set of touch times {t 0 , t 1 , . . . , t N }, where t i can take discrete values, i.e., t i = k i t, with t being the bin width and k i ∈ N. Equivalently, the model output can be described by the touch train s t where s t = 1 denotes the presence of a touch while s t = 0 indicates the absence of a touch.
Every touch is the result of a decision process. We assume that an individual can perform tasks from only two categories: either a task related to a smartphone screen touch or other task such as driving a car. In each category, there can be important tasks (such as dialing an emergency number) or less important tasks (such as checking the news). So we will assume that every task can be described by its priority level which is a number between 0 and 1. Let x t ∈ [0, 1] denote the priority associated with a touch task at time t and y t ∈ [0, 1] the priority associated to the other task. If at time t the touch task associated to priority x t is executed (i.e., s t = 1), then a new touch task is considered and will be attributed a new touch priority value drawn from the touch priority distribution, i.e., x new ∼ p(x). If the touch task is not executed (s t = 0), then its priority remains the same. This can be summarized as Conversely, the dynamics for the other priority y t is such that when the screen is not touched at time t (i.e., s t = 0), then it is the other action that is executed and a new priority y new must be drawn from q(y). This is summarized as To generate a smartphone touch, two conditions need to be satisfied. First, the priority x t of the smartphone action needs to be greater than the priority y t of the other action, and second, the individual must be in a nonrefractory state. Formally, the touch variable s t is sampled from the following Bernoulli distribution: where the touching intensity λ (probability per time bin t) is given by where τ = t −t is the time since last touch (t = max t k {t k < t}) and H is the Heaviside step function which guarantees that touches can only be generated when x > y and ρ is the touching rate. r(τ ) 0 is the refractory function which includes post-touch effects (i.e., right after a touch, the touch probability can be reduced). A hard refractoriness function takes the following form: where H (τ ) is the Heaviside step function [i.e., H (τ ) = 1 if τ 0 and H (τ ) else] and the hard refractory time (i.e., minimal ITI). If we relax this strong condition and allow touches for any τ > 0 (but with reduced probability when τ 0), then we can define a relative refractoriness function as a sum of basis functions, with logarithmically spaced inverse time constants, i.e., α k = α 1 β −(k−1) . We took α −1 1 = 50 ms and set β such that α −1 n = 1000 ms. Note that the set {γ k } n k=1 has to be chosen such that the condition r(τ ) 0 is satisfied for all τ 0. If r(0) 0.5, then we define the effective time constant τ * as the time for which refractoriness is half, i.e., r(τ * ) = 0.5; see Fig. 4(c). In the rare cases where multiple solutions exists for τ * (which can occur when r(τ ) is nonmonotonic), we took the maximal value of the set of solutions.
The discrete-time model described by Eqs. (1), (2), and (3) is a latent dynamical system. Note that sampling this model is slow since the complexity of this sampling scheme scales with the number of bins. Even more critical is the learning procedure for such a latent dynamical model which can be prohibitively slow for smartphone touching data sets which typically extend over months. A much faster sampling scheme is proposed below.

B. Continuous-time model
The idea of the continuous-time model is to directly sample the intervals τ instead of sampling the touch variable s t at each time step. The transition to this continuous model can be done in two steps. First, we observe that when t is small, the other priorities y t constantly change (except at the rare times where s t = 1), i.e., Eq. (2) can be approximated as y t ∼ q(y). This means that the priorities y t are independent of time and therefore, the probability of generating a touch can be marginalized over y t : where the average touching intensityλ is given bȳ and π (x) is the probability of having x > y for a given x, In the second step, we take the limit t → 0 and therefore, the intertouch interval distribution conditioned on x can be expressed as (see also Ref. [21]): The unconditioned ITI distribution is obtained by averaging the conditioned ITI distribution over the touch priority distribution p(x): So samples of the continuous-time model can be simply obtained in a two-step procedure. First, x is sampled from p(x), then τ is sampled from p(τ |x) given by Eq. (10). For this second step, one can use the time rescaling theorem [21]. Note that this continuous-time model describes a renewal process and hence the sampling complexity scales with the number of touches N.
In the absence of refractoriness [i.e., r(τ ) = 1], the sampling procedure is even simpler. First, a Poisson rateλ = ρx can be drawn from a distribution of rates p(λ) with maximal rateλ max = ρ, then τ is sampled from an exponential distribution p(τ |λ) =λ exp(−λτ ) and the ITI can be expressed as which is precisely the ITI one would get from an heterogeneous Poisson model [12]. This shows the equivalence between the priority-based model and the Poisson-based models.

A. Invariance of the model
Before giving a parametric form for all distributions, let us first note an invariant property of the model. In particular, it can be shown (see Appendix A 1) that the ITI distribution given by Eq. (11) remains unchanged if the pair of priority distributions [p(x), q(y)] is replaced by [p(x),q(y)] given bỹ (13) where φ is a differentiable and strictly monotonously increasing function with boundary conditions φ(0) = 0 and φ(1) = 1. This invariance can be understood intuitively by noting that the notion of priority contains some arbitrariness. Indeed, the only element which is relevant in the decision process is whether x is larger or smaller than y [see Eq. (4)]. If we define a new priority x = φ(x) (with the above conditions on φ), then we observe that the ordering remains unchanged, i.e., x > y ⇒ φ(x) > φ(y). This observation can also be made more formally with a change of variable in Eq. (11) (see Appendix A 1). Second, this invariance property of the model means that without loss of generality, we can set one distribution and rescale the other one. For example, without loss of generality, we can set q(y) = 1. For the touch priority distribution, we will assume that it is given by a β distribution: where With the above choice of q, the ITI distribution in Eq. (11) can be rewritten in a simpler form

B. Scale-free intertouch interval distribution
For short timescales (τ < α −1 n ), the ITI distribution is governed by the refractory function r [see Fig. 2(a)]. However, for longer timescales (τ α −1 n ), the ITI distribution follows a power-law distribution. This can be seen in two steps. First, in the limit of large τ , we have r(τ ) → 1. Second, in the limit of large τ , we know from Eq. (15) that the ITI distribution is only sensitive to the touch priority distribution in the vicinity of x = 0 that we denote as p 0 (x). Note that p 0 (x) is not normalized.
where (z) = ∞ 0 x z−1 e −x dx is the function. Therefore, the power-law exponent is given by a + 1 [see Fig. 2(b)]. So when a is large, touching tasks are more important-in the sense that there are few low-priority touch tasks-and consequently there are less large ITI.

A. Model fitting
For each subject, we fitted six versions of the priority-based model (M 1 -M 6 ) and three benchmark models (M 7 -M 9 ). Those six versions are obtained by fixing a subset of parameters to specific values in the case of no refractoriness (M 1 -M 2 ), hard refractoriness (M 3 -M 4 ) as well as relative refractoriness (M 5 -M 6 ) (see Table I). The three benchmark models (M 7 -M 9 ) are further detailed in Appendix A 4. So, in total, we fitted nine models per subject: (1) Model M 1 is the simplest model and contains only two parameters: θ = (a, ρ). It is assumed that b = 1 and that there is no refractoriness ( = 0).
(2) Model M 2 is the same as model 1 except that the touch priority distribution has 2 free parameters: a and b. Overall, it contains three parameters: θ = (a, b, ρ).
(4) Model M 4 is the same as model 3 except that the touch priority distribution is described by both a and b. It contains four parameters: θ = (a, b, ρ, ).
(6) Model M 6 is the same as model M 5 , but b is not constrained to be equal to one. The models contains therefore n + 3 parameters: θ = (a, b, ρ, γ 1 , . . . , γ n ).
For each model and for each subject, the model parameters θ are fitted from the set To do so, we relied on the continuoustime model which massively simplifies the expression of the log-likelihood. Indeed, the detailed model can be seen as a dynamical latent variable model (where the latent variables are x and y) which can be fitted through EM type algorithm but is known to be very slow. Here, because of the analytical expression of the ITI for the continuous-time model [see Eq. (15)], we can express the following objective function: which is the log-likelihood ] minus a regularization term on the coefficients γ k to prevent overfitting. This regularization term (with λ = 1000) is only used in models 5 and 6. Note that this objective function can be seen as the log-posterior with a Gaussian prior (with variance 1/2λ) on the coefficients γ k and a flat prior for the other parameters.
Because the refractory kernel must remain positive for all time, i.e. r(τ ) 0, ∀τ 0, the optimization task can be expressed as However, the difficulty of the optimization problem defined in Eq. (18)

B. Fitting results
We recorded smartphone touches from 84 individuals for an average duration of 36.5 days (see Appendix B 1 for details on data collection). The average number of smartphone screen touches per day ranged from 285 to 9915 with a median value of 2540 touches per day.
For each individual, the six different models were fitted according to the procedure described above. In particular, we first fitted the models without refractoriness (M 1 and M 2 ) and the models with hard refractoriness (M 3 and M 4 ). We found that the likelihood can be drastically improved by adding the of the data, the power-law relationship extends over 5 decades (from 10 3 to 10 8 ms).
The fitted refractory kernel [see Fig. 4(c)] shows a strong reduction of touching rate during the first few hundreds of milliseconds after the last touch. For one subject, it even displays a small increase in touching rate about 300 ms after the last touch [see Fig. 4(d)]. This smooth transition from short ITI to longer ITI removes the need to define an arbitrary onset of the power-law distribution [20].
The fitted touch priority distribution [see Fig. 4(e)] [assuming that the other priority distribution is given by q(y) = 1] diverges for small priorities (which is the case when a < 1). We repeated this fitting procedure for the 84 subjects. The population results are displayed on Figs. 4(b), 4(d) and 4(e). We found that over the population the priority parameter a is fairly scattered around a median value of a = 0.53 (for model M 6 ) and of a = 0.49 (for model M 5 ). The large interindividual differences is also highlighted in Fig. 4(g) which displays a broad distribution of touching rate ρ over the population.

C. Model comparison
To compare the different models (see Table I) for each individual, we can use the Bayesian information criterion (BIC) which is well suited for large data sets (i.e., large number of touching intervals N) which is precisely our case. BIC is given by BIC = log(N )|θ | − 2L(θ * ), where |θ | is the number of parameters and L(θ * ) is the objective function given by Eq. (17) and is evaluated at the MAP parameter θ * . given by First, we found that the three benchmark models i.e., the gamma, Weibull, and log-normal distributions are outperformed by all the priority-based models (Fig. 5). This is interesting since the exact shape of interevent distribution, although in the context of other data sets such as email correspondence, has been an object of dispute [4,22,23]; see also the review from Karsai et al. [3]. In particular, it has been argued that for email correspondence a log-normal distribution can better capture interevent distribution [22] but see also Ref. [23]. In our case, it is actually not a surprise that log-normal distribution is outperformed by a power-law distribution. Indeed, as we have seen in Fig. 4, the power-law exponent is on the order of α = a + 1 1.5, whereas the log-normal distribution has a power-law decay for large τ with an exponent of α = 1 which is significantly different from 1.5.
Second, we found that the simplest priority-based models without any refractoriness (M 1 and M 2 ) or with hard refractoriness (M 3 and M 4 ) are outperformed by models with relative refractoriness (M 5 and M 6 ); see Fig. 5. Indeed, despite their relative large number of parameters which penalizes the BIC, the models with relative refractoriness have a better (i.e., lower) BIC than the other models since they better describe short intervals. In particular, we found that the overall best model is M 6 with n = 21 basis functions. When the priority parameter b is set to one, then the best model of M 5 is when n = 20. Note that the difference between the difference in BIC between the best model M 6 and the best model M 5 is BIC pop = −3.8 × 10 4 which is highly significant.

D. Short versus long intervals
Given the fairly broad distribution of fitted power-law exponent a [ Fig. 5 We then asked whether the effective refractory time constant τ * is correlated with the power-law exponent across different subjects. Note that from the way the model is constructed, those two parameters are a priori unrelated since the refractory affects short intertouch intervals and the power-law exponent affects longer intervals (see Fig. 2). We found that a and τ * are indeed inversely correlated [ Fig. 6(a)] with an explained variance of 40% for model M 5 and 22% for model M 6 [ Fig. 6(b)]. This indicates that subjects that have a fast motor control (i.e., have a small τ * ) also put a higher priority on their smartphone (i.e., higher a) in the sense that they have less low priority smartphone tasks.

V. DISCUSSION
We proposed a generalized priority-based model which is both flexible and tractable. The flexibility comes from the set of basis functions which describe refractory effects at short intertouch intervals, while the tractability stems from the simplified structure of the generative model in continuoustime which enables a fast fitting procedure. The flexibility is essential to capture interindividual differences in touching behavior while the tractability is crucial for fitting large data sets. We found that the intertouch intervals are better captured by the priority-based model than by other reference distributions such as log-normal, gamma, or Weibull distributions. We also found that the interindividual differences in low level motor control ability (reflected by the effective refractory time constant) can be partially explained by the higher-level cognitive processes which attributes priority to specific tasks (reflected by the priority parameter a). (c) The distribution of the difference a n (s) − a n−1 (s) across subjects s = 1, . . . , 84 and across basis function numbers n = 10, . . . , 30 (black, σ fit = 0.022) is much narrower than the distribution of the differences a n (s) − a n−1 (s ) where the indices s are randomly permuted (white, σ rand = 0.28). (d) Same as panel (c) but for model M 6 (σ fit = 0.043, σ rand = 0.23). (e) Similarly as in panel (d), the distribution of the difference of effective refractory time constants τ * n (s) − τ * n−1 (s) (black, σ fit = 88.7 ms) much narrower than the one obtained when subjects are permuted for each n (white, σ rand = 458 ms). (f) Same as in panel (e) but for model M 6 (σ fit = 100 ms, σ rand = 366 ms).
Other models describing smartphone activity have been proposed. However, they are aimed at addressing different questions (mostly sleep related) and use different type of data. The model of Cuttone et al. [24] aims at predicting sleep patterns and rely on app launch timings binned over 15-min duration and therefore lack the possibility to describe refractory effects on the tens of milliseconds resolution. Abdullah et al. [25] only considered screen-on and screen-off events to predict sleep patterns.
More generally, circadian rhythms have received recently a lot of attention [26][27][28] and the question has been addressed whether circadian rhythms could explain heavy tail distribution [14] has argued that a cascade of Poisson processes can give rise to power-law distribution. Interestingly, Jo et al. [19] showed that even if the data is deseasoned (i.e., the circadian and weekly patterns are removed from the time series of mobile phone events), the heavy-tails remain.
Closer to the present study, the priority model proposed by Refs. [2,4,7] already predicts a power-law distribution. However, our model deviates significantly from their approach in several aspects.
First, and most importantly, the priority-based model from Barabási is aimed at predicting the response time distribution (e.g., the time interval between the reception of a letter and its response), whereas in our case, our priority-based model predicts interevent times (e.g., the time interval between two consecutive posted letters-or in our case intertouch intervals). It is, however, interesting to note that even though the response-time and the interevent-time metrics are clearly different, they can still be seen within the same umbrella. On the formal level, the response-time priority-based model with a list of L = 2 is equivalent to the discrete-time emission model described in Sec. II A. Furthermore, on the interpretation level, since a new smartphone touching task arrives just after every smartphone touching event [see Eq. (1)], one could still reinterpret the following interevent time as a response time. For a further discussion on the relationship between responsetime distribution and interevent distribution, see Ref. [3]. Second, on the conceptual level, those authors stress the universality of the various behaviors. Their claim is that the interevent intervals for certain activities such as browsing the Web, sending emails or loaning books fall into a specific universality class with power-law exponent of α = 1 while the response time for other activities such as writing mails follow another universality class with exponent of α = 3/2. Actually, more recently, Formentin et al. [29] elegantly showed that various response-time datasets (on sms [30], email [29], and mail [1]), with apparently different power-law exponents can be cast within the same α = 3/2 universality class provided that the events are properly reclocked [31]. Here, we found that the power-law exponent (averaged over the population) is α = a + 1 1.56 ± 0.16 which is indeed close to the rational exponent of 3/2. However, it should be noted that the powerlaw exponent of individuals are fairly spread ranging from α = 1.31 ± 0.04 to α = 2.19 ± 0.04 which are clearly different from α = 3/2. Capturing those nonuniversal exponents is possible in our model since the power-law exponent is given by a + 1 where a can take any real positive value. In contrast, in the work of Oliveira et al. [7], the exponent is determined by the length of the list of tasks [32].
The third difference w.r.t. the studies of Refs. [2,4,7] is that our model has been actually fitted to the whole set of event times (the touch times) such that we did not neglect small interevent intervals by defining a (somewhat artificial) onset of the power-law distribution [20]. This is possible in our model since short intervals are captured by the refractory kernel. Note that even though refractory kernels have been used in other fields (e.g., in spiking neuron models, the probability of generating a spike just after a first one is also modulated by a refractory kernel [33][34][35]), the particularity here is that the specific form of the refractory kernel is such that its integral can be computed analytically which boosts the computational efficiency.
Finally, fitting our model to the touching data has been possible because we considered the continuous-time priority model. Indeed, the marginal likelihood can be expressed analytically for the continuous-time model (and not for the discrete time model) which makes the maximum-likelihood parameter learning extremely efficient.
A separate line of research based on biological signals has also encountered scale-invariant relationships referred as 1/ f pink noise [36,37]. However, those studies compute the power-spectrum density and not the interevent distribution. Actually, if we do compute the power-spectrum density for the smartphone touching model, then we find that in the limit of large frequencies, the power-spectrum density remains constant and does not decrease as 1/ f .
Here, this generalized priority-based model has been applied to smartphone touching data but could be applied to other event-based data sets which display power-law property for large interevent intervals such as surface mails, emails, or even foraging patterns.
The simulations have been performed using Matlab. The code and the smartphone touching data are available online [38].

ACKNOWLEDGMENTS
We thank Simone Carlo Surace for very helpful discussions. J.P.P. was supported by the Swiss National Science Foundation Grants No. PP00P3_150637 and No. PP00P3_179060. A.G. was supported by the Society in Science, the Branco Weiss Fellowship, and a research grant from the Holcim Foundation.
J.P.P. and A.G. designed the study, developed the model used here, and drafted the manuscript. J.P.P. analyzed the data and formulated the mathematical model, aided by A.G., who collected the data and conceived the project with J.P.P.
A.G. is an inventor of the patent-pending technology used to track touchscreen interactions in this study. A.G. and J.P.P. are cofounders of QuantActions GmbH, a company focused on quantifying human behavior through smartphone interactions.

Invariance of the model
In this section, we will show that that the ITI distribution remains unchanged if the pair of priority distribution [p(x), q(y)] is replaced by [p(x),q(y)], wherep(x) andq(y) are given by Eq. (13).
Let us consider the following change of variable: x = φ(x ). The ITI distribution can be therefore expressed as where the conditional ITI distribution p q [τ |φ(x )] depends on the other priority distribution q(y) via the instantaneous ratē λ q [φ(x ), τ ] which can be expressed as whereq is given by Eq. (13). Note that the dependence on q is included only here for the clarity of the argument, but is omitted otherwise for the simplicity of the notation. Therefore, the ITI distribution is invariant under the change of variable φ for both x and y. Indeed, we have For example, if the touch priority distribution is given by p(x) = β(x; a, 1) and the other priority distribution is given by q(y) = β(y; a , 1), then the function φ(x) = x k allows us to generate a family of equivalent pairs of priority distributions [p(x),q(y)] = [β(x; ka, 1), β(y; ka , 1)]. Therefore, the ITI remains unchanged as long as a/a remains constant. Note that this argument can be generalized to arbitrary smooth distributionp andq. LetÃ(x) denote the ratio of the logarithm of both cumulative density functionsP( Sincep andq are smooth, when x → 0, we can express those priority distribution as p(x) c 1 x a−1 and q(y) = c 2 y a −1 . Also, since φ(x) is smooth, it can be approximated as φ(x) = x k in the vicinity of x = 0. Now we can show that theÃ(0) is independent of the change of variable function φ. Indeed, is independent of k.

Log-likelihood gradient
For the models with relative refractoriness, we fitted the parameters θ = (a, b, c, γ 1 , . . . , γ n ) by performing maximum likelihood with a suitable regularization for the parameters γ i . Note that for a practical implementation, it is easier to learn c = log(ρ) instead of ρ itself. For a set of intertouch intervals , the log-likelihood can be expressed as where the expectation · is w.r.t. p(x) = β(x; a, b) and the function E i (x) is given by and R(τ i ) is given by By noting that we can compute the log-likelihood gradient w.r.t. a: By symmetry, the gradient of L w.r.t. to b yields The gradient of L w.r.t. c is given by Finally, the gradient of L w.r.t. γ k can be expressed as For the models with hard refractoriness, the integral over the refractory kernel is given by when 0 < τ min where τ min = min i τ i is the minimal intertap interval. Under this condition, the log-likelihood gradient yields which is positive as long as τ min . When > τ min , the log-likelihood goes to −∞. Therefore, the optimal refractory time constant is * = τ min .

Computing the integrals
Both the log-likelihood L as well its gradient w.r.t. to the parameters θ contain integrals that are delicate to evaluate. Indeed, the integrand of all those integrals depend on the beta distribution β(x; a, b) which can diverge at x = 0 or x = 1 depending on the parameters a and b. So whenever possible, we compute those integrals analytically. This can be done for the following integrals: where ψ (z) = d log (z)/dz is the digamma function and B(a, b) = (a) (b)/ (a + b) is the beta function. By symmetry, we have By Taylor expanding the exponential in the expression of E i (x), the integral xE i (x) a,b can be expressed as where 1 F 1 is the hypergeometric function defined as and (a) k = k−1 i=0 (a + k) for k 1 [and (a) 0 = 1] is the rising factorial (also called Pochhammer function). Similarly, can be expressed as (A20) When it is not possible to compute the integrals analytically, the idea is to express the integral as a sum of two integrals where the first one is well suited for a numerical integration and the second one can be performed analytically. For example, xE i (x) log(1 − x) a,b can be computed as where the first term of the r.h.s can be computed numerically and the second term can be computed with Eq. (A17). Finally, it should be noted that the integral xE i (x) log(x) a,b can be computed numerically straightforwardly since the integrand does not diverges when x = 0 nor when x = 1.

Benchmark distributions
We benchmarked the generalized priority-based model is benchmarked with three simple models: the gamma distribution, the Weibull distribution, and the log-normal distribution. The gamma distribution is given by where κ is the shape parameter and ϑ is the scale parameter. The log-likelihood is therefore expressed as a function of a two-dimensional parameter θ = (κ, ϑ ) and is given by where the the expectation symbol here denotes the average over the data points, i.e., f (τ ) τ ≡ N −1 N i=1 f (τ i ). The maximum-likelihood parameters θ * can be found by maximising Eq. (A23). This can be done easily by taking advantage of the fact that the optimal scale parameter has to obey ϑ * = τ τ /κ. So the shape parameter can be found by evaluating the gradient of the log-likelihood w.r.t. to κ at ϑ * and setting the gradient to zero. This amounts to finding the root of where ψ (k) is the digamma function. The Weibull distribution is given by where κ is the shape parameter and λ is the scale parameter. The log-likelihood is therefore given by with θ = (κ, λ). By observing that the maximising scale parameter is given by λ * = k τ κ τ , we can obtain the shape parameter k by computing the root of Finally, the log-normal distribution is given by where the maximum likelihood parameters are simply given by (σ * ) 2 = (log(τ ) − μ * ) 2 τ . (A30) APPENDIX B

Smartphone data collection
A custom-designed software application called Touchometer (now available as TapCounter, QuantActions, Lausanne) that could record the touchscreen events with a maximum error of 5 ms [39] was installed on each participant's phone.
To determine this accuracy, controlled test touches were done at precisely 150, 300, and 600 ms while the Touchometer recorded at 147, 301, and 600 ms, respectively, with standard deviations less than 15 ms (interquartile range less than 5 ms). The app posed as a service to gather the timestamps of touchscreen events that were generated when the screen was in an unlocked state. The operation was verified in a subset of phones by using visually monitored tactile events. The data were stored locally and transmitted by the user at the end of the study via secure email. One subject was eliminated as the app intermittently crashed after a software update. The smartphone data were processed by using MATLAB (MathWorks, USA).