2.1 Marginal log-likelihood

Let $f(x|\theta , \boldsymbol {\tau })$ be the probability of observed response x given ability parameter $\theta $ and item parameter vector $\boldsymbol {\tau }$ . For example, when the two-parameter logistic model (2PL; Birnbaum, Reference Birnbaum, Lord and Novick1968) is used, the function becomes

(1) $$ \begin{align} \begin{aligned} P(\theta, \boldsymbol{\tau}) = \Pr(x=1|\theta, \boldsymbol{\tau}) = \frac{\exp(a(\theta-b)}{1+\exp(a(\theta-b))}, \end{aligned} \end{align} $$

(2) $$ \begin{align} \begin{aligned} f(x|\theta, \boldsymbol{\tau}) = P(\theta, \boldsymbol{\tau})^{x}(1-P(\theta, \boldsymbol{\tau}))^{1-x}, \end{aligned} \end{align} $$

where $x\in \{0,1\}$ , $\boldsymbol {\tau }=\left [a,b\right ]^{'}$ , and a and b are item discrimination and difficulty parameters, respectively. Then, the likelihood of the jth test taker ( $j=1,2,\ldots ,N$ ) can be expressed as follows:

(3) $$ \begin{align} \begin{aligned} L_{j}(\theta_{j}, \mathbf{T}|\mathbf{x}_j)=\prod_{i=1}^{I}{f(x_{ji}|\theta_{j}, \boldsymbol{\tau} _{i})}, \end{aligned} \end{align} $$

where i denotes items ( $i=1,2,\ldots ,I$ ), $\mathbf {T}$ is item parameter matrix, and $\mathbf {x}_j=[x_{j1}, \ldots , x_{ji}, \ldots , x_{jI}]^{\text {T}}$ is the vector of the jth test taker’s item responses.

The MML-EM procedure uses a quadrature scheme to facilitate the EM algorithm. Although other quadrature schemes can be adopted, such as adaptive quadrature (e.g., Haberman, Reference Haberman2006), equidistant points are typically used to approximate the integral. For example, we can focus on the $\theta $ range from $-6.05$ to $6.05$ to include the density mass as much as possible when assuming that the mean and standard deviation of the latent density are 0 and 1, respectively. Dividing the range into 121 equal width intervals of 0.1, the quadrature points are then set at the midpoints of these intervals, resulting in points at $-6.0, -5.9, -5.8, \ldots , 6.0$ . Using the latent density $g(\theta , \boldsymbol {\xi })$ with the density parameter vector $\boldsymbol {\xi }$ , the marginal log-likelihood of the N test takers becomes

(4) $$ \begin{align} \begin{aligned} \log{L} &= \log{L(\mathbf{T}, \boldsymbol{\xi}|\mathbf{X})}\\ &= \log{\left\{\prod_{j=1}^{N}\left[ \int_{-\infty}^{\infty}{L_{j}(\theta, \mathbf{T}|\mathbf{x}_j)}g(\theta, \boldsymbol{\xi}) \, d\theta \right]\right\}}\\ &=\sum_{j=1}^{N}\log{\left[ \int_{-\infty}^{\infty}{L_{j}(\theta, \mathbf{T}|\mathbf{x}_j)}g(\theta, \boldsymbol{\xi}) \, d\theta \right]} \\ &\approx \sum_{j=1}^{N}\log{\left[ \sum_{q=1}^{Q}L_{j}(\theta_{q}^{*}, \mathbf{T}|\mathbf{x}_j)A(\theta_{q}^{*}, \boldsymbol{\xi}) \right]}, \end{aligned} \end{align} $$

where $\mathbf {X}$ is the item response matrix with $\mathbf {X}_{j,i}=x_{ji}$ . The subscript j on $\theta $ is discarded by the independent and identically distributed ( $iid$ ) assumption. Approximating the integral by summation using a quadrature scheme, $\theta _{q}^{*}$ denotes the qth quadrature point ( $q=1,2,\ldots ,Q$ ) and $A(\theta _{q}^{*})=\frac {g(\theta _{q}^{*})}{\sum ^{Q}_{q=1}{g(\theta _{q}^{*})}}$ is a normalized latent distribution at $\theta _{q}^{*}$ . The superscript $*$ is introduced to emphasize the discretization of the continuous variable $\theta $ .

2.2 Likelihood of latent density

In the expectation step (E-step) of the MML-EM procedure, the quantity $\gamma _{jq}$ is calculated using Bayes’ theorem (Baker & Kim, Reference Baker and Kim2004):

(5) $$ \begin{align} \gamma_{jq} = \Pr_{j}(\theta^{*} = \theta_{q}^{*}) = \frac{L_{j}(\theta_{q}^{*}, \mathbf{T}|\mathbf{X})A(\theta_{q}^{*}, \boldsymbol{\xi})}{\sum_{q=1}^{Q}L_{j}(\theta_{q}^{*}, \mathbf{T}|\mathbf{X})A(\theta_{q}^{*}, \boldsymbol{\xi})}, \end{align} $$

which is the posterior probability of the jth test taker at $\theta ^{*} = \theta _{q}^{*}$ .

Using the results from the E-step, the marginal likelihood of the M-step can be expressed as follows:

(6) $$ \begin{align} \log{L(\mathbf{T}, \boldsymbol{\xi}|\mathbf{X})} &\approx \sum_{j=1}^{N}\log{\left[ \sum_{q=1}^{Q}L_{j}(\theta_{q}^{*}, \mathbf{T}|\mathbf{X})A(\theta_{q}^{*}, \boldsymbol{\xi}) \right]}\notag\\ &= \sum_{j=1}^{N}\sum_{q=1}^{Q}\gamma_{jq}\log{\left[ L_{j}(\theta_{q}^{*}, \mathbf{T}|\mathbf{X})A(\theta_{q}^{*}, \boldsymbol{\xi}) \right]} - \sum_{j=1}^{N}\sum_{q=1}^{Q}\gamma_{jq}\log{\gamma_{jq}} \notag\\ &= \sum_{j=1}^{N}\sum_{q=1}^{Q}\gamma_{jq}\log{\left[ L_{j}(\theta_{q}^{*}, \mathbf{T}|\mathbf{X}) \right]} + \sum_{j=1}^{N}\sum_{q=1}^{Q}\gamma_{jq}\log{\left[A(\theta_{q}^{*}, \boldsymbol{\xi}) \right]} - \sum_{j=1}^{N}\sum_{q=1}^{Q}\gamma_{jq}\log{\gamma_{jq}} \notag\\ &= \qquad \ell_{item} \qquad + \qquad \ell_{density} \qquad - \qquad \text{constant}\,. \end{align} $$

The separation of the likelihood demonstrates the statistical independence of the item parameters $\mathbf {T}$ and density parameters $\boldsymbol {\xi }$ . Consequently, the density parameters depend only on $\ell _{density}$ .

2.3 Existing methods

Several latent density estimation methods have been proposed to flexibly model the distribution of latent traits in IRT. These methods vary in how they represent the latent distribution, ranging from simple empirical estimates to flexible nonparametric approaches. While these approaches aim to relax the normality assumption underlying the MML estimation, they differ in estimation accuracy, interpretability, computational burden, and availability of model evaluation criteria.

EHM (Bock & Aitkin, Reference Bock and Aitkin1981; Mislevy, Reference Mislevy1984) is a nonparametric approach that directly utilizes the expected frequencies calculated in the E-step of the EM algorithm. Specifically, it estimates the latent density at each quadrature point as the normalized sum of posterior probabilities across respondents. This method is computationally efficient and easy to implement, since results from the E-step can be directly transformed into a density estimate without additional computation. However, because it does not apply any smoothing, EHM is sensitive to sampling noise, which can result in erratic density shapes and suboptimal item parameter estimates (Li, Reference Li2021; Woods, Reference Woods, Reise and Revicki2015; Woods & Lin, Reference Woods and Lin2009).

DCM (Woods & Lin, Reference Woods and Lin2009) and RCM (Woods, Reference Woods2006) both aim to improve upon EHM by smoothing the noisy expected frequencies. These methods differ in the functional forms they employ for density estimation, but share a common estimation strategy and model selection procedure. DCM models the latent distribution as a squared polynomial multiplied by a standard normal density. The order of the polynomial, denoted by h, determines the complexity of the distribution: when $h = 1$ , the model reduces to the normal distribution, while larger values of h allow for increasingly flexible density shapes. To ensure that the resulting function is a proper probability distribution, the polynomial coefficients are constrained via a polar coordinate transformation. On the other hand, RCM uses basis splines, where the order of the splines and the number of breaks are the hyperparameters that control the smoothness of the density. In the estimation of density parameters of RCM, it employs a normal prior to improve estimation stability. However, the standard deviation of the prior is often chosen through trial and error. Both methods require the selection of smoothing parameter(s), which is typically achieved by fitting multiple models and selecting the best one based on the Hannan–Quinn (HQ) criterion (Hannan & Quinn, Reference Hannan and Quinn1979). This criterion, which is similar to information criteria, such as the Akaike information criterion (AIC) and the Bayesian information criterion (BIC), balances model flexibility and overfitting. For RCM, however, visual inspection of the candidate models is also necessary in addition to the HQ criterion.

LLS (Casabianca & Lewis, Reference Casabianca and Lewis2015) is another approach that aims to regularize the noisy frequency estimates from EHM by fitting a log-linear model to the expected posterior counts. It expands the log-density as a linear combination of polynomials, and its hyperparameter controls the smoothness of the density by determining the highest degree of the polynomials. However, unlike DCM, the model selection strategy of LLS has not been clearly established to determine the optimal hyperparameter.

The 2NM (Li, Reference Li2021; Mislevy, Reference Mislevy1984) is a parametric alternative that extends the standard normality assumption by modeling the latent distribution as a weighted sum of two normal distributions. This added flexibility allows 2NM to capture deviations from normality, such as skewness or bimodality. The density parameters are estimated using a secondary EM algorithm nested within the MML-EM framework. While the parametric form of 2NM offers advantages, such as interpretability and analytic simplicity, it also imposes constraints on flexibility, making it less capable of capturing complex or highly irregular latent distributions compared to nonparametric or semi-nonparametric approaches.

3.1 Kernel density estimation

KDE is a nonparametric density estimation method that assigns a kernel function to each observed data point (Silverman, Reference Silverman1986). Let the vector $\mathbf {z}=\left [Z_{1}, \, Z_{2}, \, \dots , \, Z_{N}\right ]^{\text {T}}$ denote N observed data points. The kernel functions formulate the density as follows:

(7) $$ \begin{align} g(z|h)= \frac{1}{Nh}\sum_{j=1}^{N}{K\left(\frac{z-Z_{j}}{h}\right)}, \end{align} $$

where $K(\cdot )$ is the kernel function and h is the hyperparameter called bandwidth. The bandwidth controls the degree of density smoothing by determining the width of the kernel function. For example, when the Gaussian kernel is applied, the bandwidth h represents the standard deviations of the kernels. Among the available kernel functions, such as Epanechnikov, biweight, triangular, Gaussian, and uniform, the Gaussian kernel is most commonly adopted due to the minimal influence of kernel choice on the resulting density estimates and the differentiability of the Gaussian form (Gramacki, Reference Gramacki2018; Silverman, Reference Silverman1986). Following this convention, the present study focuses on the Gaussian kernel and its applications.

In Figure 1, for illustrative purposes, the standard normal distribution (dotted lines) is estimated using four different bandwidth values (i.e., 0.4, 0.7, 1.0, and 1.3). Using the seven data points generated from the standard normal distribution, seven kernel functions (dashed lines) are formulated. Mathematically, kernel functions are Gaussian distributions with their means as the seven data points and a common standard deviation of h, and are scaled by the total sample size: $\frac {1}{N}=\frac {1}{7}$ . The estimated densities (solid lines) are summations of the kernel functions. Among the four density estimates, the standard normal distribution is best recovered when $h=0.7$ . In comparison, the function is under-smoothed when $h=0.4$ , resulting in a bimodal distribution. When $h=1.0$ and $h=1.3$ , the flatter shapes compared to the standard normal distribution indicate that they were over-smoothed. Considering that the optimal bandwidth is $h=0.7177$ for this example, it is shown that the density estimate is most accurate when the bandwidth h was closest to the optimal value (see Gramacki, Reference Gramacki2018 for the formula for calculating the optimal bandwidth when the true density is known).

Figure 1

Estimated densities according to different h values.

Importantly, the bandwidth h of KDE determines the degree of smoothing, without necessitating an estimation of additional distribution parameters. This characteristic of KDE obviates the need to include multiple model-fitting procedures and a model selection process to accurately estimate the latent density.

3.2 Bandwidth estimation and approximate mean integrated square error (AMISE)

As bandwidth estimation (or optimal bandwidth selection) itself becomes the density estimation of KDE, the search for an optimal bandwidth is crucial to obtain accurate density estimates. To this end, existing literature on density estimation has focused on the AMISE criterion for bandwidth estimation (see Gramacki, Reference Gramacki2018; Silverman, Reference Silverman1986; Wand & Jones, Reference Wand and Jones1995). Consequently, within this context, the optimal bandwidth can be defined as the one that minimizes the AMISE criterion.

AMISE is the approximation of the mean integrated squared error (MISE) obtained from a Taylor expansion (Jones, Reference Jones1991; Sheather, Reference Sheather2004; Silverman, Reference Silverman1986; Wand & Jones, Reference Wand and Jones1995), where $\text {MISE}(\hat {g})=E\left [\int {\left \{\hat {g}(\theta )-g(\theta )\right \}^{2}\, dx}\right ]$ is an appropriate index for evaluating the density estimation results (Jones, Reference Jones1991). AMISE can be expressed as

(8) $$ \begin{align} \text{AMISE}(\hat{g}_{h})= \frac{1}{Nh}R(K)+\frac{h^4}{4}\mu_2(K)^2R(g"), \end{align} $$

where $\hat {g}_{h}$ is the density estimate using the bandwidth h, $R(K)=\int {K(t)^2\, dt}$ , $\mu _2(K)=\int {t^2K(t)\, dt}$ , and $R(g")=\int {g"(t)^2\, dt}$ . When the Gaussian kernel is used for the kernel K, it can be obtained from a simple algebra that $R(K)=\left (2\sqrt {\pi }\right )^{-1}$ and $\mu _2(K)=1$ .

Another equation can be obtained by differentiating Equation (8) with respect to h and equating it to $0$ :

(9) $$ \begin{align} \hat{h}_{\text{AMISE}}=\left[\frac{R(K)}{N\mu_2(K)^2\psi_4}\right]^{\frac{1}{5}}, \end{align} $$

where $\hat {h}_{\text {AMISE}}$ denotes the bandwidth estimate from the AMISE criterion, $\psi _d=E_{t}\left [g^{(d)}(t)\right ]$ , and $\psi _4=\int {g^{(4)}(t)g(t)\, dt}=\int {g"(t)^2\, dt}=R(g")$ . Unfortunately, Equation (9) is not a closed-form solution to estimate h as it depends on the information from an unknown density g, which is $\psi _4=R(g")$ . Therefore, one of the widely used strategies is to calculate the estimate $\hat {\psi }_4$ and plug it into Equation (9), and consequently, the accurate estimation of $\psi _4$ determines the overall performance of KDE (Gramacki, Reference Gramacki2018; Sheather, Reference Sheather2004; Sheather & Jones, Reference Sheather and Jones1991; Wand & Jones, Reference Wand and Jones1995).

3.3 Plug-in method

Among the three types of bandwidth selection methods—rule-of-thumb, cross-validation, and plug-in—the plug-in method is known to produce a bandwidth estimator that has an efficient convergence rate to the parameter and small bias and variance (Gramacki, Reference Gramacki2018; Harpole et al., Reference Harpole, Woods, Rodebaugh, Levinson and Lenze2014; Sheather, Reference Sheather2004; Sheather & Jones, Reference Sheather and Jones1991; Wand & Jones, Reference Wand and Jones1995). Due to its precision and efficiency in bandwidth estimation, the plug-in method has been a recommended option for density estimation in statistical software programs, such as R and SAS (R Core Team, 2025; SAS Institute Inc., 2020).

Direct plug-in (DPI) and solve-the-equation (STE) are two types of plug-in methods. Within the overall scheme of the $\psi _{4}$ estimation, both plug-in methods utilize $\psi _{d+2}$ in estimating $\psi _{d}$ , where an even number d is the order of differentiation and $p_d$ is a pilot bandwidth introduced to estimate $\psi _{d}=E\left (g^{(d)}(x)\right )$ (Gramacki, Reference Gramacki2018; Sheather, Reference Sheather2004; Wand & Jones, Reference Wand and Jones1995). That is, we start the algorithm by setting d to a specific value (e.g., $d=8$ ), then iteratively estimate $\psi _{d-2}$ using $\psi _{d}$ until we get the estimate of $\psi _4$ . In this study, we focus on STE (Sheather & Jones, Reference Sheather and Jones1991) between the two, which is the most widely used method in unidimensional density estimation due to its small bias and variance in bandwidth estimation (Gramacki, Reference Gramacki2018; Harpole et al., Reference Harpole, Woods, Rodebaugh, Levinson and Lenze2014; Jones et al., Reference Jones, Marron and Sheather1996; Sheather, Reference Sheather2004; Sheather & Jones, Reference Sheather and Jones1991). STE is a suggested option in R (R Core Team, 2025) and is set as the default in SAS (SAS Institute Inc., 2020). Note that, at the time of writing this article, STE is not fully extended to multidimensional cases. In comparison, DPI showed more accurate and stable performance in estimating the multidimensional bandwidth compared to rule-of-thumb and cross-validation (Chacón & Duong, Reference Chacón and Duong2010; Duong & Hazelton, Reference Duong and Hazelton2003).

To implement STE, some additional equations are used:

(10) $$ \begin{align} p_{d} = \left[ \frac{2K^{(d)}(0)}{-\mu_{2}(K)\psi_{d+2}(p_{d+2})N} \right]^{\frac{1}{d+3}}, \end{align} $$

(11) $$ \begin{align} \psi_{d}^{N} = \left[ \frac{(-1)^{d/2}d!}{(2\hat{\sigma})^{d+1}\left(\frac{d}{2}\right)!\pi^{1/2}} \right]^{1/9}, \end{align} $$

and

(12) $$ \begin{align} p(h) = \left[ \frac{2K^{(4)}(0)\mu_2(K)\hat{\psi}_4(p_4)}{-\hat{\psi}_6(p_6)R(K)} \right]^{1/7}h^{5/7}. \end{align} $$

Above, Equation (10) provides the unbiased estimator of $p_d$ using the information of $\psi _{d+2}$ , Equation (11) is the formula to compute $\psi _d$ assuming that the distribution is normal, and Equation (12) expresses a pilot bandwidth $p_4$ as a function of h (Gramacki, Reference Gramacki2018; Sheather & Jones, Reference Sheather and Jones1991).

The algorithm can be summarized as follows (Gramacki, Reference Gramacki2018):

1. Step 1: Calculate $\psi _{6}^{N}$ and $\psi _{8}^{N}$ from Equation (11).

2. Step 2: Calculate $p_4$ and $p_6$ from Equation (10) using $\psi _{6}^{N}$ and $\psi _{8}^{N}$ .

3. Step 3: Estimate $\psi _{4}(p_4)$ and $\psi _{6}(p_6)$ using $p_4$ and $p_6$ , respectively.

4. Step 4: Estimate $\psi _{4}$ using the estimator $\hat {\psi }\left (p(h)\right )$ , where $p(h)$ is calculated from Equation (12) using $\hat {\psi }_4(p_4)$ and $\hat {\psi }_6(p_6)$ obtained from Step 3.

5. Step 5: The estimated bandwidth $\hat {h}$ is the numerical solution of Equation (9) when $\hat {\psi }\left (p(h)\right )$ is plugged into the equation.

Note that the illustrated algorithm above is the two-step algorithm where “two-step” refers to the process of starting from $d=8$ and terminating at $d=4$ . In the same way, the three-step algorithm starts from $d=10$ , then sequentially moves to $8$ , $6$ , and $4$ . Although increasing the number of steps reduces bias, the two-step algorithm is generally preferred due to the trade-off between bias and variance.

3.4 Application of KDE to IRT

As with other IRT latent density estimation methods, KDE can be applied to the latent density estimation of IRT by treating $\hat {f}_{q}=\sum _{j=1}^{N}\gamma _{jq}$ as observed data points, where $\hat {f}_{q}$ is the expected frequency of test takers at $\theta ^{*}=\theta ^{*}_{q}$ . For the ease of programming, $\hat {f}$ is rounded in this study. Suppose, for illustration, that we are in the ith iteration of the MML-EM procedure, where “i” temporarily denotes the sequence of the EM iteration, not related to the subscript i used for items. Using $\hat {f}_{q}$ obtained from the E-step of the ith EM iteration, the STE algorithm in Section 3.3 is implemented to get $\hat {h}$ during the M-step of the ith EM iteration. The density estimate $g(\theta |\hat {h})$ (see Equation (7)) is first scaled to have a mean of 0 and a standard deviation of 1 using linear interpolation and extrapolation, which is also applied to DCM to identify the latent scale (Woods & Lin, Reference Woods and Lin2009). Then, it becomes the estimated latent distribution of the ith EM iteration, which, in turn, is used in the E-step of the $(i+1)$ th EM iteration. When the MML-EM procedure converges, $g(\theta |\hat {h})$ from the last EM iteration is taken as the estimated latent distribution.

Unlike other methods that employ additional density parameters according to the hyperparameter h, only the bandwidth h is involved during the application of KDE. As a result, KDE can yield a smoothed latent density estimate within a single MML-EM procedure.

4.1 Compared methods

To evaluate the accuracy of the KDE method in parameter estimation using computer simulation techniques, the performance of KDE is assessed under various simulation conditions by comparing it with three existing approaches: the normality-assumption method (NM), EHM, and DCM.

NM refers to the conventional approach that does not estimate the latent distribution but instead assumes normality. It serves as a baseline model—similar to a null model—against which the advantages of more flexible approaches can be evaluated. In a strict sense, NM differs from the other approaches because it does not attempt to estimate the latent density.

EHM provides a baseline among density estimation methods, as it relies on normalized histograms without smoothing. DCM further smooths the non-smoothed density estimate produced by EHM at every EM iteration. DCM has shown competitive accuracy among existing methods and can offer a meaningful benchmark for assessing new techniques (Woods & Lin, Reference Woods and Lin2009). In addition, its principled model selection procedure based on the HQ criterion makes it particularly suitable for fair comparisons.

Other methods, such as RCM, LLS, and 2NM, are omitted from the comparison due to challenges in evaluating their performance on equal grounds. RCM and LLS both involve subjective decisions in the model selection process: RCM requires trial and error in setting the standard deviation of the prior distribution and often relies on visual inspection, while LLS lacks a clear model selection guideline. The 2NM method depends heavily on the choice of the true density in the simulation design due to its parametric nature. This dependency is misaligned with the focus of the present study on density estimation methods applicable when little or no information about the true density is available, which is the motivation for adopting nonparametric approaches.

Accordingly, NM, EHM, and DCM are chosen to represent the null, lower, and upper bounds of existing approaches in terms of flexibility and performance.

4.2 Simulation design

The simulation conditions are formulated by three factors that can potentially affect the results: the shape of the latent distribution, the length of the tests, and the number of test takers. Figure 2 shows the four shapes of the latent distribution used for the simulation study: normal, skewed, bimodal, and skewed bimodal. The three non-normal distributions are formulated by a mixture of two normal components, and their mean and variance are $0$ and $1$ , respectively (see Li, Reference Li2021 for reparameterization of the two-component normal mixture distribution to fix the overall mean and variance). Setting the standard normal distribution as a reference, the skewed and bimodal distributions represent the two types of normality violation, and the skewed bimodal distribution represents the combination of the two sources of violation. The skewness values of the normal and bimodal distributions are both 0. In contrast, the skewed and skewed bimodal distributions have skewness values of approximately 0.6 and 0.5, respectively. These values were chosen to approximate the skewness observed in the IRT scale score distributions of the U.S. state-level testing programs (Ho & Yu, Reference Ho and Yu2015). For example, the median skewness of the 24 IRT scale score distributions in Colorado was approximately $-0.6$ , whereas in New York, it was approximately $0.7$ , ignoring the sign of skewness. The shape of the skewed bimodal distribution is designed to resemble the example of bimodal test score distributions presented in Sibbald (Reference Sibbald2014). Although a symmetric bimodal distribution may not frequently occur in practice as a latent distribution, it is included in the simulation conditions to examine the effect of bimodality independently of skewness. Test lengths of 8, 15, 30, and 60 items can be viewed as representing very short, short, medium, and long tests, respectively. Very short or short tests include brief scales, such as the 8-item version of the Center for Epidemiologic Studies Depression (CES-D8) scale (Missinne et al., Reference Missinne, Vandeviver, Van de Velde and Bracke2014) and end-of-chapter quizzes in online learning platforms, where the latent distribution has a relatively strong influence on parameter estimation. Examples of long tests are large-scale educational assessments, such as the Scholastic Aptitude Test and the American College Testing. Lastly, the number of test takers is varied—set at 500, 1,000, and 2,000—to examine the performance of the methods with respect to sample size.

Figure 2

Normal, skewed, bimodal, and skewed bimodal distributions used for the simulation study.

For each of the 48 simulation conditions (4 distributions $\times $ 4 test lengths $\times $ 3 sample sizes), 200 item response datasets are generated using the 2PL (Birnbaum, Reference Birnbaum, Lord and Novick1968) and fitted by NM, EHM, DCM, and KDE. For the random sampling of item parameters, the a parameter is generated from a uniform distribution from $0.8$ to $2.5$ , and the b parameter is generated from the standard normal distribution truncated at $\pm 2$ . These parameters serve as true values in calculating the evaluation criteria in Equations (14) and (15). The IRTest package (Li, Reference Li2024; R Core Team, 2025)—an R package specialized in latent density estimation of IRT—is used for both data generation and model fitting. The source code for the simulation study is available at https://github.com/SeewooLi/KDE_IRT_code. The KDE option in IRTest has been developed particularly for this study. An MML-EM procedure is considered to be successfully converged when the maximum parameter change drops below 0.0001 within the specified maximum number of EM iterations. The maximum number of iterations is set to 200 for NM, DCM, and KDE, and to 2,000 for EHM to accommodate its slower convergence. These settings are consistent with previous studies that used similar criteria to determine convergence of the MML-EM procedure (Woods, Reference Woods2006; Woods & Lin, Reference Woods and Lin2009). An AMD Ryzen 7 5700G processor is used for the study.

The accuracy of the density estimation is evaluated by the integrated squared error (ISE). In particular, ISE from NM can indicate the degree to which the true density deviates from the normal distribution, because NM assumes the latent distribution to be normal and, thus, does not estimate it. The ISE criterion can be written as follows:

(13) $$ \begin{align} \begin{aligned} \text{ISE}(\hat{g})&=\int^{\infty}_{-\infty}{\left[ \hat{g}(\theta) - g(\theta) \right]^2\,d\theta}\\ &\approx \sum^{121}_{q=1}{\frac{1}{10}\left[ \hat{g}(\theta^{*}_{q}) - g(\theta^{*}_{q}) \right]^2}. \end{aligned} \end{align} $$

The integration is approximated by the quadrature scheme of using 121 equally spaced grids from $-6$ to 6.

Lastly, the root mean square error (RMSE) evaluates the accuracy of the estimation of the ability parameter $\theta $ , item characteristic curve (ICC), and the item parameters (a and b):

(14) $$ \begin{align} \begin{aligned} \text{RMSE}(\hat{\text{ICC}})=\sqrt{\frac{1}{I}\sum^{I}_{i=1}{\sum^{121}_{q=1}{ \left(\hat{P}_{i}(\theta^{*}_{q}) - P_{i}(\theta^{*}_{q})\right)^2 \hat{A}(\theta^{*}_{q}) }}}, \end{aligned} \end{align} $$

(15) $$ \begin{align} \begin{aligned} \text{RMSE}(\hat{\tau})=\sqrt{\frac{1}{I}\sum^{I}_{i=1}{\left(\hat{\tau} - \tau\right)^2}}, \end{aligned} \end{align} $$

and

(16) $$ \begin{align} \begin{aligned} \text{RMSE}(\hat{\theta})=\sqrt{\frac{1}{N}\sum^{N}_{i=1}{\left(\hat{\theta} - \theta\right)^2}}, \end{aligned} \end{align} $$

where $\hat {P}_{i}(\theta ^{*}_{q})$ in Equation (1) is the estimated ICC of Item i evaluated at $\theta ^{*}_{q}$ , $\hat {A}(\theta ^{*}_{q})$ is the estimated normalized density (see Section 2.1) at the qth grid, and $\tau $ represents either the a or b parameter. Expected a priori (EAP) scores are used for $\hat {\theta }$ . Based on 200 data replications, the average RMSE and ISE values for each method are reported in Section 4.3.