A Semiparametric Gaussian Mixture Model for Chest CT-based 3D Blood Vessel Reconstruction

Let $(Y_{i},X_{i})$ be the observation collected from the $i$th subject, with $Y_{i}\in\mathbb{R}^{1}$ being the univariate response of interest and $X_{i}\in\mathbb{D}=[0,1]^{3}$ being the associated voxel position in a three-dimensional Euclidean space. We assume that $X_{i}$ is uniformly distributed on $\mathbb{D}$. Furthermore, we assume that for each $i$ a latent class label $Z_{i}\in\{1,\cdots,M\}$, where $M$ is the total number of classes. We then assume that $P(Z_{i}=m|X_{i})=\pi_{m}(X_{i})$ for $1\leq m\leq M$. Clearly, we should have $\sum_{m=1}^{M}\pi_{m}(x)=1$ for any $x\in\mathbb{D}$. This suggests that we can define $\pi_{M}(x)=1-\sum_{m=1}^{M-1}\pi_{m}(x)$ for any $x\in\mathbb{D}$ throughout this article. Conditional on $Z_{i}=m$ and $X_{i}$, we assume that $Y_{i}$ is normally distributed with a mean $\mu_{m}(X_{i})$ and variance $\sigma_{m}^{2}(X_{i})$. Here, we assume that class prior probability function $\pi_{m}(x)$, mean function $\mu_{m}(x)$, and variance function $\sigma_{m}^{2}(x)$ are all smooth functions in $x\in\mathbb{D}$ and vary between classes. Next, we consider how to consistently estimate them.

We collect the observed voxel positions and CT values into $\mathbb{X}=\{X_{i}:1\leq i\leq N\}$ and $\mathbb{Y}=\{Y_{i}:1\leq i\leq N\}$, respectively. Furthermore, we collect the unobserved latent class labels into $\mathbb{Z}=\{Z_{i}:1\leq i\leq N\}$. For a given voxel position $x$, define $\theta(x)=\{\pi^{\top}(x),\mu^{\top}(x),\sigma^{\top}(x)\}^{\top}\in\mathbb{R}^{3M-1}$, where $\pi(x)=\{\pi_{1}(x),\cdots,\pi_{M-1}(x)\}^{\top}\in\mathbb{R}^{M-1}$, $\mu(x)=\{\mu_{1}(x),\cdots,\mu_{M}(x)\}^{\top}\in\mathbb{R}^{M}$, and $\sigma(x)=\{\sigma_{1}(x),\cdots,\sigma_{M}(x)\}^{\top}\in\mathbb{R}^{M}$. Define $\Theta=\{\theta(x):x\in\mathbb{D}\}$. To estimate $\theta(x)$, we develop a novel KEM method. We begin with a highly simplified case with $\pi_{m}(x)=\pi_{m}$, $\mu_{m}(x)=\mu_{m}$, and $\sigma_{m}^{2}(x)=\sigma_{m}^{2}$ for some constants $\pi_{m}>0$, $\mu_{m}$, and $0<\sigma_{m}<+\infty$. We then have the log-likelihood function for an interior point $x$ in $\mathbb{D}$ as

where $f(x,y|\Theta)$ stands for the joint probability density function of $(X_{i},Y_{i})$ evaluated at $(X_{i},Y_{i})=(x,y)$; and $f(y|m,X_{i},\Theta)$ is the marginal probability density function of $Y_{i}$ evaluated at $Y_{i}=y$ conditional on $Z_{i}=m$ and $X_{i}$. Moreover, $f(x|\Theta)$ is the probability density function of $X_{i}$ evaluated at $X_{i}=x$. Since $X_{i}$ is assumed to be uniformly generated on $\mathbb{D}$, we have $f(x|\Theta)=1$ for $x\in\mathbb{D}$ and $f(x|\Theta)=0$ for $x\not\in\mathbb{D}$. Moreover, the function $\phi(y)=\exp(-y^{2}/2)/\sqrt{2\pi}$ is the probability density function of a standard normal random variable.

Nevertheless, class prior probability function $\pi_{m}(x)$, mean function $\mu_{m}(x)$, and variance function $\sigma_{m}^{2}(x)$ in our case are not constant. Instead, they should vary in response to the value of $x$ (i.e., the voxel position within the 3D CT space) in a fully nonparametric manner. However, for an arbitrarily given $x$ position, we should expect $\pi_{m}(x)>0$, $\mu_{m}(x)$, and $0<\sigma_{m}(x)<+\infty$ to be locally approximately constant. This is indeed a reasonable assumption as long as these functions are sufficiently smooth with continuous second-order derivatives. Thus, according to Proposition 3.10 of Lang (2012), we then know that the second-order derivatives of $\pi_{m}(x)$, $\mu_{m}(x)$, and $\sigma_{m}(x)$ are uniformly upper bounded from infinity on the compact set $\mathbb{D}$ for any $1\leq m\leq M$. Similarly, we know that $\sigma_{m}(x)$ is uniformly lower bounded on $\mathbb{D}$ by a constant $\sigma_{\min}$ for $1\leq m\leq M$. Therefore, we draw inspiration from the concept of the local constant estimator (Nadaraya, 1965; Watson, 1964) and local maximum likelihood estimation (Fan et al., 1998) and introduce the subsequent locally weighted log-likelihood function based on the observed data $(\mathbb{X},\mathbb{Y})$:

where $x=(x_{1},x_{2},x_{3})^{\top}$ stands for a fixed interior point in $\mathbb{D}$. Note that $\theta=(\pi^{\top},\mu^{\top},\sigma^{\top})^{\top}\in\mathbb{R}^{3M-1}$ is a set of working parameters with $\pi=(\pi_{1},\cdots,\pi_{M-1})^{\top}\in\mathbb{R}^{M-1}$, $\mu=(\mu_{1},\cdots,\mu_{M})^{\top}$ $\in\mathbb{R}^{M}$, and $\sigma=(\sigma_{1},\cdots,\sigma_{M})^{\top}\in\mathbb{R}^{M}$. For class $M$, we set $\pi_{M}=1-\sum_{m=1}^{M-1}\pi_{m}$. Moreover, the three-dimensional kernel function, $\mathbb{K}(\cdot)$, is assumed to be $\mathbb{K}(x)=\mathcal{K}(x_{1})\mathcal{K}(x_{2})\mathcal{K}(x_{3})$, where $\mathcal{K}(t)$ with $t\in\mathbb{R}^{1}$ is a continuous probability density function symmetric about 0. Here, $h>0$ is the associated bandwidth. What distinguishes (2) from the classical log-likelihood function (1) is that information on the peer positions of $x$ is also blended in for local estimations. For convenience, we refer to $\mathcal{L}_{x}(\theta)$ as a locally weighted log-likelihood function. Then, the local maximum likelihood estimators can be defined as $\hat{\theta}(x)=\mathop{\mbox{argmax}}\limits_{\theta}\mathcal{L}_{x}(\theta)$. For convenience, we refer to $\hat{\theta}(x)$ as the kernel maximum likelihood estimators (KMLE). We next consider how to optimize $\mathcal{L}_{x}(\theta)$ so that $\hat{\theta}(x)$ can be computed.

Recall that the locally weighted log-likelihood function based on the observed data $(\mathbb{X},\mathbb{Y})$ is already defined in (2). Ever since Dempster et al. (1977), a typical way to optimize (2) is to develop an EM algorithm based on the so-called complete log-likelihood function. That is the log-likelihood function obtained by assuming that the latent class label $Z_{i}$ is observed. Thus, if we have access to the complete data $(\mathbb{X},\mathbb{Y},\mathbb{Z})$, we can define a complete log-likelihood function for the given voxel position $x$ as follows:

where $I(Z_{i}=m)$ is an indicator function. In theory, any location of interest can be taken as $x$. In practice, a location of interest is often taken at the recorded voxel positions of CT data. The objective here is to consistently estimate the nonparametric functions $\pi_{m}(x)$, $\mu_{m}(x)$, and $\sigma_{m}(x)$ for $1\leq m\leq M$. We start with a set of initial estimators as $\hat{\pi}_{m}^{(0)}(x)=1/M$ and $\hat{\mu}_{m}^{(0)}(x)=\hat{\mu}_{m}$, where $\hat{\mu}_{m}$ is an initial estimator obtained by (for example) a standard $k$-means algorithm (MacQueen et al., 1967). Furthermore, we set $\hat{\sigma}_{m}^{(0)}(x)=\sigma^{(0)}$, where $\sigma^{(0)}$ can be some pre-specified constant. We write $\hat{\pi}_{m}^{(t)}(x)$, $\hat{\mu}_{m}^{(t)}(x)$, and $\hat{\sigma}_{m}^{(t)}(x)$ as the estimators obtained in the $t$th step.

Next, we study the asymptotic properties of KMLE $\hat{\theta}(x)$. First, we define some notations. Recall that $\mathcal{L}_{x}(\theta)$ is the locally weighted log-likelihood function defined in (2). We define $\nu_{k,m}=\int t^{k}\mathcal{K}^{m}(t)\mathrm{d}t$ for $0\leq k\leq 2$ and $1\leq m\leq 2$. The following technical conditions are then required:

• (C1)

  (Kernel Function) Assume $|\nu_{k,m}|<+\infty$ for $1\leq m\leq 2$ and $0\leq k\leq 2+\delta$ with some $\delta>0$.

• (C2)

  (Bandwidth and Sample Size) Assume $h=C_{h}N^{-1/7}$ for some constant $C_{h}>0$.

As we mentioned before, we assume that kernel function $\mathbb{K}(t)$ is the product of three univariate kernel density functions symmetric about 0. By Condition (C1), each univariate kernel density function should be a well-behaved probability density function with various finite moments. By Condition (C2), we require the bandwidth $h$ to converge to zero at the speed of $N^{-1/7}$ as the sample size, $N$, goes to infinity. As a consequence, the locally effective sample size, denoted by $Nh^{3}$, diverges toward infinity as $N\rightarrow+\infty$. In the meanwhile, the bandwidth $h$ is well constructed so that the resulting estimation bias should not be too large. Both Conditions (C1) and (C2) are fairly standard conditions, which have been widely used in the literature (Fan and Gijbels, 1996; Ullah and Pagan, 1999; Li and Racine, 2007; Silverman, 2018).

We define $I\{\theta(x)\}=\int\partial\ln f(x,y|\Theta)/\partial\theta(x)\times\partial\ln f(x,y|\Theta)/{\partial\theta(x)^{\top}}\times f(x,y|\Theta)\mathrm{d}y$. Denote the Euclidean norm as $\|z\|=\sqrt{z^{\top}z}$ for any vector $z$. Then, we have the following Theorem 2.1, with the technical details provided in Appendix A and Appendix B of the Supplementary Materials.

To ensure the asymptotic properties of the KMLE, we must carefully select the optimal bandwidth $h$. From Condition (C2), we know that the optimal bandwidth should satisfy $h=C_{h}N^{-1/7}$. The question is how to make an optimal choice for $C_{h}$. To do so, an appropriately defined optimality criterion is required. One natural criterion could be out-of-sample forecasting accuracy. Let $(X^{*},Y^{*})$ be an independent copy of $(X_{i},Y_{i})$. Recall that $\hat{\theta}(x)=\{\hat{\pi}_{1}(x),\cdots,\hat{\pi}_{M-1}(x),\hat{\mu}_{1}(x),\cdots,$ $\hat{\mu}_{M}(x),\hat{\sigma}_{1}(x),\cdots,\hat{\sigma}_{M}(x)\}^{\top}$ are the estimators obtained from data $\{(X_{i},Y_{i}):1\leq i\leq N\}$. Note that $E(Y^{*}|X^{*},\Theta)=\sum_{m=1}^{M}\pi_{m}(X^{*})\mu_{m}(X^{*})$. Therefore, a natural prediction for $Y^{*}$ can be constructed as $\hat{Y}^{*}=\sum_{m=1}^{M}\hat{\pi}_{m}(X^{*})\hat{\mu}_{m}(X^{*})$. Its square prediction error (SPE) can then be evaluated as $E(Y^{*}-\hat{Y}^{*})^{2}$. Using a standard Taylor expansion-type argument, we can obtain an analytical formula for SPE using the following theorem:

The technical details of Theorem 2.2 are provided in Appendix C in the Supplementary Materials. By optimizing the leading term of $E(\hat{Y}^{*}-Y^{*})^{2}$ from Theorem 2.2, we know that the optimal bandwidth constant is given by $C_{h}^{*}=(3C_{2}/C_{1})^{1/7}$. However, both the critical constants, $C_{1}$ and $C_{2}$, are extremely difficult to compute without explicit assumptions on the concrete forms of $\pi(x)$, $\mu(x)$, and $\sigma(x)$ with respect to $x$. To solve this problem, a cross-validation (CV) type method is used to compute the optimal bandwidth. To this end, we need to randomly partition all voxel positions into two parts. The first part contains approximately 80% of the total voxels used for training. The remaining 20% is used for testing. For convenience, the indices of the voxels in the training and testing dataset are collected as $\mathcal{I}_{0}$ and $\mathcal{I}_{1}$, respectively. Next, for an arbitrary testing sample, $i\in\mathcal{I}_{1}$, we have $E(Y_{i}|X_{i})=\sum_{m=1}^{M}\pi_{m}(X_{i})\mu_{m}(X_{i})$. Therefore, we are inspired to predict the $Y_{i}$ value using $\hat{Y}_{i}=\sum_{m=1}^{M}\hat{\pi}_{m}(X_{i})\hat{\mu}_{m}(X_{i})$, where the estimates $\hat{\pi}_{m}(X_{i})$ and $\hat{\mu}_{m}(X_{i})$ are computed based on the training data with a given bandwidth constant. Let $\mathbb{C}_{\text{CV}}=\{C_{h}^{(g)}:1\leq g\leq G_{\text{CV}}\}$ be a set of tentatively selected pilot-bandwidth constants, where $G_{\text{CV}}$ is the total number of pilot-bandwidth constants. Therefore, an estimator for SPE can be constructed as $\hat{\mathscr{L}}_{\text{SPE}}(C_{h}^{(g)})=|\mathcal{I}_{1}|^{-1}\sum_{i\in\mathcal{I}_{1}}(Y_{i}-\hat{Y}_{i})^{2}$. A CV-based estimator for the optimal bandwidth constant can then be defined as $\hat{C}_{h}^{\text{CV}}=\mathop{\mbox{argmin}}\limits_{C_{h}^{(g)}\in\mathbb{C}_{\text{CV}}}\hat{\mathscr{L}}_{\text{SPE}}(C_{h}^{(g)})$; that is, the bandwidth constant in $\mathbb{C}_{\text{CV}}$ with the lowest SPE value is chosen. Consequently, a CV-based estimator for the optimal bandwidth is given by $\hat{h}^{\text{CV}}=\hat{C}_{h}^{\text{CV}}\times N^{-1/7}$.

Even though the above CV idea is intuitive, its practical computation is expensive. This is because, for an accurate estimation of the optimal bandwidth constant $C_{h}^{*}$, we typically need to evaluate a large number of pilot-bandwidth constants. To reduce the computation cost, we develop a novel regression method as follows. From Theorem 2.2, we know that $E\{\hat{\mathscr{L}}_{\text{SPE}}(C_{h})\}\approx\sigma_{y}^{2}+(C_{h}^{4}C_{1}/4+C_{h}^{-3}C_{2})N^{-4/7}$. Note that $E\{\hat{\mathscr{L}}_{\text{SPE}}(C_{h})\}$ is approximately a linear function in both $C_{1}$ and $C_{2}$. The corresponding weights are given by $C_{h}^{4}N^{-4/7}/4$ and $C_{h}^{-3}N^{-4/7}$, where $N$ is the given total sample size, and $C_{h}$ is the tentatively selected pilot-bandwidth constant. This immediately suggests an interesting regression-based method for estimating the two critical constants, $C_{1}$ and $C_{2}$, as follows. We define $\mathbb{C}_{\text{REG}}=\{C_{h}^{(g)}:1\leq g\leq G_{\text{REG}}\}$ as a set of carefully selected pilot-bandwidth constants. Note that $G_{\text{REG}}$ determines the total number of pilot-bandwidth constants to be tested. It must not be too large so that the computation cost can be reduced. For example, we fix $G_{\text{CV}}=25$ and $G_{\text{REG}}=5$ in our subsequent numerical studies. We define a pseudo response, $\mathbb{Y}=\{\mathscr{L}_{\text{SPE}}(C_{h}^{(g)})-\bar{\mathscr{L}}_{\text{SPE}}:1\leq g\leq G_{\text{REG}}\}^{\top}\in\mathbb{R}^{G_{\text{REG}}}$, where $\bar{\mathscr{L}}_{\text{SPE}}=G_{\text{REG}}^{-1}\sum_{g=1}^{G_{\text{REG}}}\mathscr{L}_{\text{SPE}}(C_{h}^{(g)})$. We define $X^{(g)}=(N^{-4/7}C_{h}^{(g)4}/4,N^{-4/7}/C_{h}^{(g)3})^{\top}\in\mathbb{R}^{2}$. We further define a pseudo-design matrix, $\mathbb{X}\in\mathbb{R}^{G_{\text{REG}}}$, where the $g$th row is $(X^{(g)}-\bar{X})^{\top}$ and $\bar{X}=G_{\text{REG}}^{-1}\sum_{g=1}^{G_{\text{REG}}}X^{(g)}$. Then, we have an appropriate regression relationship as $\mathbb{Y}=\mathbb{X}\mathcal{C}$, where $\mathcal{C}=(C_{1},C_{2})^{\top}\in\mathbb{R}^{2}$. Subsequently, an ordinary least squares estimator is obtained for $\mathcal{C}$ as $\hat{\mathcal{C}}=(\mathbb{X}^{\top}\mathbb{X})^{-1}(\mathbb{X}^{\top}\mathbb{Y})=(\hat{C}_{1},\hat{C}_{2})^{\top}$. This leads to a regression-based estimator for the optimal bandwidth constant $C_{h}^{*}$ as $\hat{C}_{h}^{\text{REG}}=(3\hat{C}_{2}/\hat{C}_{1})^{1/7}$. For convenience, we abbreviate this method as the REG method.

Extensive simulation studies are conducted to validate the finite sample performance of the proposed KEM method. The entire experiment is designed so that the simulation data can mimic real chest CT data as much as possible. Specifically, we first fix a total of $M=3$ mixture components, which represent the background, bone tissue, and lung tissue, respectively. Next, we need to specify the class prior probability function, $\pi_{m}(x)$. To make the experiment as realistic as possible, we use the LIDC-IDRI dataset, which is probably the largest publicly available chest CT database (Setio et al., 2017). It consists of 888 human chest CT files from seven medical institutions worldwide. After downloading the dataset, we found that six files were damaged and could not be read into memory. Thus, we use the remaining 882 CT files for the subsequent study.

Let $\mathbb{Y}=(Y_{ijk})\in\mathbb{R}^{d_{x}\times d_{y}\times d_{z}}$ be an arbitrarily selected chest CT data from the LIDC-IDRI dataset. According to the definition of the database, we should always have $d_{x}=d_{y}=512$ for every CT scan file. However, the number of slices (i.e., $d_{z}$) may vary across different CT scan files, but it should fall within the range of $d_{z}\in[95,764]$. Next, we consider how to specify the values of $\pi_{m}(X_{ijk})$ for every possible voxel position $X_{ijk}=(i/d_{x},j/d_{y},k/d_{z})^{\top}\in\mathcal{V}$, where $\mathcal{V}\subset[0,1]^{3}$ contains all the available voxel positions. Next, we define a class label tensor $\mathbb{Z}=(Z_{ijk})\in\mathbb{R}^{d_{x}\times d_{y}\times d_{z}}$ with $Z_{ijk}\in\{1,2,3\}$. Specifically, we define $Z_{ijk}=2$ if its voxel position, $X_{ijk}$, belongs to the lung-mask data provided by the LIDC-IDRI dataset. Otherwise, we define $Z_{ijk}=3$ if $Y_{ijk}>400$. In this case, the $(i,j,k)$th voxel position should be the bone tissue (Molteni, 2013). Finally, we define $Z_{ijk}=1$ if the previous two criteria are not met.

For the $(i,j,k)$th voxel position $x=(i/d_{x},j/d_{y},k/d_{z})^{\top}=(x_{1},x_{2},x_{3})^{\top}\in\mathcal{V}$, we define $x$ for a local neighborhood as $\mathcal{N}(x)=\big{\{}(x_{1}^{\prime},x_{2}^{\prime},x_{3}^{\prime})^{\top}:|(x_{1}^{\prime}d_{x},x_{2}^{\prime}d_{y},x_{3}^{\prime}d_{z})^{\top}-(x_{1}d_{x},x_{2}d_{y},x_{3}d_{z})^{\top}|_{\max}<3\big{\}}$, where $|x|_{\max}=\mathop{\max}\limits_{1\leq j\leq 3}|x_{j}|$. Then, we define $\pi_{m}(x)=\{\sum_{x^{\prime}\in\mathcal{N}(x)}I(Z_{i^{\prime}j^{\prime}k^{\prime}}=m)\}/|\mathcal{N}(x)|$, where $x^{\prime}=(i^{\prime}/d_{x},j^{\prime}/d_{y},k^{\prime}/d_{z})^{\top}=(x_{1}^{\prime},x_{2}^{\prime},x_{3}^{\prime})^{\top}\in\mathcal{V}$. Next, we redefine $\pi_{m}(x)=\{\pi_{m}(x)+0.6\}/2.8$ for $1\leq m\leq M$. By doing so, we can guarantee that $0<\pi_{m}(x)<1$ for every $1\leq m\leq M$ and $x\in\mathcal{V}$. Next, we define $\mu_{m}=\{\mathop{\sum}\limits_{x\in\mathcal{V}}I(Z_{ijk}=m)Y_{ijk}\}/\{\mathop{\sum}\limits_{x\in\mathcal{V}}I(Z_{ijk}=m)\}$ and $\sigma_{m}^{2}=\mathop{\sum}\limits_{x\in\mathcal{V}}\{I(Z_{ijk}=m)(Y_{ijk}-\mu_{m})^{2}\}/\{\mathop{\sum}\limits_{x\in\mathcal{V}}I(Z_{ijk}=m)\}$ for $2\leq m\leq M$. The $m=1$ case is set as the background case. To differentiate the background from the other classes, $\mu_{1}$ is set to 1, and $\sigma_{1}$ is set to $\sigma_{3}$. Then, we set $\mu_{m}(x)=\mu_{m}+0.25\times\sin(8\pi x_{1})\times\sin(8\pi x_{2})\times\sin(8\pi x_{3})$ and $\sigma_{m}(x)=\sigma_{m}+\sigma_{m}\times\sin(8\pi x_{1})\times\sin(8\pi x_{2})\times\sin(8\pi x_{3})$. Thus, we allow the mean function $\mu_{m}(x)$ and variance function $\sigma_{m}^{2}(x)$ to vary within a reasonable range. For intuitive understanding, $\pi_{m}(x)$, $\mu_{m}(x)$, and $\sigma_{m}(x)$ with $1\leq m\leq 3$ generated from an arbitrarily selected CT are graphically displayed in Figure 1(a)–(i). Throughout this article, the depth index for the displayed 2D slices is consistently set to 100. For any voxel position $x\in\mathcal{V}$, once $\mu_{m}(x)$ and $\sigma_{m}(x)$ are given, a normal random variable can be generated. These normal random variables are then combined across different components according to a multinomial distribution with probability $\pi_{m}(x)$ for class $m$ with $1\leq m\leq M$. This leads to a final response $Y_{ijk}$.

Although the KEM algorithm developed in equations (5)–(8) is theoretically elegant, its computation is not immediately straightforward. If one implements the algorithm faithfully, as equations (5)–(8) suggest, the associated computational cost would be unbearable. Consider, for example, the computation of the denominator of equation (6), which is $\sum_{i=1}^{N}\mathbb{K}(X_{i}/h-x/h)$. The resulting computational complexity is $O(N)$, where $N=|\mathcal{V}|=d_{x}\times d_{y}\times d_{z}$ is the total number of voxel positions. Note that this procedure should be performed for each voxel position in $\mathcal{V}$. Then, the overall computational cost becomes $O(N^{2})$ order. Consider, for example, a chest CT of size $512\times 512\times 201$; then, we have $N=5.27\times{10}^{7}$ and $N^{2}=2.78\times{10}^{15}$. Thus, the computational cost is prohibitive, and alleviating the computational burden becomes a critical problem.

To address this problem, we develop a novel solution. We specify $\mathcal{K}(t)=\exp(-t^{2}/2)/\sqrt{2\pi}$ as the Gaussian kernel and $\mathbb{K}(x)=\mathcal{K}(x_{1})\mathcal{K}(x_{2})\mathcal{K}(x_{3})$. Accordingly, the weight assigned to its tail decays toward 0 at an extremely fast rate. Therefore, instead of directly using the full kernel weight $\mathbb{K}(x)$, we consider its truncated version as $\mathbb{K}^{*}(x)=\mathbb{K}(x)I(|(x_{1}d_{x},x_{2}d_{y},x_{3}d_{z})^{\top}|_{\max}<s)$, where $s>0$ represents the filter size. We define $\mathcal{N}_{s}(x)=\{x^{\prime}=(x_{1}^{\prime},x_{2}^{\prime},x_{3}^{\prime})^{\top}\in\mathcal{V}:|(x_{1}^{\prime}d_{x},x_{2}^{\prime}d_{y},x_{3}^{\prime}d_{z})^{\top}-(x_{1}d_{x},x_{2}d_{y},x_{3}d_{z})^{\top}|_{\max}<s\}$ as a cubic space locally around $x$ with volume $|\mathcal{N}_{s}(x)|=s^{3}$. Instead of computing $\sum_{i=1}^{N}\mathbb{K}(X_{i}/h-x/h)$, we can compute $\sum_{i=1}^{N}\mathbb{K}^{*}(X_{i}/h-x/h)=\sum_{1\leq i\leq N}^{X_{i}\in\mathcal{N}_{s}(x)}\mathbb{K}(X_{i}/h-x/h)$. Thus, the computational cost for one voxel position can be significantly reduced from $O(N)$ to $O(s^{3})$, which is approximately the size of set $\{X_{i}\in\mathcal{N}_{s}(x):1\leq i\leq N\}$. This is typically a much-reduced number compared with $N$. Furthermore, this operation can be easily implemented on a GPU device using a standard 3D convolution operation (e.g., the Conv3D function in Keras of TensorFlow). This substantially accelerates the computational speed. To this end, the bandwidth must be appropriately specified according to the filter size of the convolution operations (i.e., $s$). Ideally, the filter should not be too large. Otherwise, the computation cost owing to $|\mathcal{N}_{s}(x)|$ could be extremely high. However, the filter size, $s$, cannot be too small compared to bandwidth $h$. Otherwise, the probability mass of $\mathbb{K}(x)$, which falls into $\mathcal{N}_{s}(x)$, could be too small.

For the REG method, we consider the filter size as $s^{(g)}=2g+1$, where $1\leq g\leq G_{\text{REG}}$ and $G_{\text{REG}}=5$. Accordingly, we set the bandwidth to $h^{(g)}=s^{(g)}/512$. Thus, only odd numbers are considered, so the filter used in the 3D convolution operations can be symmetric about its central voxel position. We wish the filter to be as large as possible. However, as previously discussed, an unnecessarily large filter size would result in prohibitive computational costs. Therefore, we set the largest filter size as 11. For the CV method, following the REG method, we generate $(s^{(g)},h^{(g)})$ for $1\leq g\leq G_{\text{REG}}$. Subsequently, each $h^{(g)}$ is multiplied by 0.3, 0.4, 0.5, 0.6, or 0.7. The bandwidth constant of the two methods is set as $C_{h}^{(g)}=h^{(g)}\times N^{1/7}$, where $N$ represents the sample size.

With the pre-specified filter sizes, both the CV and REG methods introduced in Sections 2.4 and 2.5 can be used to select the optimal bandwidth. Recall that 80% of the voxels are used for training and the remaining 20% are used for testing. Truncated Gaussian kernel and GPU-based convolutions are used to speed up the practical computation. Both the CV and REG methods are used to select the optimal bandwidth constant for each of the 100 randomly selected patients in the LIDC-IDRI dataset. The optimal bandwidth constants selected by the two methods are then boxplotted in Figure 2(a). We find that the REG method tends to select a larger $C_{h}$ than the CV method, leading to relatively larger filter sizes for convolution. The time cost of the two methods for selecting the optimal bandwidth for each patient is boxplotted in Figure 2(b). Evidently, the REG method has a much lower computational cost than the CV method.

To show the superiority of the proposed KEM method, we also compared it with two traditional methods. The two traditional methods in concern are, respectively, $k$-means (MacQueen et al., 1967) and GMM (McLachlan and Peel, 2000). We then compare the aforementioned methods in terms of testing accuracy and estimation accuracy. Testing accuracy reflects the ability to identify the class labels on the testing set, while estimation accuracy reflects the ability to recover the true parameter values. When comparing the testing accuracies, we train the model on the training set and evaluate the testing accuracy on the testing set. See Figure 2(c) for the results of the testing accuracy. The medians of the testing accuracy results for the KEM method with the bandwidth selected by the CV method (KEM-CV), the KEM method with the bandwidth selected by the REG method (KEM-REG), the $k$-means method, and the GMM method are, respectively, 0.9883, 0.9950, 0.9721, and 0.9719. We can see that both the traditional methods suffer from worse testing accuracy as compared with the proposed KEM methods. When comparing the estimation accuracies, all voxel positions of the CT data are used. We can compute $\hat{\pi}_{m}(x)$, $\hat{\mu}_{m}(x)$, and $\hat{\sigma}_{m}(x)$ using the KEM-CV, KEM-REG, $k$-means and GMM methods for every $1\leq m\leq M$ and every $x\in\mathcal{V}$. Because this is a simulation study with the true parameter values specified in Section 3.1, which are known to us, we can evaluate the accuracy of the resulting estimates using the following RMSE criteria: