Amortized Inference with User Simulations

The inverse modeling process ascertains the best model parameters for a given dataset. After inversion, one can approach the model’s parameters as hypotheses about the latent characteristics of interest. Where the model represents the user’s behavior via a simple closed-form expression, researchers can, from the parameters, calculate the likelihood of any given observation or a particular discrepancy between a prediction and what actually happened. Fitts’ law [6, 7, 17, 42, 63], the diffusion model for binary decisions [64, 65], the classic model for visual perception of speed [71], and the integrated eye-movement model [67] are some of the many traditional user models with parameters determined through MLE or least-squares estimation. This approach does not permit the simulation models’ inversion, however: since most user simulation models involve complex relationships without a known likelihood function, likelihood-based fitting methods (MLE etc.) cannot invert them. Hence, many previous studies adopt values from the literature or tune parameters manually.

ABC has emerged as a popular likelihood-free approach for inverting simulation models [4, 11, 27]. It offers a systematic way of finding the parameters that replicate a given dataset well in the model’s parameter space and thus revealing the parameters’ posterior distribution from the observed data. Kangasrääsiö et al. [34] first attempted to apply ABC to fit an RL-based simulation of menu selection. However, its application to user simulation models is plagued by the time-efficiency problem. Running the multiple simulations required with various parameter candidates is time-consuming when the simulation model insists on optimizing its control policy anew for each candidate. Moon et al. [49] introduced a way of addressing this. They showed that the control policy of user simulation models can be generalized across varying model parameters by means of multi-task RL. Consequently, the ABC process, which had previously required days of computation [34], may be reduced to one or two hours. Nonetheless, as is hardly surprising, individual-level fitting based on ABC is normally restricted to smaller numbers of users (e.g., 5 [34] or 20 [49]). Though one recent study [18] of task-interleaving behavior employed ABC to fit an RL-based model to the data of 211 individuals, fitting this model to an individual user took 1.5 hours of computation time. Our motivation in delving into amortized inference was to see whether it could overcome these issues, thereby opening the door to larger-scale and real-time applications of inverse modeling in HCI.

Amortized inference consists of “investing upfront computation to support rapid online inference” [72]. With modern ML methods, a complex probability distribution can be estimated through variational inference [5, 38, 77], which uses optimization of tractable and parameterized distributions to approximate the target distribution. One can amortize this variational inference by building a neural proxy model, a conditional density estimator that learns an effective mapping from a given observation (e.g., an individual’s behavior) to an approximate distribution (e.g., a posterior of model parameters) [38]. Once trained through upfront computation, the conditional density estimator enables rapid inference of the posterior of the parameters conditioned on the given observation.

Researchers have examined many ML methods in efforts to solve this conditional density estimation problem  [11, 19, 45, 55]. One critical factor is the neural networks’ representation capacity – they must be able to learn the complex probability distribution. Early methods combined Gaussians, to estimate the posterior distribution, with trained neural-network models (e.g., mixture density networks), to predict the means and covariance of the Gaussians [57, 62]. Though these methods performed sufficiently well for several inverse problems, it shows limited ability to represent more complex posterior distributions (multimodal ones especially). As a emerging approach, normalizing flows [58, 66] approximate complex distributions through sequential bijective transformation steps, starting with a simple normal distribution. Here, invertible neural networks (INNs) model the bi-directional conversion between the target distribution and a tractable one. This invertibility constitutes the unique advantage of normalizing flows over other variational inference models (e.g., a variational autoencoder [38]). A tractable distribution that can be converted to the target distribution in a bijective manner renders it possible to calculate that target’s posterior density precisely and provides for easy sampling of the parameters from it. Both properties (tractable density and ease of sampling) are beneficial in the training and inference processes for density estimators. Therefore, recent research into amortized inference has welcomed normalizing flows [2, 21, 23, 40, 58]. BayesFlow [61] is a good example of recent work demonstrating the utility of normalizing flow models (RealNVP [14]) for amortized inference. Non-HCI applications of amortized inference have exploited its flexibility; it is easy to apply to most simulation models yet does not necessitate critical assumptions [13, 22, 60], but the HCI field has been devoid of attempts to amortize inference for simulation models, with the exception of Murray-Smith et al.’s successful inference of a finger’s 3D pose from sensor data via a conditional variational autoencoder [51]. We fill this gap with the first presentation of a general workflow for applying amortized inference to the HCI field’s simulation models.

We express the operation of a user simulation model as y = f(x;θ), where y represents the model’s output (i.e., user behavior), x represents the input (i.e., task state), and θ represents the parameters (i.e., a user’s cognitive-physiological characteristics). The prior of model parameters is denoted by p(θ), and their posterior when behavior data are given is p(θ|y). We use the hat symbol (∧) to highlight variables for data estimated or predicted through models (e.g., \(\boldsymbol {\hat{\theta }}\) for the estimated model parameters and \(\boldsymbol {\hat{y}}\) for the predicted user behavior). Finally, the subscript “o” refers to data observed from real users, as in y_o for the observed user-behavior data.

With user simulation models, we can generate a training dataset for the density estimator, consisting of plausible sets of θ and simulated y under them, in the form of (θ, y) pairs. Recent simulation models reproduce user behavior (y) through the sequential decision-making of an RL agent: the agent obtains partial (or noisy) information about the task environment’s state in keeping with an individual’s cognitive-physiological bounds (dictated by mechanisms’ limits [47, 67, 71] and variable parameters, θ) and chooses an action that achieves maximum utility. In RL terms, the simulated user’s internal decision-making process is formulated as a partially observable Markov decision process (POMDP). The agent’s decision-making function (i.e., control policy) can be implemented as a neural-network model that, from the observed state, ascertains the action, and deep RL techniques such as the deep Q-network (DQN) algorithm [15, 30, 48] afford optimizing it. One problem that most current simulation models face here is that their control policy is usually optimized within fixed bounds (i.e., with fixed θ) [10, 15, 30]. To collect the training data under various θ values, one must adjust the control policy for each θ. Rather than take the time-consuming path of optimizing a control policy anew for each θ, this workflow follows recent multi-task RL approaches; i.e., a single generalized control policy is trained for a family of tasks (solving for a POMDP family parameterized by θ). Here, a neural-network-based policy model is implemented that receives θ as external input alongside the observed task state. Modulated Q-network [49] and optimal control ensembles [41] exemplify such multi-task RL approaches’ recent use. By optimizing the policy model subject to the multi-task environment (i.e., episodes with various θ), we are able to achieve a generalized optimal policy. In consequence, the training dataset can be simulated, where θ ∼ p(θ), y = f(x;θ).

The conditional density estimator’s objective is rapid estimation of the posterior (p(θ|y)) from given behavior data (y). For the general case, it has to meet two crucial requirements. Firstly, it must deal with various data types present in user-behavior observations (e.g., a simple summary of multiple task trials or full details of each) and extract meaningful features from them as a fixed-size vector. Secondly, when given the extracted features, it should be capable of representing even highly complex distributions of p(θ|y), including multiple modes (peaks). To address both requirements, we designed the density estimator to consist of an encoder network and a conditional INN (see Figure 2). The estimator learns bi-directional conversion between a latent variable (denoted as z) from a normal distribution \(\mathcal {N}(\boldsymbol {0}, \boldsymbol {{I}})\) and θ from a complex distribution p(θ|y), when given y. Therefore, θ = g_ϕ(z;y) and \(\boldsymbol {z}=g^{-1}_{\boldsymbol {\phi }}(\boldsymbol {\theta };\boldsymbol {y})\) , where g_ϕ represents the forward conversion through the density estimator with neural-network weights ϕ. Because the latent variable (z) follows a normal distribution, it is easy to sample z and calculate its exact density. Therefore, with the bijective transformation between z and θ, we can achieve easy sampling and density estimation for p(θ|y) as well.

With the trained density estimator, we can rapidly infer the posterior of model parameters (θ) from real users’ observed behavior (y_o), as illustrated in Figure 3(b). Sampling the latent variable z from \(\mathcal {N}(\boldsymbol {0}, \boldsymbol {{I}})\) and computing θ via g_ϕ(z;y_o) is equivalent to sampling θ from its true posterior p(θ|y_o) if one assumes perfect convergence of density-estimator training [61]. Therefore, obtaining a plausible set of θ from the given y_o comes at very low computational cost, only that of passing the sampled z and y_o through the trained density-estimator network, which can be performed on the order of milliseconds. With the set of plausible θ candidates, we can approximate the full posterior distribution p(θ|y_o) by means of well-known methods, such as kernel density estimation (KDE), and determine the estimated values of model parameters \(\boldsymbol {\hat{\theta }}\) , typically using maximum a posteriori (MAP) estimation [34].

Should distribution information be unnecessary, one can simplify the estimator model’s structure by seeking only point-estimate values of the parameters. This can yield the benefit of decreased size of the neural-network model, thereby reducing the computation costs for training and inference. The workflow presented above supports training a point estimator accordingly. Instead of a conditional INN in a density estimator, one can build a point estimator with any feedforward network (e.g., an MLP) that directly outputs a point estimate from a feature vector extracted from the encoder network. With the same training dataset, the process would use a general regression loss (e.g., mean squared error) rather than the training loss of Equation 1. The inference process would consist only of feeding forward the given observation data, without any need for sampling of the latent variable z. Further on (in Section 8.3), we report the quantitative results from testing how using point estimators as opposed to density estimators affects inference performance.

With the workflow described in Section 3, we applied amortized inference as presented below.

We evaluated the performance of amortized inference from the following two perspectives: how well the conditional density estimator can recover ground-truth model parameters from the simulated data (parameter recovery) and how well the parameter values inferred from the empirical user data can actually predict the user’s behavior (behavior prediction).

To judge the former, we randomly sampled 100 sets of θ from p(θ) and simulated y from them. Then, we measured the coefficient of determination (R²) between the ground-truth θ and the density estimator’s prediction ( \(\boldsymbol {\hat{\theta }}\) ). Here, R² represents the proportion of the variation in the predicted \(\boldsymbol {\hat{\theta }}\) that can be explained by variation in the true θ. Hence, higher values are better; the range is 0 to 1, so R² = 1 denotes a perfect prediction for θ.

Next, we assessed the model’s behavior-prediction performance with real users’ data when using the fitted parameters with amortized inference. We investigated the two scenarios for applying amortized inference: group-level and individual-level fitting. For group level, one model parameter set was estimated from the all-users dataset. In the individual-level fitting, on the other hand, model parameter sets were estimated for each individual user’s behavior data. We measured the prediction accuracy as the distance between the simulated behavior under the estimated parameters ( \(\boldsymbol {\hat{y}^\prime }\) ) and the actual behavior of the same users ( \(\boldsymbol {y_o^{\prime }}\) ).³ We considered three distinct distance metrics: 1) mean difference, 2) KLD,⁴ and 3) maximum mean discrepancy (MMD) [24]. Mean difference captures the absolute difference between the mean values from the behavior features in the trials in each case (ground-truth observation vs. model prediction). In contrast, KLD and MMD estimate the difference between distributions. The closer KLD and MMD are to 0, the closer the behavior-feature distributions of \(\boldsymbol {\hat{y}^\prime }\) and \(\boldsymbol {y_o^{\prime }}\) are.

We measured the distance with regard to four features of menu-search behavior: each trial’s completion time and number of fixations, when the target is present and when it is absent (2 × 2). We compared the resulting fit (group- and individual-level) from amortized inference with the baseline: use of the earlier model parameters fitted via ABC [34]. This baseline represents applying conventional non-amortized inference to fit the parameters at group level.

The inference process, run a single time for each dataset, took about 10 ms with our trained density estimator (i.e., both group- and individual-level fitting). Accordingly, individual-level inference for 18 users consumed only about 200 ms in total. This is a significant improvement over previous (ABC) methods, which took 37 h for group-level fitting for the same user dataset [34].

Figure 5(a) depicts parameter-recovery performance with our trained density estimator, showing the correlation between the ground-truth values and our predicted values for each model parameter. The results show that the trained estimator predicts most parameters with a high coefficient of determination (R²=0.74 for d_fix, 0.87 for d_sel, and 0.82 for p_rec). It was more difficult to estimate p_sem from the given behaviors (R²=0.33).

Figure 5(b) shows the estimated posterior distribution from the actual behavior measured from all users. Our estimated posterior distribution closely matches the baseline model parameters, which lie within its high-probability region. Table 1 shows the distance between the observed behavior and the three model predictions. While the baseline, group-level, and individual-level model fits with amortized inference all replicate the observations reasonably well (see Figure 6 for the full histograms of observed and simulated behaviors over trials), from among the three it was our individual-level fitting that made the most accurate prediction, with the smallest distance by nine out of the 12 metrics (4 behavior features × 3 distance metrics). On the other hand, there was no noticeable difference between the baseline and our group-level fit, which showed the smallest distance by one and two metrics, respectively; that is, the amortizing approach offered inference performance comparable to ABC’s (in the same group-fitting cases) while vastly improving the computation cost (from 37 h to 10 ms).

Most aspects of the implementation for amortized inference were the same as in Case 1. We implemented the control policy of the point-and-click model with a modulated Q-network, using the prior work’s network structure and optimization method [49]. This generalization of the control policy enables single-model simulation of point-and-click behavior y under the various values of θ. Supplement C.1 details the control-policy structure and optimization.

The simulation’s prior distributions p(θ) too were determined from the earlier work [49]: the prior for σ_v is a log-uniform distribution where (min, max) = (0.069, 0.415); that for n_v is log-uniform with (min, max) = (0.145, 0.413); and that for c_σ is a log-uniform one where (min, max) = (0.055, 0.400). The previous report did not cite any prior knowledge of the range for T_{h, max}, so we used a simple uniform prior for it where (min, max) = (0.5 s, 2.5 s). We sampled 262K distinct θ values from p(θ) and simulated 32 trials for each sampled θ. In total, 8M simulated point-and-click trials were generated for density-estimator training. With 16 CPU cores, the simulation process took seven hours.

In contrast against Case 1’s y, consisting of only summary features across trials, Case 2’s consists of multiple trials’ data. The number of observed trials can be varied across users (i.e., y size can be varied). Accordingly, we implemented the attention-based encoder network to extract fixed-size features ( \(\boldsymbol {\tilde{y}}\) ) from variable-size observations (y). One million gradient descent steps (with batch size 32) were performed with the simulated training dataset. The total training time, approximately 30 h, was longer than Case 1’s because of the more complex observation data and network structure involved. Supplement C.2 elaborates on the density-estimator model implementation and training.

We evaluated the performance of our amortized inference from the same angles considered in Case 1: parameter recovery and behavior prediction (detailed in Section 5.2). With regard to the first, we assessed how well the trained density estimator could infer 100 sets of ground-truth parameter values from the simulated behaviors. For Case 2, the investigation additionally considered the estimator’s ability to recover real users’ measured cognitive-physiological capabilities from each user’s point-and-click behavior.

In this case too, we assessed behavior-prediction performance with regard to baseline, group-level, and individual-level fit. We measured the distance between the point-and-click behavior simulated via the fitted model parameters ( \(\boldsymbol {\hat{y}^\prime }\) ) and actual user behavior ( \(\boldsymbol {y_o^{\prime }}\) ) in terms of three features: 1) trial-completion time; 2) the cursor’s final click location relative to the target, with the target’s radius as the unit (hence, a value below 1 indicates that the click fell within the target); and 3) the cursor’s total travel distance in each trial. The baseline case for comparison came from simulated behavior where the participant-specific cognitive-physiological capability values [49] furnished the model’s parameters; therefore, the baseline here represents obtaining each individual user’s unique characteristics to inform the simulation, through laborious measurement processes.

Since each dataset’s inference process took 125 ms with our trained density estimator, individual-level inference for 20 participants’ datasets used just 2.5 s. In contrast, the previous study reported 30 h for the same individual-level inference (using ABC alongside a policy-modulation technique for rapid model simulation).

Figure 7 presents the parameter-recovery performance with the density estimator, for both simulated behavior data and actual users’ behavior data. With the former, the density estimator predicts the model parameters with high (R²=0.94 for σ_v, 0.67 for c_σ, and 0.96 for T_{h, max}) and moderate (R²=0.38 for n_v) coefficients of determination. From the actual user behavior, it was possible to predict their cognitive-physiological capability values moderately well (R² = 0.37 for n_v and 0.61 for c_σ) and at low levels (R²=0.23 for σ_v). This inference performance was comparable with that of the ABC-based prior technique [49], which estimated the same participants’ capability values σ_v and c_σ moderately well (with R² values of 0.50 and 0.62, respectively) but failed to infer n_v. The R² values for the inferred σ_v are slightly lower in our amortized inference case, but we succeeded in inferring n_v and achieved a huge speed improvement via amortized inference (to 2.5 s from 30 h).

Figure 8 plots the estimated posterior distributions of each model parameter from using our density estimator on the aggregated behavior data of all participants. The averaged values of the participants’ measured capabilities (σ_v, n_v, and c_σ) lie well within the high-probability region of the estimated posterior distributions. Table 2 shows the distance between the actual observation and the three model predictions: baseline, group-level, and individual-level fits (Figure 9 presents the full histogram for each simulated behavior). From among the three models, the individual-level fit (with amortized inference) showed the results closest to actual user behavior by the majority of metrics (5 out of the 9). That is, the individual-level model fit outperformed the baseline method, in which the parameters were determined in line with each individual user’s measured capabilities.

The control policy of the original touchscreen typing model, composed of two distinct network models (actor and critic networks), was obtained by proximal policy optimization (PPO) [70]. We generalized this control policy to operate under a range of model parameters, in the same manner as for the modulated Q-network used in Case 1 and 2. We implemented and trained both networks to take given model parameters (θ) in addition to the observed task state as inputs. Such a method of conditioning for actor-critic networks on the basis of given parameters has been verified [41] under the name “optimal control ensembles.” Supplement D.1 gives further details. The prior distributions of the model parameters were set such that the one for p_obs is a uniform distribution where (min, max) = (0, 1); the prior for α is a truncated normal where (mean, std, min, max) = (0.6, 0.3, 0.4, 0.9); and that for k is a truncated normal where (mean, std, min, max) = (0.12, 0.08, 0.04, 0.20). We sampled 66K distinct θ values from p(θ) and simulated 16 trials for each one sampled. Therefore, 1M simulated typing trials were obtained for density-estimator training. The simulation process took 14 hours for 16 CPU cores.

Since y consists of multiple trials’ data (as in Case 2), we implemented our attention-based encoder network to extract fixed-size features from the variable-size observations; for the training, we used 600K gradient descent steps (batch size = 100), and training took 15 h in all. Supplement D.2 provides details.

As with the other cases, we assessed parameter recovery (the accuracy of inferring 100 sets of ground-truth model parameters from the given simulated behavior) and behavior prediction (the accuracy of replicating real users’ behavior through simulation of group- and individual-level fitted models). For evaluating the latter, we measured the distance between the real user’s behavior and the model’s prediction, considering the following four behavior features: 1) WPM, 2) error rate, 3) backspace presses, and 4) KSPC. We used behavior simulated on the basis of the previously reported model parameters [30] as a baseline for comparison. In that the original study borrowed the parameter values from prior literature without fitting the model to empirical datasets, the baseline in Case 3 represents common practice in the HCI field. With Case 3, we showcase one of the promising areas for applying high-computational-efficiency amortized inference: large-scale analysis of user data. With reference to the large-scale dataset [56], we can report on how the inferred cognitive parameters of individuals are distributed in a population (N=1,057), in terms of two demographic factors: age and gender.

The inference process took about 15 ms per given behavior dataset. Accordingly, the time for the individual-level inference of 1,057 users’ data with amortized inference was only 20 seconds in total.

Figure 10(a) depicts the parameter-recovery performance from the simulated behavior with our trained density estimator. Our trained estimator predicted all three model parameters with a high level of coefficient determination (R²=0.82 for p_obs, 0.71 for α, and 0.96 for k). Figure 10(b) plots the estimated posterior distributions from the aggregated behavior data of all 1,057 users. The baseline model parameters (adopted from the literature) fall within the estimated posterior distribution but not its high-probability region. We interpret this as related to a limitation of using values imported from the literature: they are not always entirely representative of the dataset, and other confounding factors may be present. Our inferred parameters’ superiority to the literature-derived baseline ones for behavior prediction may support this interpretation. Table 3 attests that, by most metrics, our group-level fitted model parameter values replicate real user behavior better than the baseline parameter values do. Furthermore, fitting the model at individuals’ level yielded even better prediction performance, with the predictions best approaching real user behavior. The full set of behavior histograms is reproduced in Figure 11.

Figure 12 shows how the parameter values inferred for each individual participant are distributed population-wise in terms of both age and gender. We detected a significant correlation between the users’ age and inferred k, which represents the motor resources behind a user’s finger (Pearson’s r=0.18, with p<0.001). Therefore, data from typing behavior might reveal that younger people show more extensive finger-related motor resources (i.e., lower k). This result is consistent with previous work [69]. We found no significant age correlation for the other two parameters, p_obs (probability of noticing an error from finger movement) and α (the movement’s speed–accuracy bias), with p-values of 0.580 and 0.088, respectively. On the other hand, regarding gender, we discovered a significant difference between male and female user groups’ inferred k values (independent-samples t-test: t=3.20, with p<0.001). This result revealed finger motor resources to be stronger in female than in male users, even though the latter were older (the female users’ mean age being 29.75 years and males’ 25.10). The other two parameters did not differ significantly in either respect between the gender groups (p=0.151 for p_obs; 0.457 for α).

Testing with out-of-distribution data – data subject to distributions different from the training data’s – can clarify the robustness to unseen user behavior. To test the density estimator with out-of-distribution data, we assumed a situation wherein differences in the task variables given to a user in each trial cause different user behavior to emerge. We split the dataset used with the point-and-click model (Case 2) in half in the task-variable space consisting of the given pointing target’s radius and speed. Three distinct experiment setups were prepared, full data, in which the estimator was trained on the entire dataset and tested on each split dataset; same split, with the estimator trained on and tested on the same subset; and crossed split, with testing on one portion of the two-part split and testing on the other. Table 4 shows the results for averaged parameter-recovery performance (i.e., R² values) for both split datasets in each condition. The R² values for both parameters (n_v and T_{h, max}) were lower in the crossed split condition than in full data or same split. Despite the performance decreases, however, the density estimator in crossed split was able to recover all parameters, with high (σ_v, c_σ, and T_{h, max}) or low (n_v) R² levels, from the data of the unseen task.

Having strong priors of model parameters can enhance inference performance. In inverse modeling problems in HCI, the priors are the scientifically plausible distributions of a user’s cognitive-physiological characteristics. Therefore, they are often set to have a peak at the most plausible external values validated by the literature [34]. To measure the impact of priors, we compared the following two conditions, for ways in which a practitioner may set the priors: Literature-prior, with distributions weighted at plausible values adopted from the literature (as in the case studies), and Uniform-prior, using flattened distributions from the literature prior, without any weighting for external values. We evaluated the density estimator trained with each prior for all three cases. The results for the touchscreen-typing task (Case 3), which is the most challenging scenario, are presented below (Supplement E presents the results for cases 1 and 2). Each estimator’s parameter-recovery performance degrades when it grapples with synthetic data from the other priors (see Table 5). In fitting of the model to real user data at the level of individuals, using the literature-based prior outperformed use of the uniform prior by most behavior-prediction metrics (see Table E5 in Supplement E). However, we discovered that even the naive Uniform-prior density estimator could outperform previous work’s fitting baselines (ABC- or hand-tuning-based).

As researcher needs dictate, one can train a point estimator instead of a density estimator, using the same workflow (see Section 3.5 for details). We compared the two estimator types’ inference performance, with all three case studies. Tables 6 shows the results for Case 3 (Supplement F presents the results for the other two cases). We implemented a point estimator by replacing the conditional INN with an MLP that has two hidden layers (of 512 units each), which directly output the optimal parameter values. There was no noticeable difference in parameter recovery between the two estimator types when tested on a simulated dataset. In addition, when fitted to the actual user dataset at individuals’ level, they were comparable in performance, both proving superior to the baseline (See Table F5 in Supplement F). One notable difference was in the computation cost per inference: the point estimator, taking about 5 ms, was three times faster than the density estimator.

The complexity of the behavioral data behind inference (y) affects computational efficiency. In Case 1’s conditions, y consists of summary features, 1-D data with only a feature dimension. The other cases had a trial dimension added (since all data, from multiple trials, were used). Case 2 involved a time dimension additionally (because each trial’s data included a time-series trajectory), so it had the most complex y (3-D data, with trial, time, and feature dimensions). This necessitated a more complicated encoder network, which, in turn, increased the inference time. Nonetheless, inference time did not grow with the number of observed trials in y (Figure 13(a)), thanks to the parallel computation of neural networks (i.e., multiple trials’ computations can be processed as a batch).

We discovered that parameter-recovery performance may differ not only between models but also within the parameters of a model. While our trained density estimator could achieve decent recovery performance for most parameters, a few parameters showed lower R² values (e.g., p_sem in Case 1 and n_v in Case 2). One influential factor is the identifiability of the model parameters. If there are multiple solutions for a certain parameter that could produce the same model output, that parameter can be difficult to isolate from the given data. This is sometimes due to the model itself (e.g., model sloppiness [26]) or small quantities of given data. One can put amortized inference to use for rapidly revealing gaps in the given data as one examines the convergence of parameter-recovery performance with increasing dataset sizes. For example, from Figure 13(b) we can identify the number of observed trials at which performance starts to converge; we need at least 16 trials to recover k with nearly optimal accuracy and 128 trials for α in Case 3. The approximation gap and amortization gap are two additional elements that could influence parameter-recovery performance [12]. The approximation gap refers to error resulting from limits to the density estimator’s representation power; that is, it results from our attempt to approximate the true posterior distribution by using a surrogate (parameterized) distribution. The amortization gap refers to error caused by training a density estimator for the entire dataset, as opposed to optimizing an approximate distribution for each individual user’s data.

Amortized inference fitted the parameters better at the level of individuals than at that of the whole user group. Across all three cases, parameters fitted at individuals’ level showed the greatest accuracy. That said, there were exceptions (e.g., individual-level fitting was not superior to group-level fitting for backspacing in Case 3). This might be a result of constraints in the behavior space that the simulation model can express; for example, some behavioral features reproduced by the model may display inherent mutual dependencies. If the constraints are not aligned with real users’ behavior space, it will be impossible to replicate their behavior features in full. Then, one must find the closest possible point instead.

We were able to infer model parameters for 1,057 individuals by means of amortized inference, whereas the previous highest number was 211, with ABC [18]. The inference process for all users took only 20 seconds in total. This level of efficiency is remarkable and signals that the method scales up well. It affords analysis that captures even smaller effects in datasets and can connect them with parameters of a simulation model. The previously unreported difference we revealed between male and female users’ motor resources in finger movement (k) is a case in point (Figure 14).

The density estimator showed some robustness to distribution shifts between training and test data, especially in dealing with data from unseen task conditions (e.g., the pointing target’s speed in Case 2) or different priors. Still, to achieve better inference performance against empirical data, one needs to hone the priors from the literature and make the training data as comprehensive as possible. Cases of other types can inform testing of density-estimator robustness. Questions remain with regard to, for instance, training data generated by mechanisms that do not perfectly match real-world data.

Our findings suggest that the two estimator types manifest a tradeoff relationship. A point estimator can further reduce the computation cost of inference, yet a density estimator might be more beneficial in several cases: 1) rigorously checking the reliability of the inferred values, 2) identifying multiple plausible explanations for observations, and 3) using the inferred posterior as a new prior to perform multi-round inference that is expected to achieve higher accuracy (e.g., in sequential neural posterior estimation [23, 57]).

Together, the robustness and low computation cost brought by amortized inference offer an exciting vista for future application of simulation models in 1) adaptive user interfaces, 2) recommendation systems, 3) diagnostic tools for user modeling, and 4) large-scale behavior analysis tools.

For the first of these, parameters inferred for a user could form the basis for adapting an interface optimally for the individual. For example, simulation models aid in optimizing a keyboard layout for people with such impairments as dyslexia or tremor [68]. Adaptation of this sort via simulation models has been handled mainly on a non-real-time basis thus far [68, 73, 74]. Our approach could open the door to developing adaptive user interfaces in a new way, enabling repetitive processes of inference from observations and optimization in real time. Secondly, amortized inference could advance interactive recommendation systems through inferring a user’s temporal intention and interest. For instance, future simulation models could assist in better disentangling various underpinnings of click data (individuals’ preference, intention, curiosity, etc.). As for developing user models with theory-based mechanisms, parameter inference enables a proposed model to be evaluated with actual user data, and amortized inference can speed up the investigation. Also, parameter identifiability is often of interest to modelers; the shape of the inferred posterior distribution informs how well each parameter can be distinguished via what is observable. With this approach, one can ensure the validity of datasets in future experiments and identify the number of observations required for valid inferences. Finally, low-cost inverse modeling enables applying simulation models to identify humans’ latent capabilities from datasets orders of magnitude larger than feasible ever before. For example, by inversely modeling an individual’s behavior, one could recover the attributes of the user’s motor and cognitive systems – muscle activation, memory recall, and many more.