# Fast Bootstrapping Nonparametric Maximum Likelihood for Latent Mixture
Models

Shijie Wang, Minsuk Shin, and Ray Bai Manuscript received January 7, 2024;
revised February 28, 2024. The work of S. Wang and R. Bai was supported by the
U.S. National Science Foundation under Grant DMS-2015528. The authors S. Wang
and R. Bai are with the Department of Statistics, University of South
Carolina, SC 29208 USA (e-mail<EMAIL_ADDRESS>rbai@mailbox.sc.edu). M.
Shin is with Gauss Labs, CA 94301 USA (e-mail: minsuk000@gmail.com). This
paper supersedes a previously circulated technical report by S. Wang and M.
Shin (https://arxiv.org/pdf/2006.00767v2.pdf).

###### Abstract

Estimating the mixing density of a latent mixture model is an important task
in signal processing. Nonparametric maximum likelihood estimation is one
popular approach to this problem. If the latent variable distribution is
assumed to be continuous, then bootstrapping can be used to approximate it.
However, traditional bootstrapping requires repeated evaluations on resampled
data and is not scalable. In this letter, we construct a generative process to
rapidly produce nonparametric maximum likelihood bootstrap estimates. Our
method requires only a single evaluation of a novel two-stage optimization
algorithm. Simulations and real data analyses demonstrate that our procedure
accurately estimates the mixing density with little computational cost even
when there are a hundred thousand observations.

###### Index Terms:

bootstrap/resampling, deep neural network, generative process, mixing density
estimation, nonparametric maximum likelihood estimation, two-stage algorithm

##  I Introduction

Nonparametric maximum likelihood estimation (_NPMLE_) is a popular methodology
in signal processing applications such as pattern recognition [1, 2], channel
estimation [3], signal recovery [4], and positron emission tomography [5, 6].
NPMLE also has applications in empirical Bayes and regression modeling [7, 8,
9]. Suppose that we observe $\bm{y}=(y_{1},\ldots,y_{n})$ where

$y_{i}\mid\theta_{i}\sim
f(y_{i}\mid\theta_{i}),~{}~{}\theta_{i}\sim\pi(\theta),~{}~{}i=1,\ldots,n,$
(1)

and the density $f$ is known but $\pi$ is unknown. NPMLE aims to estimate the
mixing density (or the prior) $\pi$.

The NPMLE estimator $\widehat{\pi}$ in the latent mixture model (1) [10, 11]
was introduced by [12] and solves the optimization,

$\displaystyle\widehat{\pi}$
$\displaystyle=\underset{\pi\in\mathrm{\Pi}}{\mathrm{argmax}}\sum_{i=1}^{n}\mathrm{log}[\mathbb{E}_{\pi}\\{f(y_{i}\mid\theta_{i})\\}]$
$\displaystyle=\underset{\pi\in\mathrm{\Pi}}{\mathrm{argmax}}\sum_{i=1}^{n}\mathrm{log}\left\\{\int_{\Theta}f(y_{i}\mid\theta_{i})d\pi(\theta_{i})\right\\},$
(2)

where $\mathrm{\Pi}$ denotes the family of all probability distributions on
the parameter space $\Theta$. Under mild regularity conditions, [13] showed
that the solution $\widehat{\pi}$ of (I) exists and is unique, and even if the
true $\pi$ is continuous, $\widehat{\pi}$ is almost surely discrete with a
support of at most $n$ points. Exploiting this property, there exist many
algorithms to solve (I), including the EM algorithm [14, 15] and convex
optimizers [16, 17].

Nevertheless, a discrete NPMLE estimator $\widehat{\pi}$ is unsatisfactory
when the latent distribution is reasonably assumed to be continuous [5, 18].
As a result, smoothed variants of NPMLE have been proposed [5, 7, 2, 19].
However, these methods typically require careful tuning of smoothing
parameters such as bandwidth [19, 2], roughness penalty term [5, 7], and/or
spline degrees of freedom [2, 7].

As an alternative to smoothing, bootstrapping has been shown to be an
effective way to simulate from both prior and posterior continuous densities
[20, 21, 22, 23]. However, bootstrapping is seldomly employed for NPMLE. This
may be because the bootstrap traditionally involves resampling the data with
replacement and repeatedly optimizing weighted objective functions. In the
context of NPMLE, the standard bootstrap procedure requires one to
repetitively optimize

$\widehat{\pi}^{(b)}=\operatornamewithlimits{argmax}_{\pi}\sum_{i=1}^{n}w_{i}^{(b)}\mathrm{log}[\mathbb{E}_{\pi}\\{f(y_{i}\mid\theta_{i})\\}],~{}~{}\bm{w}^{(b)}\overset{\text{iid}}{\sim}\mathbb{P}_{\bm{w}},$
(3)

for some probability measure $\mathbb{P}_{\bm{w}}$. For example, if
$\bm{w}\sim\text{Multinomial}(n,\mathbbm{1}_{n}/n)$, where $\mathbbm{1}_{n}$
denotes an $n$-dimensional one vector, we have the nonparametric bootstrap
[24]. If $\bm{w}\sim n\times\text{Dirichlet}(n,\mathbbm{1}_{n})$, we have the
weighted likelihood bootstrap [22]. Repeatedly solving (3) can be time-
consuming.

To resolve this computational bottleneck and broaden the applicability of
bootstrapping, we build upon the generative bootstrap sampling (GBS) approach
of [25]. We introduce a new generative framework called Generative
Bootstrapping for NPMLE (_GB-NPMLE_). Our contributions are as follows:

  1. 1.

We propose GB-NPMLE for estimating a continuous prior density via
bootstrapping. Different from GBS which only requires weights $\bm{w}$ as
inputs [25], GB-NPMLE needs _additional_ noise inputs $\bm{z}$ to capture the
latent representation of $\pi$ in (1).

  2. 2.

To optimize GB-NPMLE, we introduce a novel two-stage algorithm. In contrast to
traditional bootstrapping, GB-NPMLE requires only a _single_ evaluation of
this two-stage algorithm instead of repetitive evaluations of (3).

##  II Proposed Method

###  II-A GB-NPMLE

Recent developments in generative learning such as variational autoencoders
(VAEs) [26] and generative adversarial networks (GANs) [27, 28, 29] provide
another perspective for NPMLE. Instead of estimating $\pi$, we can construct a
generator $T$ that matches to the target density $\pi$, i.e.
$T(\bm{z})\sim\pi$ where $\bm{z}$ $\in\mathbb{R}^{q}$ follows a known
distribution such as a Gaussian or a uniform. We then reformulate NPMLE
optimization (I) via a generative framework,

$\widehat{T}=\operatornamewithlimits{argmax}_{T\in\mathcal{T}}\sum_{i=1}^{n}\mathrm{log}[\mathbb{E}_{\bm{z}}\\{f(y_{i}\mid
T({\bm{z}})\\}],$ (4)

where ${\bm{z}}\in\mathbb{R}^{q}$ follows (for example) a standard uniform.
[30] considered (4) where the generator $T$ was constructed by a VAE [26] with
a Langevin bias-correction algorithm [31]. Unfortuantely, both the classical
NPMLE $\widehat{\pi}$ in (I) and the generative NPMLE $\widehat{T}$ in (4)
result in discrete solutions [13]. This discreteness is less attractive when
the true latent $\pi$ is continuous, as in many signal processing applications
[5, 6].

To efficiently obtain a bootstrapped NPMLE distribution as a smooth estimator
for $\pi$ in (1), we introduce GB-NPMLE. Our framework (approximately)
generates NPMLE bootstrap estimates via a generator function. The generator
$G:=G(\bm{w},\bm{z})$ takes _both_ noise $\bm{z}\in$ $\mathbb{R}^{q}$ _and_
bootstrap weights $\bm{w}\in$ $\mathbb{R}^{n}$ as inputs. We construct $G$
using a feedforward neural network (FNN) [32, 33], because of the FNN’s
universal approximation properties for a large class of functions [32, 33].

GB-NPMLE incorporates the GBS approach of [25]. However, GBS is only geared
towards uncertainty quantification of a single _point estimate_. Thus, GBS
takes _only_ weights $\bm{w}$ as inputs and _cannot_ be directly used to learn
an entire _probability density_ $\pi$ as in (3). In order to learn an unknown
density $\pi$, GB-NPMLE _also_ requires noise inputs $\bm{z}$, as in (4). In
summary, GB-NPMLE approximates the bootstrap distribution for NPMLE with the
objective function,

$\widehat{G}=\operatornamewithlimits{argmax}_{G}\
\mathbb{E}_{\bm{w}}\Big{[}\sum_{i=1}^{n}w_{i}\
\mathrm{log}[\mathbb{E}_{\bm{z}}\\{f(y_{i}\mid G({\bm{w},\bm{z}})\\}]\Big{]},$
(5)

where $G(\bm{w},\bm{z}):$ $\mathbb{R}^{n+q}\mapsto\mathbb{R}$. Once we have
optimized $G$ (i.e. the weights and biases of the FNN), it is then effortless
to generate novel NPMLE bootstrap estimates as
$\widehat{G}(\bm{w}_{\text{new}},\bm{z}_{\text{new}})$.

To optimize (5) in practice, the expectations $\mathbb{E}_{\bm{z}}$ and
$\mathbb{E}_{\bm{w}}$ are approximated by Monte Carlo averaging, similarly as
in VAEs [26] and GANs [27]. Specifically, we independently sample
$z\sim\text{Unif}(0,1)$ and $\bm{w}$ $\sim
n\times\text{Dirichlet}(n,\mathbbm{1}_{n})$ a sufficiently large number of
times and then approximate the expectations in (5) with sample averages.

Directly optimizing $G$ in (5) using stochastic gradient descent (SGD) [34] is
a nontrivial task. In practice, the Monte Carlo approximation for
$\mathbb{E}_{\bm{z}}$ introduces high variance to the approximate gradient for
SGD because the log functions in (5) greatly shatter the linearity of
$\mathbb{E}_{\bm{z}}$ [35]. Consequently, this can lead to slow convergence
when optimizing the GB-NPMLE objective function (5) [35]. We address these
estimation difficulties with a novel two-stage algorithm.

###  II-B GB-NPMLE Two-Stage Algorithm

To efficiently optimize $G$ in (5), we generalize the generator output
$\theta=G(\bm{w},\bm{z})$ to be a multi-dimensional
$\bm{\theta}\in\mathbb{R}^{l}$ where $l>1$ is the number of realizations of
$\theta$. The resulting sequence $\bm{\theta}=\\{\theta^{(j)}\\}_{j=1}^{l}$
serves as the bootstrap sample _candidates_. As we explain shortly, multiple
realizations of $\theta$ help to stabilize the gradient approximation in SGD
for the Monte Carlo estimate of $\mathbb{E}_{\bm{z}}$ in (5). Unfortunately,
this generalization also leads to correlation across the bootstrap samples.

One way to overcome the correlation issue is to take a random draw of one
$\theta^{(j)}$ from $\bm{\theta}$ according to mixing probabilities
$\bm{\tau}=(\tau_{1},\ldots,\tau_{l})$ for the entries of $\bm{\theta}$.
Repeating this procedure $B$ times produces $B$ independent bootstrap samples.
However, in addition to generator $G$, we now _also_ need to estimate
$\bm{\tau}$. Accordingly, we propose a two-stage algorithm where in Stage I,
we train $G$, and in Stage II, we estimate $\bm{\tau}$.

_GB-NPMLE Algorithm Stage I:_ Fix $\bm{\tau}=(\tau_{1},...,\tau_{l})$ where
$\tau_{1}=...=\tau_{l}=1/l$, and solve the modified GB-NPMLE objective
function,

$\widehat{G}=\operatornamewithlimits{argmax}_{G}\mathbb{E}_{\bm{w}}\Big{[}\sum_{i=1}^{n}w_{i}\
\mathrm{log}[\mathbb{E}_{\bm{z},\gamma}\\{f(y_{i}\mid\bm{e}_{\gamma}^{\top}G({\bm{w},\bm{z}})\\}]\Big{]},$
(6)

where $\gamma\sim$ Multinomial$(1,\bm{\tau})$, $\bm{e}_{\gamma}$ is the
$\gamma$-th unit vector, and
$\bm{e}_{\gamma}^{\top}G({\bm{w},\bm{z}})\in\mathbb{R}$ is the $\gamma$-th
entry of generator output $\bm{\theta}$.

To see how the modified GB-NPMLE objective (6) helps to relieve high Monte
Carlo variance, let
$\mathbb{E}_{\bm{z},\gamma}\\{f(y_{i}|\bm{e}_{\gamma}^{\top}G({\bm{w},\bm{z}})\\}$
in (6) be denoted by
$\mathbb{E}_{\bm{z}}\\{\mathbb{E}_{\gamma}\\{f(\bm{\theta}[\bm{w},\bm{z},\gamma])\\}\\}$,
and similarly, let $\mathbb{E}_{\bm{z}}\\{f(y_{i}|G(\bm{w},\bm{z}))\\}$ be
$\mathbb{E}_{\bm{z}}\\{f(\bm{\theta}[\bm{w},\bm{z}])\\}$ in (5). By the law of
total variance,
$\text{Var}\\{f(\bm{\theta}[\bm{w},\bm{z}])\\}\geq\text{Var}\\{\mathbb{E}_{\gamma}f(\bm{\theta}[\bm{w},\bm{z},\gamma])\\}$.
Hence, the generalization of the generator output dimension to $l>1$ reduces
Monte Carlo variance.

_GB-NPMLE Algorithm Stage II:_ Fix the well-trained generator $\widehat{G}$
from (6) and estimate $\bm{\tau}$ via the Monte Carlo EM (MCEM) algorithm
[36]. In particular, the MCEM updates for each $k$th entry of $\bm{\tau}$ have
the closed form,

$\tau_{k}^{(t+1)}=\frac{1}{n}\sum_{i=1}^{n}\frac{\tau_{k}^{(t)}\cdot\mathbb{E}_{\bm{w},\bm{z}}\\{f(y_{i}\mid\bm{e}_{k}^{\top}\widehat{G}({\bm{w}},\bm{z})\\}}{\sum_{k=1}^{l}\tau_{k}^{(t)}\cdot\mathbb{E}_{\bm{w},\bm{z}}\\{f(y_{i}\mid\bm{e}_{k}^{\top}\widehat{G}({\bm{w}},\bm{z})\\}}.$
(7)

Note that we _already_ obtained an efficient estimator of $\bm{\theta}$ from
Stage I, and hence, the MCEM algorithm tends to have very fast convergence. In
the literature, slow convergence of MCEM mostly lies in obtaining a sequence
of estimators $\bm{\theta}$ [30] . We circumvent this issue by _fixing_ the
generator from Stage I.

Once we have estimated both $G$ (in Stage I) and $\bm{\tau}$ (in Stage II), we
can easily generate a new bootstrap estimate by sampling
$\bm{w}_{\text{new}}$, $\bm{z}_{\text{new}}$, $\gamma_{\text{new}}$, and then
taking
$\theta_{\text{new}}=\bm{e}_{\gamma_{\text{new}}}^{\top}\widehat{G}(\bm{w}_{\text{new}},\bm{z}_{\text{new}})$.
Repeating this procedure $B$ times results in $B$ independent bootstrap
estimates. Thus, in contrast to bootstrapped NPMLE which requires $B$ total
repetitive evaluations of (3), GB-NPMLE requires only a single evaluation of
the two-stage algorithm. The entire two-stage GB-NPMLE procedure is given in
Algorithm 1. Empirical evidence (reported in Appendix -C) demonstrates
sufficient convergence of our algorithm.

Figure 1: Results from one replication of Simulations (i)-(iii). In addition
to GB-NPMLE (solid red), bootstrapped NPMLE (dashed blue), and smoothed NPMLE
(solid green), we also plot the true density (solid black) and the classical
discrete NPMLE (dotted purple). Algorithm 1 GB-NPMLE Two-stage Algorithm

1:  _Set_ $l=100$, $tol=0.001$, $T=2000$, $B=1000$.

2:  _Initialize_ neural network parameters $\phi$ and
$\bm{\tau}=(1/l,\ldots,1/l)^{\top}$

3: procedure: Stage I

4: for $t$ in $1,\ldots,\text{Epoch }T$ do

5: Sample $\bm{w}\sim n\times\text{Dirichlet}(n,\mathbbm{1}_{n})$

6: Sample $z\sim\text{Unif}$(0,1)

7: Sample $\gamma\sim$ Multinomial$(1,\bm{\tau})$

8: Update $G_{\phi}$ by SGD using modified GB-NPMLE objective (6)

9: end for

10: return $G_{\phi}$

11: end procedure

12: procedure: Stage II

13:  _Initialize_ $\bm{\tau}^{new}=(1/l,\ldots,1/l)$ and
$\bm{\tau}^{old}=(0,\ldots,0)$

14: while $\min_{k}(|\tau_{k}^{old}-\tau_{k}^{new}|)\geq tol$ do

15: set $\bm{\tau}^{old}$ = $\bm{\tau}^{new}$

16: for $k$ in $1,\ldots,l$ do

17: Sample $\bm{w}\sim n\times\text{Dirichlet}(n,\mathbbm{1}_{n})$

18: Sample $z\sim\text{Unif}$(0,1)

19: Update $\tau_{k}^{new}$ as in (7)

20: end for

21: end while

22: return $\bm{\tau}^{new}$

23: end procedure

24: procedure: Generate $B$ bootstrap estimates

25: for $b$ in $1,\ldots,B$ do

26: Sample $\bm{w}^{(b)}\sim n\times\text{Dirichlet}(n,\mathbbm{1}_{n})$

27: Sample $z^{(b)}\sim\text{Unif}$(0,1)

28: Sample $\gamma^{(b)}\sim$ Multinomial$(1,\bm{\tau}^{new})$

29: Generate
$\theta^{(b)}=e_{\gamma^{(b)}}^{\top}\widehat{G}(\bm{w}^{(b)},z^{(b)})$

30: end for

31: return $\theta^{(1)},\ldots,\theta^{(B)}$

32: end procedure

##  III Experiments and Results

###  III-A Simulation Studies

We investigate the performance of GB-NPMLE optimized by the two-stage
algorithm from Section II-B. The generator $G$ is an FNN with two hidden
layers and 500 neurons per layer. Sensitivity analysis to this choice of
architecture is presented in Appendix -A. The expectations
$\mathbb{E}_{\bm{w}}$, $\mathbb{E}_{\bm{z}}$ and $\mathbb{E}_{\gamma}$ in the
modified GB-NPMLE objective (6) are approximated by Monte Carlo averages of
$100$ independent samples of $\bm{w}\sim
n\times\text{Dirichlet}(n,\mathbbm{1}_{n})$, $z\sim\text{Unif}$(0,1) and
$\gamma\sim$ Multinomial$(1,\bm{\tau})$, and the length of $\bm{\tau}$ is set
to be $l=100$. To train $G$, we use SGD with the Adam optimizer [37]. All
experiments are conducted on a single NVIDIA GeForce RTX 2080 Ti graphics
processing unit (GPU) with 11GB RAM.

We consider three simulation settings for the prior $\pi$ in the mixture model
(1): (i) bimodal, (ii) unimodal and bounded on $(0,1)$, and (iii) skewed
right. Additional simulations are conducted in Appendix -B. Here, we consider:

  1. (i)

Gaussian mixture model (GMM): $y\mid\theta\sim\mathcal{N}(\theta,1)$ and
$\theta=0.5\mathcal{N}(-3,2)+0.5\mathcal{N}(3,1)$;

  2. (ii)

Gamma mixture model (GaMM): $y\mid\theta\sim\text{Gamma}(10,\theta)$ and
$\theta\sim\text{Beta}(10,5)$;

  3. (iii)

Poisson mixture model (PMM): $y\mid\theta\sim\text{Poisson}(\theta)$ and
$\theta\sim\text{Gamma}(3,1)$.

We compare GB-NPMLE to bootstrapped NPMLE (3) and a smoothed version of the
discrete NPMLE (I). We use $B=10{,}000$ bootstrap estimates to estimate the
prior for GB-NPMLE and bootstrapped NPMLE. For bootstrapped NPMLE, we use the
nonparametric bootstrap, i.e. we sample $n$ instances of the data with
replacement and then optimize (I) a total of $B$ times using the R package
REBayes [38]. For the smoothed NPMLE, we use kernel smoothing of
$\widehat{\pi}$ in (I) with the KWsmooth function in REBayes [38], where we
select the optimal bandwidth from an equispaced grid of length $25$ from $0.1$
to $10$ using leave-one-out cross-validation (LOOCV).

For performance comparison of the three NPMLE methods, we generate artificial
datasets of size $n=1000$. We look at: a) Wasserstein-1 distance between $\pi$
and $\widehat{\pi}$, defined as
$W_{1}(\pi,\widehat{\pi})=\int|\mathbb{F}_{\pi}(x)-\mathbb{F}_{\widehat{\pi}}(x)|dx,$
where $\mathbb{F}$ is the cumulative density function (CDF) of $\pi$, and b)
Integrated Squared Error (ISE), defined as
$\int(\widehat{\pi}(x)-\pi(x))^{2}dx$. Lower $W_{1}(\pi,\widehat{\pi})$ and
ISE indicate better performance. Our results for simulations (i)-(iii)
averaged across 20 replications are shown in Table I.

TABLE I: COMPARISONS OF PERFORMANCE OF DIFFERENT NPMLE METHODS ON SIMULATED
DATASETS

| GB-NPMLE | Bootstrapped NPMLE | Smoothed NPMLE
---|---|---|---
Model | $W_{1}(\pi,\widehat{\pi})$ | ISE | $W_{1}(\pi,\widehat{\pi})$ | ISE | $W_{1}(\pi,\widehat{\pi})$ | ISE
GMM | 0.348 | 0.008 | 0.298 | 0.008 | 0.936 | 0.034
GaMM | 0.032 | 0.263 | 0.037 | 0.497 | 0.223 | 1.731
PMM | 0.400 | 0.045 | 0.389 | 0.036 | 1.763 | 0.059

From Table I, we first observe that GB-NPMLE has a very similar performance to
bootstrapped NPMLE. On the other hand, smoothed NPMLE (with bandwidth
parameter chosen by LOOCV) performs much worse than the bootstrapping
approaches. This demonstrates the difficulty of optimally tuning smoothing
parameters for smoothed NPMLE.

The results from one replication of simulations (i)-(iii) are plotted in Fig.
1, which also plots the true latent density $\pi$ (solid black) and the
discrete NPMLE solution (I) (dotted purple). Fig. 1 shows that the classical
(discrete) NPMLE (I) is inadequate for capturing $\pi$ when $\pi$ is truly
continuous. For example, in the GaMM model, the classical NPMLE does not
clearly indicate that $\pi$ is unimodal. On the other hand, GB-NPMLE (solid
red) and bootstrapped NPMLE (solid blue) are quite close to the true $\pi$,
capturing all inherent aspects of the true $\pi$ like bimodality, boundedness
and unimodality, and skewness respectively in the GMM, GaMM, and PMM models.
The smoothed NPMLE appears to be oversmoothed and is unable to capture
bimodality of $\pi$ in the GMM model or the unimodality and boundedness of
$\pi$ in the GaMM model. In the PMM model, the sharp peak of the skewed
density $\pi$ is also not adequately captured by smoothed NPMLE.

Figure 2: Mean computation time in log(seconds) for the three simulations vs.
sample size $n\in\\{1000,10{,}000,100{,}000\\}$. For $n=100{,}000$, GB-NPMLE
takes five minutes, whereas bootstrapped NPMLE takes almost two days.

GB-NPMLE and bootstrapped NPMLE both give similar performance, but GB-NPMLE is
much faster because it only requires a single evaluation of Algorithm 1. We
compare the average computational time of GB-NPMLE and bootstrapped NPMLE for
sample sizes $n\in\\{1000,10{,}000,100{,}000\\}$ across simulations (i)-(iii).
Our results are plotted (on the log scale) in Fig. 2. Fig. 2 shows that GB-
NPMLE is much more scalable. In particular, GB-NPMLE takes only about five
minutes to complete when $n=100{,}000$; on the other hand, bootstrapped NPMLE
requires almost two days to finish.

###  III-B Evaluation on Real Datasets

We evaluate GB-NPMLE on three count datasets where the counts are assumed to
follow a Poisson mixture model,
$y_{i}\mid\theta_{i}\sim\text{Poisson}(\theta_{i})$ and
$\theta_{i}\sim\pi,i=1,\ldots,n$. The _Norberg_ dataset consists of $1125$
group life insurance statistics in Norway [39], the _Thailland_ dataset
consists of 602 preschool children’s health condition, and the _Mortality_
dataset contains the number of deaths for women aged greater than eighty [40].

Figure 3: Results for real datasets for GB-NPMLE (dashed red, top panel) and
bootstrapped NPMLE (dotted blue, bottom panel). The pink vertical lines in
both plots represent the discrete NPMLE solution.

Fig. 3 compares GB-NPMLE (top panel) with bootstrapped NPMLE (bottom panel).
We see that GB-NPMLE approximates bootstrapped NPMLE very well. Furthermore,
both methods identify the same number of local modes (three in _Norberg_ ,
four in _Thailand_ , and two in _Mortality_).

We also examine out-of-sample prediction using the log predictive score (LPS)
$\mathbb{E}_{\theta}[-\log
f(\bm{y}_{\text{new}};\theta)\mid\bm{y}_{\text{obs}}]$. To approximate LPS, we
use $K$-fold cross validation (CV),

$\text{LPS}=K^{-1}\sum_{k=1}^{K}\sum_{i\in
I_{k}}-\log\bigg{[}B^{-1}\sum_{b=1}^{B}f(y_{i};\widehat{\theta}^{b}_{(-k)})\bigg{]},$

where $\widehat{\theta}^{b}_{(-k)}\sim\widehat{\pi}_{(-k)}$ and
$\widehat{\pi}_{(-k)}$ is the bootstrap distribution obtained by excluding the
$k$th fold $I_{k}$. We set $K=10$ and $B=500$. Table II shows that the GB-
NPMLE predictions are quite close to their standard bootstrap counterparts..
Table II also shows that GB-NPMLE is much more scalable, with only a modest
increase in computation time for 10-fold CV when $n$ increases from 72 to
1096. When $n=1096$, GB-NPMLE is about 34 times faster than bootstrapped
NPMLE.

TABLE II: COMPARISONS OF PERFORMANCE OF GB-NPMLE AND BOOTSTRAPPED NPMLE ON
REAL COUNT DATASETS

|  | GB-NPMLE | Bootstrapped NPMLE
---|---|---|---
Dataset | $n$ | LPS | Time (sec) | LPS | Time (sec)
Norberg | 72 | 20.35 | 126.1 | 20.83 | 732.8
Thailland | 602 | 156.26 | 140.9 | 156.21 | 3439.3
Mortality | 1096 | 199.30 | 153.3 | 199.31 | 5279.2

##  IV Conclusion

In this letter, we introduced GB-NPMLE, a generative framework for estimating
a mixing density through bootstrapping. A novel two-stage algorithm was
proposed to rapidly and accurately obtain bootstrap estimates. Bootstrapping
has rarely been employed for NPMLE in the past, and our work paves the way for
bootstrapping to serve as an attractive alternative to smoothing NPMLE.
Simulations and real data analyses demonstrated the effectiveness and
scalability of our approach. Codes to implement GB-NPMLE and the real datasets
analyzed in this article are available on GitHub at
https://github.com/shijiew97/GBnpmle.

Our methodology is only presented for estimating a univariate latent density.
However, it will be useful to extend GB-NPMLE to _multivariate_ density
estimation. If covariates are also observed, then GB-NPMLE can be extended to
mixture of regression models [8]. Finally, theoretical analysis of the two-
stage algorithm (e.g. the convergence rate) is needed.

## References

  * [1] H. T. Huynh and L. Nguyen, “Nonparametric maximum likelihood estimation using neural networks,” _Pattern Recognition Letters_ , vol. 138, pp. 580–586, 2020.
  * [2] Y. Li and Z. Ye, “Boosting in univariate nonparametric maximum likelihood estimation,” _IEEE Signal Processing Letters_ , vol. 28, pp. 623–627, 2021\.
  * [3] V. Bhatia and B. Mulgrew, “Non-parametric likelihood based channel estimator for Gaussian mixture noise,” _Signal Processing_ , vol. 87, no. 11, pp. 2569–2586, 2007.
  * [4] L. Feng, R. Ma, and L. H. Dicker, “Nonparametric maximum likelihood approximate message passing,” in _2017 51st Annual Conference on Information Sciences and Systems (CISS)_. IEEE, 2017, pp. 1–6.
  * [5] B. W. Silverman, M. C. Jones, J. D. Wilson, and D. W. Nychka, “A smoothed EM approach to indirect estimation problems, with particular, reference to stereology and emission tomography,” _Journal of the Royal Statistical Society. Series B (Methodological)_ , vol. 52, no. 2, pp. 271–324, 1990.
  * [6] Y. Vardi, L. A. Shepp, and L. Kaufman, “A statistical model for positron emission tomography,” _Journal of the American Statistical Association_ , vol. 80, no. 389, pp. 8–20, 1985.
  * [7] B. Efron, “Empirical Bayes deconvolution estimates,” _Biometrika_ , vol. 103, no. 1, pp. 1–20, 2016.
  * [8] H. Jiang and A. Guntuboyina, “A nonparametric maximum likelihood approach to mixture of regression,” _arXiv preprint arXiv:2108.09816_ , 2021.
  * [9] Z. Fan, L. Guan, Y. Shen, and Y. Wu, “Gradient flows for empirical Bayes in high-dimensional linear models,” _arXiv preprint arXiv:2312.12708_ , 2023\.
  * [10] S. Frühwirth-Schnatter, _Finite Mixture and Markov Switching Models_. Springer, 2006.
  * [11] G. J. McLachlan, S. X. Lee, and S. I. Rathnayake, “Finite mixture models,” _Annual Review of Statistics and Its Application_ , vol. 6, pp. 355–378, 2019\.
  * [12] J. Kiefer and J. Wolfowitz, “Consistency of the maximum likelihood estimator in the presence of infinitely many incidental parameters,” _The Annals of Mathematical Statistics_ , vol. 27, no. 4, pp. 887–906, 1956.
  * [13] B. G. Lindsay, “Mixture models: Theory, geometry and applications,” in _NSF-CBMS Regional Conference Series in Probability and Statistics_ , vol. 5, 1995, pp. 1–163.
  * [14] N. Laird, “Nonparametric maximum likelihood estimation of a mixing distribution,” _Journal of the American Statistical Association_ , vol. 73, no. 364, pp. 805–811, 1978.
  * [15] C. H. Zhang, “Compound decision theory and empirical Bayes methods,” _Annals of Statistics_ , vol. 31, no. 2, pp. 379–390, 2003.
  * [16] R. Koenker and I. Mizera, “Convex optimization, shape constraints, compound decisions, and empirical Bayes rules,” _Journal of the American Statistical Association_ , vol. 109, no. 506, pp. 674–685, 2014.
  * [17] L. Feng and L. H. Dicker, “Approximate nonparametric maximum likelihood for mixture models: A convex optimization approach to fitting arbitrary multivariate mixing distributions,” _Computational Statistics & Data Analysis_, vol. 122, pp. 80–91, 2018.
  * [18] L. Liu, M. Levine, and Y. Zhu, “A functional EM algorithm for mixing density estimation via nonparametric penalized likelihood maximization,” _Journal of Computational and Graphical Statistics_ , vol. 18, no. 2, pp. 481–504, 2009.
  * [19] R. Koenker and J. Gu, “Comment: Minimalist g-modeling,” _Statistical Science_ , vol. 34, no. 2, pp. 209–213, 2019.
  * [20] M. Knott and P. Tzamourani, “Bootstrapping the estimated latent distribution of the two-parameter latent trait model,” _British Journal of Mathematical and Statistical Psychology_ , vol. 60, no. 1, pp. 175–191, 2007.
  * [21] D. B. Rubin, “The Bayesian bootstrap,” _The Annals of Statistics_ , vol. 9, no. 1, p. 130434, 1981.
  * [22] M. A. Newton and A. E. Raftery, “Approximate Bayesian inference with the weighted likelihood bootstrap,” _Journal of the Royal Statistical Society: Series B (Methodological)_ , vol. 56, no. 1, pp. 3–26, 1994.
  * [23] M. A. Newton, N. G. Polson, and J. Xu, “Weighted Bayesian bootstrap for scalable posterior distributions,” _Canadian Journal of Statistics_ , vol. 49, no. 2, pp. 421–437, 2021.
  * [24] B. Efron, “Bootstrap methods: Another look at the jackknife,” _The Annals of Statistics_ , vol. 7, no. 1, pp. 1–26, 1979.
  * [25] M. Shin, S. Wang, and J. S. Liu, “Generative multiple-purpose sampler for weighted M-estimation,” _Journal of Computational and Graphical Statistics (to appear)_ , 2023.
  * [26] D. P. Kingma and M. Welling, “Auto-encoding variational bayes,” _arXiv preprint arXiv:1312.6114_ , 2013.
  * [27] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial nets,” _Advances in Neural Information Processing Systems_ , vol. 27, 2014.
  * [28] X. Zhou, Y. Jiao, J. Liu, and J. Huang, “A deep generative approach to conditional sampling,” _Journal of the American Statistical Association_ , vol. 118, no. 543, pp. 1837–1848, 2023.
  * [29] S. Liu, X. Zhou, Y. Jiao, and J. Huang, “Wasserstein generative learning of conditional distribution,” _arXiv preprint arXiv:2112.10039_ , 2021.
  * [30] Y. Qiu and X. Wang, “ALMOND: Adaptive latent modeling and optimization via neural networks and langevin diffusion,” _Journal of the American Statistical Association_ , vol. 116, no. 535, pp. 1224–1236, 2021.
  * [31] G. O. Roberts and R. L. Tweedie, “Exponential convergence of Langevin distributions and their discrete approximations,” _Bernoulli_ , vol. 2, no. 4, pp. 341–363, 1996.
  * [32] K. Hornik, M. Stinchcombe, and H. White, “Multilayer feedforward networks are universal approximators,” _Neural Networks_ , vol. 2, no. 5, pp. 359–366, 1989.
  * [33] K. Hornik, “Approximation capabilities of multilayer feedforward networks,” _Neural networks_ , vol. 4, no. 2, pp. 251–257, 1991.
  * [34] F. Emmert-Streib, Z. Yang, H. Feng, S. Tripathi, and M. Dehmer, “An introductory review of deep learning for prediction models with big data,” _Frontiers in Artificial Intelligence_ , vol. 3, 2020.
  * [35] X. Lian, M. Wang, and J. Liu, “Finite-sum composition optimization via variance reduced gradient descent,” in _Proceedings of the 20th International Conference on Artificial Intelligence and Statistics_ , A. Singh and J. Zhu, Eds., vol. 54. PMLR, 20–22 Apr 2017, pp. 1159–1167.
  * [36] R. A. Levine and G. Casella, “Implementations of the Monte Carlo EM algorithm,” _Journal of Computational and Graphical Statistics_ , vol. 10, no. 3, pp. 422–439, 2001.
  * [37] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [38] R. Koenker and J. Gu, “REBayes: An R package for empirical Bayes mixture methods,” _Journal of Statistical Software_ , vol. 82, no. 8, p. 1–26, 2017\.
  * [39] R. Norberg, “Experience rating in group life insurance,” _Scandinavian Actuarial Journal_ , vol. 1989, no. 4, pp. 194–224, 1989.
  * [40] D. Böhning, _Computer-Assisted Analysis of Mixtures and Applications: Meta-analysis, disease mapping and others_. CRC Press, 1999.

###  -A Sensitivity Analysis for Generator Architecture

We conduct a sensitivity analysis to study whether different choices for the
number of hidden layers $L$ or the number of hidden neurons $h$ in the
generator $G$ affects GB-NPMLE’s performance. We first fix $L=2$ and vary
$h\in\\{50,250,500,750\\}$. Next, we fix $h=500$ and vary $L\in\\{1,2,3,4\\}$.

We repeat all of the simulations described in Section III-A 20 times and
record the average performance metrics (i.e. $W_{1}(\pi,\widehat{\pi})$ and
ISE). Our results are summarized in Table III. We do not observe any
significant differences in the performance of GB-NPMLE for these different
combinations of $(L,h)$. In practice, we recommend a default choice of $L=2$
and $h=500$ for satisfactory performance.

TABLE III: SENSITIVITY ANALYSIS RESULTS

$L=2$1! | $h=50$ | $h=250$ | $h=500$ | $h=750$
---|---|---|---|---
Model | $W_{1}(\pi,\widehat{\pi})$ | ISE | $W_{1}(\pi,\widehat{\pi})$ | ISE | $W_{1}(\pi,\widehat{\pi})$ | ISE | $W_{1}(\pi,\widehat{\pi})$ | ISE
GMM | 0.334 | 0.008 | 0.326 | 0.008 | 0.348 | 0.008 | 0.329 | 0.008
GaMM | 0.032 | 0.223 | 0.032 | 0.282 | 0.032 | 0.263 | 0.032 | 0.245
PMM | 0.413 | 0.033 | 0.380 | 0.028 | 0.400 | 0.045 | 0.386 | 0.030

$h=500$ | $L=1$ | $L=2$ | $L=3$ | $L=4$
---|---|---|---|---
Model | $W_{1}(\pi,\widehat{\pi})$ | ISE | $W_{1}(\pi,\widehat{\pi})$ | ISE | $W_{1}(\pi,\widehat{\pi})$ | ISE | $W_{1}(\pi,\widehat{\pi})$ | ISE
GMM | 0.330 | 0.008 | 0.348 | 0.008 | 0.331 | 0.008 | 0.332 | 0.008
GaMM | 0.032 | 0.292 | 0.032 | 0.263 | 0.032 | 0.304 | 0.032 | 0.270
PMM | 0.400 | 0.029 | 0.400 | 0.045 | 0.398 | 0.035 | 0.393 | 0.036

###  -B More Simulation Studies

We investigate the performance of GB-NPMLE under two additional settings for
the prior $\pi$: (iv) a trimodal density, and (v) a bounded and skewed left
density. Namely, we consider:

  1. (iv)

Gaussian trimodal mixture (GMM-tri): $y\mid\theta\sim\mathcal{N}(\theta,1)$
and $\theta=0.2\mathcal{N}(-4,0.5)+0.6\mathcal{N}(0,1)+0.2\mathcal{N}(4,0.5)$;

  2. (v)

Binomial-beta mixture (BBM): $y\mid\theta\sim\text{Binomial}(10,\theta)$ and
$\theta\sim\text{Beta}(3,2)$.

Figure 4: Results from one replication of Simulations (iv)-(v). In addition to
GB-NPMLE (solid red), bootstrapped NPMLE (dashed blue), and smoothed NPMLE
(solid green), we also plot the true density (solid black) and the classical
discrete NPMLE (dotted purple).

As shown in Fig. 4, GB-NPMLE approximates the bootstrapped NPMLE quite well
for both GMM-tri and BBM. GB-NPMLE and bootstrapped NPMLE also provide density
estimates that are close to the true $\pi$, whereas the smoothed NPMLE with
bandwidth chosen from LOOCV is once again oversmoothed. A performance
comparison of GB-NPMLE, bootstrapped NPMLE, and smoothed NPMLE based on the
average of 20 replications is summarized in Table IV.

TABLE IV: COMPARISONS OF PERFORMANCE OF DIFFERENT NPMLE METHODS

| GB-NPMLE | Bootstrapped NPMLE | Smoothed NPMLE
---|---|---|---
Model | $W_{1}(\pi,\widehat{\pi})$ | ISE | $W_{1}(\pi,\widehat{\pi})$ | ISE | $W_{1}(\pi,\widehat{\pi})$ | ISE
GMM-tri | 0.213 | 0.058 | 0.273 | 0.071 | 0.908 | 0.031
BBM | 0.033 | 0.059 | 0.035 | 0.055 | 0.109 | 0.320

###  -C Convergence Analysis

To empirically check the convergence of the proposed two-stage algorithm, we
plot the GB-NPMLE loss vs. epoch number for training the generator $G$ in
Stage I. For Stage II, we plot the log-likelihood vs. MCEM iteration number
for learning the multinomial weights $\bm{\tau}$. Fig. 5 plots these solution
paths for all 40 replications of the Gaussian and Gamma mixture models (i.e.
Experiments (i) and (ii) in Section III-A, denoted as GMM and GaMM
respectively.

Fig. 5 shows that Stage I converges quickly, typically within 50 epochs.
Meanwhile, Stage II converges even faster – typically within one or two
iterations, resulting in a very flat solution path. The convergence plots for
our other numerical experiments were very similar to those in Fig. 5. Although
empirical evidence suggests that the GB-NPMLE two-stage algorithm converges
quickly, a rigorous theoretical analysis of convergence rate should be done
for future work.

Figure 5: Convergence plots for Stage I and Stage II of the two-stage
algorithm under the GMM and GaMM models (Section III-A).