# A Bayesian Nonparametric Framework for Uncertainty Quantification in
Simulation

Wei Xie111Department of Mechanical and Industrial Engineering, Northeastern
University<EMAIL_ADDRESS>Cheng Li222Department of Statistics and
Applied Probability, National University of Singapore<EMAIL_ADDRESS>Yuefeng Wu333Department of Mathematics and Computer Science, University of
Missouri St. Louis<EMAIL_ADDRESS>and Pu Zhang444Department of Industrial and
Systems Engineering, Rensselaer Polytechnic Institute<EMAIL_ADDRESS>
###### Abstract

When we use simulation to assess the performance of stochastic systems, the
input models used to drive simulation experiments are often estimated from
finite real-world data. There exist both input and simulation estimation
uncertainty in the system performance estimates. Without strong prior
information on the input models and the system mean response surface, in this
paper, we propose a Bayesian nonparametric framework to quantify the impact
from both sources of uncertainty. Specifically, nonparametric input models are
introduced to faithfully capture the important features of the real-world
data, such as heterogeneity, multi-modality and skewness. Bayesian posteriors
of flexible input models characterize the input uncertainty, which
automatically accounts for both model selection and parameter value
uncertainty. Then, the input uncertainty is propagated to outputs by using
direct simulation. Thus, under very general conditions, our framework delivers
an empirical credible interval accounting for both input and simulation
uncertainties. A variance decomposition is further developed to quantify the
relative contributions from both sources of uncertainty. Our approach is
supported by rigorous theoretical and empirical study.

Key words— Bayesian nonparametric model, stochastic simulation, input
modeling, input uncertainty, uncertainty quantification

##  1 Introduction

Simulation is widely used in many applications to assess the performance of
stochastic systems, e.g., manufacturing, supply chain and health care systems.
The input models, defined as the driving stochastic processes in simulation
experiments, are often estimated from finite real-world data. Therefore, there
exist both input estimation error, called the input uncertainty, and
simulation estimation error, called the simulation uncertainty. Ignoring
either source of uncertainty could lead to unfounded confidence in the
simulation assessment of system performance.

Various approaches have been proposed in the literature to quantify the input
and simulation uncertainties; see Barton (2012); Song et al. (2014) and Lam
(2016) for a comprehensive review. Based on methodologies developed for
characterizing the input uncertainty, they can be divided into frequentist and
Bayesian approaches. The frequentist approaches typically study the sampling
distributions of point estimators of underlying input models. Since it could
be hard to get the exact sampling distributions in many situations, the
asymptotic approximation, including the normal approximation and the
bootstrap, is often used to quantify the input uncertainty, which is valid
when the amount of real-world data is large. However, even in the current big
data world, we often face the situations where the amount of real-world data
is limited, especially for high-tech products with short life cycles. For
example, biopharma manufacturing requires 9 to 12 months from raw materials
sourcing to the finished drug products, and it requires another 2 to 3 months
for quality testing. However, the drug substances typically expire after 18 to
36 months; see Otto et al. (2014). Compared to frequentist methods, Bayesian
approaches derive the posterior distributions quantifying the input
uncertainty and they do not need a large-sample asymptotic approximation for
their validation. It is also straightforward for Bayesian approaches to
incorporate the prior information about the underlying input models; see Xie
et al. (2014) for the discussion of frequentist v.s. Bayesian approaches for
input uncertainty.

In this paper, we focus on developing a Bayesian nonparametric framework to
quantify the estimation uncertainty of system mean performance when we do not
have strong prior information on the input models and system mean response
surface. We consider univariate input models, which model independent and
identically distributed (i.i.d.) data by mutually independent input
distributions.

Many existing methods assume specific parametric families for input models
with unknown parameter values estimated from finite real-world data; see the
review in Barton (2012). The input uncertainty can be quantified by the
posteriors of input parameters. However, for most practical problems, we do
not know the underlying input families. In addition, the real-world data tend
to have rich properties, such as skewness and multi-modality, and there is no
single standard parametric model that can capture these properties well. Even
flexible parametric Johnson Transformation System (JTS) can only capture
unimodal and bimodal distributional shapes (Biller and Nelson, 2009; Biller
and Gunes, 2010). If the selected parametric families do not have sufficient
flexibility and cannot represent the underlying input models well, there
always exists the distribution family selection error which does not vanish as
the amount of real-world data becomes large. This inconsistent estimation
could lead to incorrect inference even for the moderate size of real-world
data (Hjort et al., 2011).

One possible remedy for the inconsistency of parametric approaches is to
introduce the family uncertainty, which accounts for the input model selection
error among a pre-specified pool of candidate parametric families. Chick
(2001) proposed Bayesian Model Averaging (BMA) to quantify input uncertainty
from both families and parameter values, where the family uncertainty is
characterized by the posterior probabilities of different candidate parametric
models. Since it is difficult for any standard parametric family to capture
the rich features in the real-world data, it could be challenging to select
the appropriate candidate models for BMA; see Section 14.1 in Bishop (2006).
Furthermore, if the selected families are not mutually exclusive, such as
exponential and Gamma distributions, it can potentially lead to model
identification problems.

Considering the limitation of parametric families, we explore nonparametric
input modeling in the framework that we develop to quantify the overall system
performance estimation uncertainty. Our approach is based on the following
consideration: In many situations, the real-world data represent the
variability caused by various latent sources of uncertainty. For example, in a
raw material (RM) inventory system in the bio-pharmaceutical manufacturing, a
single RM is used to satisfy the demands of various products from many
production lines, called move orders. Since the variability of demands induced
by latent sources of uncertainty could be different, the distribution of move
orders could exhibit heterogeneity. When the modes of demands from different
sources of uncertainty are separate, the distribution of move orders typically
exhibits multi-modality. In addition, contamination and cross-contamination in
the production processes could cause the shutdown of production lines and the
throwaway of batches of products, which can lead to extreme large needs of RM
for reproduction. In such cases, the distribution of move orders can exhibit
right skewness and tail. As a result, the underlying physical input
distributions characterizing the variability induced by various latent sources
of uncertainty could have the rich features, including heterogeneity, multi-
modality, skewness and tails. These properties are also observed in the real-
world data collected from other industries; see for example Wagner et al.
(2009); Ma (2011) and Akcay et al. (2011).

Flexible Bayesian nonparametric input models are presented to efficiently
capture important features in the real-world data. For discrete random
variables with finite support points, the multinomial distribution could be
used as a straightforward nonparametric estimator of the true underlying
distribution. Therefore, in this paper, we focus on the input modeling of
continuous random variables whose support is one of the following: (i)
$\Re^{+}\equiv[0,\infty)$, the nonnegative half real line; (ii)
$\Re\equiv(-\infty,\infty)$, the entire real line; or (iii) a bounded interval
$[a_{1},a_{2}]\in\Re$. Specifically, our input models are based on the
Dirichlet Processes Mixtures (DPM), a popular Bayesian method in both
statistics and machine learning communities. For details about the DPM, we
refer the readers to Ghosh and Ramamoorthi (2003). Motivated by the kernel
density estimation (KDE), the DPM with Gaussian kernel was introduced in the
statistics community; see for example West 1990, Escobar and West 1995, etc.
It is extended to other kernel functions; see for example Hanson (2006);
Kottas (2006); Wu and Ghosal (2008).

From the modeling perspective, DPM has clear advantages over standard
parametric families because the variability across different mixing components
naturally represents various latent sources of uncertainty, which makes it
straightforward to capture the important properties in the real-world data.
Different from parametric approaches, the number of active mixing components
and parameters can automatically adjust to the complexity of real data. Thus,
our empirical study demonstrates that DPM has better and more robust finite
sample performance. From the theoretical perspective, DPM is able to
consistently estimate a wide class of distributions under relatively general
conditions (Ghosal et al. 1999b, Wu and Ghosal 2008, etc.). Compared to BMA,
our approach avoids the difficulty of selecting the “appropriate” candidate
distributions. From the computational perspective, one can develop efficient
posterior samplers for DPM of popular exponential families (see our Section
3.2, Escobar and West 1995, Neal 2000, etc.).

Among frequentist approaches, the empirical distribution is the most commonly
used nonparametric approach in the simulation literature, and the bootstrap is
typically used to quantify the input estimation uncertainty; see for example
Barton and Schruben (1993); Barton (2007). Empirical distribution is simple
and easy to implement. However, DPM has some important advantages compared to
empirical distribution. First, even though the underlying true distribution is
continuous, empirical distribution is always discrete. With limited real-world
data, the empirical distribution could overlook some important properties in
the underlying input models. For example, in the presence of extreme values,
DPM with infinite mixture components can provide a better fit to the tails.
Second, the validity of using the bootstrap to quantify the input uncertainty
relies on asymptotic approximation, and therefore it requires large samples of
real-world data. As we mentioned above, the decision makers often face the
situations where the amount of real-world data is limited. As a Bayesian
approach, DPM can overcome these limitations. Our empirical study demonstrates
that DPM has better finite sample performance compared to frequentist
competitors, especially when the sample size of real-world data is limited.
Our empirical study indicates that it has the potential to detect the latent
sources of uncertainty. The selection of kernel function may impact the
performance of DPM.

Therefore, in this paper, we develop a flexible Bayesian nonparametric
framework to quantify the system mean response estimation uncertainty. We
first introduce nonparametric input models based on DPM with various kernels,
which could capture the important properties in the real-world data. The
samples drawn from posteriors of flexible input models can automatically
quantify both model selection and parameters value uncertainty. Then, the
input uncertainty is propagated to the output through direct simulation which
runs simulations at each posterior sample of input models to estimate the
system mean response. Our framework leads to a sampling procedure that
delivers a percentile empirical credible interval (CrI) quantifying the
overall uncertainty of system performance estimation. Specifically, given the
real-world input data, the estimation uncertainty of system mean response is
quantified by the distribution of a compound random variable. Since the input
uncertainty is characterized by the posterior distribution of input models and
the simulation uncertainty is characterized by the sampling distribution of
mean response estimator, we call the interval induced by the distribution of
compound random variable as empirical credible interval.

In sum, the main contributions of our paper are as follows.

  1. 1.

Considering that the real-world data represent the variability caused by
various latent sources of uncertainty in many situations, DPM provides
sufficient flexibility to capture the important features in the real-world
data. Further, DPM with Gamma, Gaussian and Beta kernels can model input data
with support on the non-negative half real line, the whole real line, and an
interval with finite length. The empirical study demonstrates that our input
models have better performance than existing approaches in both input density
estimation and system mean response assessment.

  2. 2.

Without strong prior information about the underlying true input distributions
and the system response surface, we propose a Bayesian nonparametric framework
which accounts for both input and simulation uncertainties, and delivers a
percentile empirical CrI quantifying the overall estimation uncertainty of
system mean response. Furthermore, a variance decomposition is developed to
quantify the relative contributions from input and simulation uncertainties.

  3. 3.

We provide the theoretical support for our nonparametric framework. The theory
includes the consistency of nonparametric input models and empirical CrI
accounting for both input and simulation uncertainties. As the amount of real-
world data increases, without strong prior information on the distribution
family, the posterior distributions of input models can converge to the
underlying distributions. Given the finite real-world input data, as the
simulation budget increases, this interval converges to the CrI quantifying
the impact of input uncertainty with the true mean response surface known.
Further, as the amount of real-world data and the simulation budget go to
infinity, the CrI converges to the true system performance.

The next section describes the related studies on input modeling and
uncertainty quantification. In Section 3, a Bayesian nonparametric framework
is introduced to quantify the overall uncertainty of the system performance
estimates. We then report results of finite sample behaviors on both input
models and system mean response estimation in Section 4, and we conclude this
paper in Section 5. All proofs, derivations and other supplementary studies
are included in the Appendix.

##  2 Background

Since Barton (2012); Song et al. (2014) and Lam (2016) provided the review on
input uncertainty and uncertainty quantification, in this section, we briefly
discuss existing Bayesian approaches related to our approach. When the
parametric families of input models are known, samples drawn from posteriors
of input parameters can quantify the input uncertainty. Two approaches are
typically used to propagate the input uncertainty to the output: direct
simulation and metamodeling. Zouaoui and Wilson (2003) run simulations at each
sample of input models. Then, a random effect model and a hierarchical normal
model were used to do inference on the system mean response.

Since each simulation run could be computationally expensive, an equation-
based metamodel as a function of input parameters could efficiently propagate
the input uncertainty to output. Ng and Chick (2006) developed a first-order
metamodel based on a Taylor series approximation. This local approximation is
suitable to the situations when there is a large amount of real-world data and
posterior distributions locate in a small neighborhood of the true input
parameters. To account for more general situations when the amount of real-
world data could be small, Xie et al. (2014) built a Gaussian process global
metamodel to propagate the input uncertainty to the output mean.

Chick (2001) proposed BMA approach to account for both input model and
parameter value uncertainty. Given a set of candidate parametric
distributions, Bayesian posteriors of the input families and parameters are
derived to quantify the input uncertainty. Then, the simulations are driven by
samples from the posteriors of input distributions, and the sample mean of
simulation outputs is used to estimate the system mean response. To separate
the relative contributions from input model, parameters and simulation
uncertainty, Zouaoui and Wilson (2004) developed a BMA-based simulation
replication algorithm and they assigned multiple simulation replications to
each sample of input distributions. Confidence intervals (CIs) are constructed
to quantify the overall variability for the posterior mean response. BMA is
further extended to input models with dependence by using Normal-to-Anything
(NORTA) (Biller and Corlu, 2011).

Although relatively new in the stochastic simulation community, the DPM model
has been extensively studied and widely applied in the statistics and machine
learning communities during the past decade; see Ferguson (1973), Lo (1984),
Escobar (1994), Ghosh and Ramamoorthi (2003), etc. In general, DPM has
demonstrated robust performance in terms of density estimation (Escobar and
West 1995, Görür and Rasmussen 2010, etc.). The Markov chain Monte Carlo
(MCMC) method enables efficient sampling of mixture distributions from the
posterior; see for example Escobar and West (1995), Neal (2000), Hanson
(2006), Kottas (2006), and Wang and Dunson (2011).

##  3 A Bayesian Nonparametric Framework

When we use simulation to assess the stochastic system performance, the output
from the $j$-th replication with input models, denoted by $F$, can be written
as

$Y_{j}(F)=\mu(F)+e_{j}(F)$

where $\mu(F)$ denotes the system mean response and $e_{j}(F)$ represents the
simulation error with mean zero and variance $\sigma_{e}^{2}(F)$. For
notational simplification, suppose that $F$ consists of a single univariate
model. Denote the underlying unknown true input model by $F^{c}$. We are
interested in the system mean response at the true input model, denoted by
$\mu^{c}\equiv\mu(F^{c})$. Without strong prior information on $F^{c}$ and
$\mu(\cdot)$, we introduce a Bayesian nonparametric framework to quantify the
overall estimation uncertainty of $\mu^{c}$.

Since the simulation output depends on the choice of input distribution $F$,
the input model failing to capture important features of $F^{c}$ can lead to
poor estimates of system performance. Thus, it is desirable to construct the
input model that can faithfully capture the important properties in the real-
world input data, including heterogeneity, multi-modality and skewness. In
Section 3.1, we present the flexible nonparametric DPM that can capture these
rich properties induced by various latent sources of uncertainty.

The underlying true input distribution $F^{c}$ is estimated by finite real-
world data of size $m$, denoted by
$\mathbf{X}_{m}\equiv\\{X_{1},X_{2},\ldots,X_{m}\\}$ with
$X_{i}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}F^{c}$. The posterior
distribution of the flexible input model derived from the Bayes’ rule can be
used to quantify the input uncertainty,

$p(F|\mathbf{X}_{m})\propto p(F)\cdot p(\mathbf{X}_{m}|F)$

where $p(F)$ characterizes our prior belief about the true input model $F^{c}$
and $p(\mathbf{X}_{m}|F)$ denotes the likelihood function of data
$\mathbf{X}_{m}$ under a generic input model $F$. Since the DPM model does not
have closed form distributions for analytical posterior analysis, we describe
Gibbs samplers in Section 3.2 to efficiently draw posterior samples of input
models,
$\mathcal{F}_{B}\equiv\\{\widetilde{F}^{(1)},\widetilde{F}^{(2)},\ldots,\widetilde{F}^{(B)}\\}$,
from $p(F|\mathbf{X}_{m})$ quantifying the input uncertainty. Then, we show
the asymptotic consistency of $p(F|\mathbf{X}_{m})$ in Section 3.3 and study
its finite sample performance in Section 4. In this paper, the notation
$\widetilde{\cdot}$ denotes posterior samples or random variables
characterizing our belief on input model or parameters.

When the Bayesian nonparametric approach is used to quantify the input
uncertainty, this poses difficulties in constructing an appropriate metamodel
as a functional of nonparametric input model. Thus, the direct simulation is
used to propagate the input uncertainty to the output. At each sample
$\widetilde{F}^{(b)}\sim p(F|\mathbf{X}_{m})$ with $b=1,2,\ldots,B$, we
generate $n_{b}$ replications and obtain the outputs
$\mathbf{Y}_{b}\equiv\\{Y_{1}(\widetilde{F}^{(b)}),Y_{2}(\widetilde{F}^{(b)}),\ldots,Y_{n_{b}}(\widetilde{F}^{(b)})\\}$.
Then, we estimate the mean response with sample mean
$\bar{Y}_{b}\equiv\bar{Y}(\widetilde{F}^{(b)})=\sum_{j=1}^{n_{b}}Y_{j}(\widetilde{F}^{(b)})/n_{b}$
and quantify the simulation uncertainty with the sampling distribution of
$\bar{Y}_{b}|\widetilde{F}^{(b)}$.

The overall uncertainty of system mean response estimation is characterized by
the conditional distribution of the compound random variable
$U\equiv\bar{Y}(\widetilde{F})$, denoted by $F_{U}(\cdot|\mathbf{X}_{m})$,
given the information obtained from the real-world data $\mathbf{X}_{m}$. In
Section 3.2, we propose a sampling procedure to construct a
$(1-\alpha^{*})100\%$ percentile empirical CrI quantifying the overall
estimation uncertainty of $\mu^{c}$, denoted by
$\mbox{CrI}=\left[\bar{Y}_{\left(\lceil(\alpha^{*}/2)B\rceil\right)},\bar{Y}_{\left(\lceil(1-\alpha^{*}/2)B\rceil\right)}\right],$based
on the order statistics
$\bar{Y}_{(1)}\leq\bar{Y}_{(2)}\leq\ldots\leq\bar{Y}_{(B)}$.

This empirical CrI accounts for both input and simulation uncertainties. We
study the asymptotic properties of this interval in Section 3.4. Define the
random variable $W\equiv\mu(\widetilde{F})$ with $\widetilde{F}\sim
p(F|\mathbf{X}_{m})$, which is the true system mean response evaluated at the
posterior sample of input model. Denote the conditional Cumulative
Distribution Function (c.d.f.) of $W$ by $F_{W}(\cdot|\mathbf{X}_{m})$. Let
$q_{W}(\gamma|\mathbf{X}_{m})\equiv\inf\\{q:F_{W}(q|\mathbf{X}_{m})\geq\gamma\\}$
be the conditional $\gamma$-quantile of $W$. We prove that given the input
data $\mathbf{X}_{m}$, as the simulation budget increases, the interval
$\left[\bar{Y}_{\left(\lceil(\alpha^{*}/2)B\rceil\right)},\bar{Y}_{\left(\lceil(1-\alpha^{*}/2)B\rceil\right)}\right]$
converges to the true underlying CrI quantifying the impact of input
uncertainty,
$[q_{W}(\alpha^{*}/2|\mathbf{X}_{m}),q_{W}(1-\alpha^{*}/2|\mathbf{X}_{m})]$.
We also show that as the size of real-world data and the simulation budget go
to infinity, it converges to the true mean system response $\mu^{c}$. In
addition, if the interval
$\left[\bar{Y}_{\left(\lceil(\alpha^{*}/2)B\rceil\right)},\bar{Y}_{\left(\lceil(1-\alpha^{*}/2)B\rceil\right)}\right]$
is too wide, the decision maker needs to know if the additional simulation
could improve the estimation accuracy of $\mu^{c}$. For this practical
consideration, we derive a variance decomposition to estimate the relative
contributions from input and simulation uncertainties, and further study the
asymptotic property of the corresponding variance components in Section 3.5.

###  3.1 Input Modeling by Dirichlet Process Mixtures

According to Lo (1984) and Barrios et al. (2013), given a kernel density
function $h(\cdot)$, the input density defined on $\mathfrak{X}$, the sample
space, from DPM can be represented as $f(x)=\int h(x|\boldsymbol{\psi})dG,$
where $\boldsymbol{\psi}$ denotes the parameters of kernel density function,
and an infinite mixture distribution $G$ on the parameter space of
$\boldsymbol{\psi}$ follows a Dirichlet process (DP), $G\sim
DP(\alpha,G_{0})$, with $G_{0}$ denoting the base distribution and $\alpha$
denoting the dispersion parameter. We say that $DP(\alpha,G_{0})$ assigns
probability on $\mathscr{G}$, the space of all mixing distribution $G$.
According to the definition of DP in Ferguson (1973), the random distribution
$G$ over any finite measurable partitions, $A_{1},\ldots,A_{r}$, of the space
of $\boldsymbol{\psi}$ follows a Dirichlet distribution,
$\big{(}G(A_{1}),\ldots,G(A_{r})\big{)}\sim\mbox{Dirichlet}\big{(}\alpha
G_{0}(A_{1}),\ldots,\alpha G_{0}(A_{r})\big{)}$. Thus, the data $X_{i}$ drawn
from DPM can be represented as

$X_{i}|\boldsymbol{\psi}_{i}\sim
h(\cdot|\boldsymbol{\psi}_{i}),\quad\boldsymbol{\psi}_{i}{\mid G}\sim G,\quad
G\sim DP(\alpha,G_{0}).$ (1)

According to Neal (2000), by integrating over $G$, we have the conditional
prior distribution for $\boldsymbol{\psi}_{i}$,

$\boldsymbol{\psi}_{i}|\boldsymbol{\psi}_{1},\ldots,\boldsymbol{\psi}_{i-1}\sim\frac{1}{i-1+\alpha}\sum_{i^{\prime}=1}^{i-1}\delta(\boldsymbol{\psi}_{i^{\prime}})+\frac{\alpha}{i-1+\alpha}G_{0}$
(2)

where $\delta(\boldsymbol{\psi})$ is the distribution concentrated at
$\boldsymbol{\psi}$.

DPM is specified by three key components: the dispersion parameter $\alpha$,
the kernel density $h(\cdot)$, and the base distribution $G_{0}$. The
parameter $\alpha$ represents the complexity of input model and the dispersion
of data. DPM with a larger value of $\alpha$ tends to generate samples of the
input density $f(\cdot)$ with higher complexity and more distinct active
components; see the explanation in Section 3.2. The appropriate value of
$\alpha$ can be inferred from the real-world data. Therefore, differing from
parametric approaches, the density generated by DPM can adapt its complexity
to the data.

The choice of the kernel density $h(\cdot)$ is based on the support of
$F^{c}$, and meanwhile it should account for the feasibility of implementation
in posterior computation. We present DPM models with three kernel densities,
including Gamma, Gaussian and Beta, which account for the real-world data that
are supported on $\mathfrak{X}$, where $\mathfrak{X}$ could be the half real
line $\Re^{+}$, the whole real line $\Re$, or a finite interval
$[a_{1},a_{2}]$ with $-\infty<a_{1}<a_{2}<\infty$. The scaled version of DPM
with Beta kernel is applicable to model continuous distributions with a known
finite support interval $[a_{1},a_{2}]$ through the transformation
$X_{i}=(X_{i}^{\prime}-a_{1})/(a_{2}-a_{1})$, where $X_{i}^{\prime}$ denotes
the raw data. Since Gamma, Gaussian and Beta distributions belong to the
exponential family and allow conjugate priors, we derive efficient samplers to
generate posterior samples of input model quantifying the input uncertainty.
Notice that even though these three kernels allow us to model many input
models commonly used in the simulation applications, we can also select other
kernels, such as exponential kernel density function.

To simplify the posterior inference and sampling, we consider the conjugate
prior $G_{0}$ for the parameters $\boldsymbol{\psi}$. For DPM with Gamma
kernel, we let $\boldsymbol{\psi}=\left(V,u\right)^{\top}$ with $V$ and $u$
denoting the shape and mean parameters. Motivated by the study on Gamma
mixture distributions in Wiper (2001), we consider a conditional conjugate
prior for $V$ and $u$,

$V\sim\mbox{Exponential}(\theta)\mbox{ and }u\sim\mbox{Inv-Gamma}(r,s).$ (3)

Equation (3) specifies $G_{0}(V,u)$ with the hyper-parameters
$\boldsymbol{\theta}_{G}=(\theta,r,s)$.

For DPM with Gaussian kernel, we let $\boldsymbol{\psi}=(u,\sigma^{2})^{\top}$
with $u$ and $\sigma^{2}$ denoting the mean and variance parameters. Following
Gelman et al. (2004), we choose the conjugate prior,

$u|\sigma^{2}\sim\mathcal{N}(u_{0},\sigma^{2}/m_{0})\mbox{ and
}\sigma^{2}/\sigma_{0}^{2}\sim\mbox{Inv-
Gamma}\left(\frac{v_{0}}{2},\frac{1}{2}\right).$ (4)

Equation (4) specifies $G_{0}(u,\sigma^{2})$ with hyper-parameters
$\boldsymbol{\theta}_{G}=(u_{0},m_{0},v_{0},\sigma_{0})^{\top}$.

For DPM with Beta kernel, we let
$\boldsymbol{\psi}=\left(\omega,\beta\right)^{\top}$ with $\omega$ and $\beta$
denoting the two shape parameters. Since the Beta distribution belongs to the
exponential family, we choose the conjugate prior,

$\omega,\beta|\lambda_{0},\lambda_{1},\lambda_{2}\propto\exp\left\\{-\lambda_{1}\omega-\lambda_{2}\beta-\lambda_{0}\log\left[\frac{\Gamma(\omega)\Gamma(\beta)}{\Gamma(\omega+\beta)}\right]\right\\}.$
(5)

Equation (5) specifies $G_{0}(\omega,\beta)$ with the hyper-parameters
$\boldsymbol{\theta}_{G}=(\lambda_{0},\lambda_{1},\lambda_{2})^{\top}$.

###  3.2 Gibbs Sampler for DPM and Uncertainty Quantification for $\mu^{c}$

For the real-world data $\mathbf{X}_{m}\equiv\\{X_{1},X_{2},\ldots,X_{m}\\}$,
each observation $X_{i}$ has the associated parameters
$\boldsymbol{\psi}_{i}$. According to Equation (2), parameters
$\boldsymbol{\psi}_{i}$ and $\boldsymbol{\psi}_{i^{\prime}}$ with
$i^{\prime}\neq i$ could take the same value. Let $K_{0}$ denote the number of
distinct values in $\\{\boldsymbol{\psi}_{i}\\}_{i=1}^{m}$, and represent the
distinct parameter values as
$\\{\boldsymbol{\psi}^{\star}_{1},\ldots,\boldsymbol{\psi}^{\star}_{K_{0}}\\}$.
Notice that $K_{0}$ is bounded by $m$. Following Muller et al. (1996), we
introduce the latent indicator variables
$\mathbf{c}=(c_{1},c_{2},\ldots,c_{m})$ associating the data
$\\{X_{1},X_{2},\ldots,X_{m}\\}$ to
$\\{\boldsymbol{\psi}^{\star}_{1},\ldots,\boldsymbol{\psi}^{\star}_{K_{0}}\\}$,
where $c_{i}=j$ if and only if
$\boldsymbol{\psi}_{i}=\boldsymbol{\psi}^{\star}_{j}$ for $i=1,2,\ldots,m$ and
$j=1,\ldots,K_{0}$. Thus, the real-world data $\mathbf{X}_{m}$ are grouped to
$K_{0}$ active components with parameters
$\\{\boldsymbol{\psi}^{\star}_{1},\ldots,\boldsymbol{\psi}^{\star}_{K_{0}}\\}$.
The parameters $\boldsymbol{\psi}$ are generated by $G_{0}$. A mixture
component specified by $\boldsymbol{\psi}$ is called active if it has at least
one observed data point from $\mathbf{X}_{m}$ associated with. Otherwise, it
is called inactive.

Since the DPM model (1) does not have closed form distributions for analytical
posterior analysis, we describe a Gibbs sampler to generate the posterior
samples of input model quantifying the input uncertainty. We first derive the
conditional posteriors required in the sampling procedure. According to Neal
(2000), by setting $X_{i}$ to be the last observation, the conditional prior
of $c_{i}$ can be derived based on Equation (2),

$\mbox{P}\big{(}c_{i}=j|\mathbf{c}^{-i})=\left\\{\begin{array}[]{ll}\frac{m_{j}^{-i}}{m+\alpha-1}\mbox{
if }\exists c_{q}=j\mbox{ for all }q\neq i\\\ \frac{\alpha}{m+\alpha-1}\mbox{
otherwise }\end{array}\right.$

for $i=1,2,\ldots,m$, where $\mathbf{c}^{-i}$ denotes all the latent variables
except $c_{i}$, and $m_{j}^{-i}$ is the number of latent variables with
$c_{q}=j$ for all $q\in\\{1,2,\ldots,m\\}$ and $q\neq i$. Then, given the
component parameters
$\Psi^{\star}\equiv\\{\boldsymbol{\psi}^{\star}_{1},\ldots,\boldsymbol{\psi}^{\star}_{K_{0}}\\}$
and $X_{i}$, by applying the Bayes’ rule, the conditional posterior is

$p(c_{i}=j|\mathbf{c}^{-i},\boldsymbol{\psi}^{\star}_{j},\alpha,X_{i})=\left\\{\begin{array}[]{ll}b_{0}\frac{m_{j}^{-i}}{m+\alpha-1}h(X_{i}|\boldsymbol{\psi}^{\star}_{j})\mbox{
if }\exists c_{q}=j\mbox{ for all }q\neq i\\\
b_{0}\frac{\alpha}{m+\alpha-1}\int h(X_{i}|\boldsymbol{\psi})dG_{0}\mbox{
otherwise }\end{array}\right.$ (6)

where $b_{0}$ is the normalizing constant.

The posterior for dispersion parameter $\alpha$ conditional on the number of
active components is $p(\alpha|K_{0})\sim p(\alpha)p(K_{0}|\alpha)$. We impose
a prior, $p(\alpha)=\mbox{Gamma}(a,b)$, on $\alpha$, with shape $a>0$ and
scale $b>0$. Thus, the hyper-parameters for $\alpha$ are
$\boldsymbol{\theta}_{\alpha}=(a,b)^{\top}$. To simplify the sampling
procedure for $p(\alpha|K_{0})$, following Escobar and West (1995), we
introduce a new random variable $\nu$ and generate $\alpha$ from
$p(\alpha|K_{0})$ by

$\displaystyle\nu|\alpha,K_{0}$ $\displaystyle\sim\mbox{Beta}(\alpha+1,m)$ (7)
$\displaystyle\alpha|\nu,K_{0}$
$\displaystyle\sim\tau\mbox{Gamma}(a+K_{0},b-\log(\nu))+(1-\tau)\mbox{Gamma}(a+K_{0}-1,b-\log(\nu)).$

where $\tau$ is defined by $\tau/(1-\tau)=(a+K_{0}-1)/[m(b-\log(\nu))]$.

Therefore, given the real-world input data $\mathbf{X}_{m}$, we provide a
sampling procedure to generate the samples of compound random variable
$U=\bar{Y}(\widetilde{F})$ and further build a percentile empirical CrI
accounting for both input and simulation estimation uncertainties. First,
based on the support of input model $F^{c}$, choose an appropriate kernel
density function $h(\cdot)$, and then specify the hyper-parameters for both
$G_{0}$ and $\alpha$ in Step 1; see Section 4.1 for the values of hyper-
parameters used in our empirical study. Then, motivated by Neal (2000), in
Step 2, we propose the Gibbs sampling procedure to generate samples from the
posterior of input distribution, $\widetilde{F}^{(b)}\sim p(F|\mathbf{X}_{m})$
with $b=1,2,\ldots,B$, accounting for the input uncertainty. At each
$\widetilde{F}^{(b)}$, run the simulations with $n_{b}$ replications, obtain
simulation outputs $\mathbf{Y}_{b}$, and record the sample mean $\bar{Y}_{b}$
in Step 3. The simulation uncertainty is characterized by the sampling
distribution $\bar{Y}_{b}|\widetilde{F}^{(b)}$, with mean
$\mu(\widetilde{F}^{(b)})$ and variance
$\sigma^{2}_{e}(\widetilde{F}^{(b)})/n_{b}$. Thus, the samples
$\\{\bar{Y}_{1},\bar{Y}_{2},\ldots,\bar{Y}_{B}\\}$ of
$U=\bar{Y}(\widetilde{F})$ with $\widetilde{F}^{(b)}\sim p(F|\mathbf{X}_{m})$
quantify both input and simulation uncertainties. We further construct a
$(1-\alpha^{*})100\%$ percentile empirical CrI,
$\left[\bar{Y}_{\left(\lceil(\alpha^{*}/2)B\rceil\right)},\bar{Y}_{\left(\lceil(1-\alpha^{*}/2)B\rceil\right)}\right]$,
quantifying the overall uncertainty of system mean performance estimation in
Step 4. Notice that the approaches proposed to improve the Gibbs sampling
efficiency for DPM through a collapse of the state space of the Markov chain
in Maceachern (1994, 1998), MacEachern and Muller (2000) could be incorporated
into our nonparametric Bayesian framework.

1\. Based on the support of $F^{c}$, choose an appropriate kernel density
function $h(\cdot)$. Then, specify hyper-parameters $\boldsymbol{\theta}_{G}$
and $\boldsymbol{\theta}_{\alpha}$ for the base distribution $G_{0}$ and the
dispersion parameter $\alpha$;

2\. Generate the posterior sample $\widetilde{F}^{(b)}\sim
p(F|\mathbf{X}_{m})$ through the Gibbs sampling for $b=1,\ldots,B$:

(2.1) Initialization:

(a) Generate $\widetilde{\alpha}\sim\mbox{Gamma}(a,b)$;

(b) Set $c_{i}=i$ for $i=1,\ldots,m$ and $\widetilde{K}_{0}=m$;

(c) Generate $\widetilde{\boldsymbol{\psi}}^{\star}_{j}\sim
p({\boldsymbol{\psi}}^{\star}_{j}|X_{j})$ with
$p(\boldsymbol{\psi}^{\star}_{j}|X_{j})\propto
p(\boldsymbol{\psi}^{\star}_{j})p(X_{j}|\boldsymbol{\psi}^{\star}_{j})$ for
$j=1,\ldots,\widetilde{K}_{0}$ by using the sampling procedure described in
Appendix A.1;

(2.2) In each Gibbs sampling iteration, there are three main steps described
as follows:

(a) For $i=1\mbox{ to }m$, generate the parameters associated with the
inactive component as
$\widetilde{\boldsymbol{\psi}}^{\star}_{\widetilde{K}_{0}+1}\sim G_{0}$, and
then sample $\widetilde{c}_{i}$ from the conditional posterior
$p(c_{i}=j|{\widetilde{\mathbf{c}}}^{-i},\widetilde{\boldsymbol{\psi}}^{\star}_{j},{\widetilde{\alpha}},X_{i})$
for $j=1,\ldots,\widetilde{K}_{0}+1$ by applying (6); see the detailed
sampling procedure for DPM with Gamma, Gaussian and Beta kernel densities in
Appendix A.1. Then, remove inactive components and update the number of active
components $\widetilde{K}_{0}$;

(b) For the $j$-th active component with $j=1,2,\ldots,\widetilde{K}_{0}$,
generate the $r$-th parameter in $\boldsymbol{\psi}^{\star}_{j}$, denoted by
$\widetilde{\psi}^{\star}_{jr}$, from the conditional posterior
$p(\psi^{\star}_{jr}|{\widetilde{\boldsymbol{\psi}}^{\star}_{j,-r}},\mathbf{X}^{j})$,
where $\boldsymbol{\psi}^{\star}_{j,-r}$ denotes the remaining parameters in
$\boldsymbol{\psi}^{\star}_{j}$ and $\mathbf{X}^{j}$ denotes all the data
associated to the $j$-th component; see the sampling procedure in Appendix
A.1;

(c) Generate $\widetilde{\alpha}$ from the posterior
$p(\alpha|\widetilde{K}_{0})$ by using Equations (7);

3\. At $\widetilde{F}^{(b)}$ for $b=1,\ldots,B$, generate input variates by
using Equation (10), run simulations with $n_{b}$ replications, and obtain the
outputs $\mathbf{Y}_{b}$. Then, record the sample mean
$\bar{Y}_{b}=\sum_{j=1}^{n_{b}}Y_{j}(\widetilde{F}^{(b)})/n_{b}$;

4\. Report a $(1-\alpha^{*})100\%$ two-sided percentile empirical CrI for
$\mu^{c}$
$\mbox{CrI}=\left[\bar{Y}_{\left(\lceil(\alpha^{*}/2)B\rceil\right)},\bar{Y}_{\left(\lceil(1-\alpha^{*}/2)B\rceil\right)}\right]$
(8) with the order statistics
$\bar{Y}_{(1)}\leq\bar{Y}_{(2)}\leq\ldots\leq\bar{Y}_{(B)}$.

Algorithm 1 The Nonparametric Bayesian Framework for Uncertainty
Quantification

Here, we present the detailed Gibbs sampling for generating posterior samples
of indicator variables $\mathbf{c}=(c_{1},c_{2},\ldots,c_{m})$, component
parameters $\Psi^{\star}$ and dispersion parameter $\alpha$. For the
initialization in Step (2.1), we generate $\widetilde{\alpha}$ from the prior
$\mbox{Gamma}(a,b)$ and assign the observed data points
$\mathbf{X}_{m}=\\{X_{1},X_{2},\ldots,X_{m}\\}$ to distinct components, which
follows the idea of hierarchical agglomerative clustering; see Chapter 15 of
Maimon and Rokach (2005). Then, each Gibbs sampling iteration includes three
main parts. In Step (2.2.a), the conditional posterior in (6) is used to
generate the sample of $\mathbf{c}$. For each observation $X_{i}$ with
$i=1,2,\ldots,m$ in the real-world data $\mathbf{X}_{m}$, since the
integration $\int h(X_{i}|\boldsymbol{\psi})dG_{0}$ in Equation (6) is often
intractable, we first sample the parameters for inactive component
$\widetilde{\boldsymbol{\psi}}_{\widetilde{K}_{0}+1}\sim G_{0}$. Then, we
update the latent indicator $c_{i}$ conditional on all the other parameters,
remove the inactive components, and update the number of active distinct
components $\widetilde{K}_{0}$. In Step (2.2.b), for each active component,
given the data associated to it, we update its parameters
$\boldsymbol{\psi}^{\star}_{j}$ with $j=1,2,\ldots,\widetilde{K}_{0}$. The
conditional posteriors and sampling procedure of $\mathbf{c}$ and
$\boldsymbol{\psi}^{\star}_{j}$ for Gamma, Gaussian and Beta kernel densities
are described in Appendix A.1. In Step (2.2.c), we update the dispersion
parameter $\alpha$ conditional on the number of active components and generate
the posterior sample by applying Equations (7).

After the convergence of Gibbs sampling, we can record a posterior sample of
input distribution $\widetilde{F}$ specified by the sequence of parameters
corresponding to real-world data, denoted by
$\widetilde{\Psi}\equiv\\{\widetilde{\boldsymbol{\psi}}_{1},\ldots,\widetilde{\boldsymbol{\psi}}_{m}\\}$
with
$\widetilde{\boldsymbol{\psi}}_{i}=\widetilde{\boldsymbol{\psi}}^{\star}_{\widetilde{c}_{i}}$
for $i=1,2,\ldots,m$ and the dispersion parameter $\widetilde{\alpha}$. In our
empirical study, we use 500 burn-in iterations to generate stable posterior
samples. Then, we can record $B$ posterior samples of input models after
convergence,
$\\{\widetilde{F}^{(1)},\widetilde{F}^{(2)},\ldots,\widetilde{F}^{(B)}\\}$,
quantifying the input uncertainty. The density function of
$\widetilde{F}^{(b)}$ with $b=1,2,\ldots,B$ can be represented as

$\widetilde{f}^{(b)}(x)=\frac{1}{m+\widetilde{\alpha}^{(b)}}\sum_{i=1}^{m}h(x|\widetilde{\boldsymbol{\psi}}^{(b)}_{i})+\frac{\widetilde{\alpha}^{(b)}}{m+\widetilde{\alpha}^{(b)}}\int
h(x|\boldsymbol{\psi})dG_{0}.$ (9)

Notice that the number of active components, $\widetilde{K}_{0}$, can vary at
different samples of input model.

In Step 3, to estimate the system mean response at the posterior sample of
input model $\widetilde{F}^{(b)}$ with $b=1,2,\ldots,B$, given the parameters
$\widetilde{\Psi}^{(b)}$ of active components obtained by the Gibbs sampling
procedure, we can generate input variates,

$\boldsymbol{\psi}|\widetilde{\Psi}^{(b)}\sim\frac{1}{m+\widetilde{\alpha}^{(b)}}\sum_{i=1}^{m}\delta\Big{(}\widetilde{\boldsymbol{\psi}}^{(b)}_{i}\Big{)}+\frac{\widetilde{\alpha}^{(b)}}{m+\widetilde{\alpha}^{(b)}}G_{0}\quad\mbox{and}\quad
X|\boldsymbol{\psi}\sim h(\cdot|\boldsymbol{\psi}),$ (10)

to drive the simulation and estimate the mean response by $\bar{Y}_{b}$.

We need $B$ to be reasonably large to accurately estimate the percentile CrI.
In the empirical study, we set $B=1000$ (Xie et al., 2014). Without any prior
information about the mean response $\mu(\cdot)$, in this paper, we assign
equal replications to all samples of input distribution
$\\{\widetilde{F}^{(1)},\widetilde{F}^{(2)},\ldots,\widetilde{F}^{(B)}\\}$.
Since each simulation run can be computationally expensive, a sequential
design of experiments could efficiently use the computational budget and
reduce the impact of simulation estimation uncertainty on the system
performance by finding the optimal setting for $(B,n_{1},n_{2},\ldots,n_{B})$
(Yi and Xie, 2017).

###  3.3 Posterior Consistency of Input Models

In the Bayesian paradigm, a very basic requirement is the posterior
consistency at the true input distribution (Ghosal et al., 1999b, 1995). It
means that as the amount of real-world data increases, the posterior becomes
more and more concentrated around $F^{c}$ with probability approaching 1. The
posterior consistency for DPM is studied in the statistics literature, such as
Ghosal et al. (1999b), Tokdar (2006), Wu and Ghosal (2008), etc. Given the
prior set as in equations (1), and the base measures $G_{0}$ for Dirichlet
Process set in Equations (3), (4) or (5), Theorem 1 summarizes posterior
consistency results on DPM with Gamma, Gaussian and Beta density functions as
kernels for input distributions supported on $[0,\infty)$, $\Re$ and fixed
interval $[a_{1},a_{2}]$.

The posterior consistency in Theorem 1 is stated in the sense of weak
consistency. The posterior distribution $p(\cdot\mid\mathbf{X}_{m})$ is said
to be weekly consistent at $F^{c}$ (or $f^{c}$), if with
$P_{f^{c}}$-probability 1, $p(U\mid\mathbf{X}_{m})\to 1$ for all weak
neighborhoods $U$ of $f^{c}$. We defer the detailed introduction and
definitions of weak neighborhood, weak consistency, and other related concepts
in classic Bayesian nonparametric theory to Appendix B. We also refer the
readers to Chapter 4 of Ghosh and Ramamoorthi (2003) and Chapter 7 of Ghosal
and van der Vaart (2017) for detailed discussions on posterior consistency.

######  Theorem 1.

Let $\mathbf{X}_{m}\equiv\\{X_{1},X_{2},\ldots,X_{m}\\}$ with
$X_{i}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}F^{c}$ for $i=1,2,\ldots,m$.

  * (i)

(Modified from Wu and Ghosal 2008 Theorem 14) Suppose the DPM with Gamma
kernel has the prior specified as Equation (3). Let $f^{c}$ be a continuous
and bounded density with support on $[0,\infty)$ satisfying the following
conditions: (a) $f^{c}(x)>0$ for all $x\in(0,+\infty)$ and $f^{c}(x)\leq
C_{f}<\infty$ for some finite constant $C_{f}$ for all $x\in[0,+\infty)$; (b)
$|\int_{0}^{\infty}f^{c}(x)\log f^{c}(x)dx|<\infty$; (c)
$\int_{0}^{\infty}f^{c}(x)\log\frac{f^{c}(x)}{\phi_{\delta}(x)}dx<\infty$ for
some $\delta>0$, where $\phi_{\delta}(x)=\inf_{[x,x+\delta)}f^{c}(t)$ if
$0<x<1$ and $\phi_{\delta}(x)=\inf_{(x-\delta,x]}f^{c}(t)$ if $x\geq 1$; (d)
there exists $\zeta>0$ such that
$\int_{0}^{\infty}\max(x^{-\zeta-2},x^{\zeta+2})f^{c}(x)dx<\infty$. Then, the
posterior $p(F|\mathbf{X}_{m})$ from DPM with Gamma kernel is weakly
consistent at $F^{c}$.

  * (ii)

(Tokdar 2006 Theorem 3.3) Suppose the DPM with Gaussian kernel has the prior
specified as Equation (4). Let $F^{c}$ (and the density $f^{c}$) be supported
on $\Re$ and assume that it satisfies the following conditions: (a)
$\left|\int_{-\infty}^{+\infty}f^{c}(x)\log f^{c}(x)dx\right|<+\infty$; (b)
there exists an $\eta\in(0,1)$, such that
$\int_{-\infty}^{+\infty}|x|^{\eta}f^{c}(x)dx<+\infty$; (c) there exist
constants $\sigma_{1}>0$, $c_{1}\in(0,\eta)$, $c_{2}>c_{1}$, $b_{1},b_{2}>0$,
such that for the base measure $G_{0}(u,\sigma)$ and for all large $x>0$:

$\displaystyle\max\left\\{G_{0}\left([x-\sigma_{1}x^{\eta/2},+\infty)\times[\sigma_{1},+\infty)\right),G_{0}\left([0,+\infty)\times(x^{1-\eta/2},+\infty)\right)\right\\}\geq
b_{1}x^{-c_{1}};$
$\displaystyle\max\left\\{G_{0}\left((-\infty,-x+\sigma_{1}x^{\eta/2}]\times[\sigma_{1},+\infty)\right),G_{0}\left((-\infty,0]\times(x^{1-\eta/2},+\infty)\right)\right\\}\geq
b_{1}x^{-c_{1}};$ $\displaystyle
G_{0}\left((-\infty,x)\times(0,e^{x^{\eta}-1/2})\right)>1-b_{2}x^{-c_{2}};\quad
G_{0}\left((-x,+\infty)\times(0,e^{x^{\eta}-1/2})\right)>1-b_{2}x^{-c_{2}}.$

Then, the posterior $p(F|\mathbf{X}_{m})$ from DPM with Gaussian kernel is
weakly consistent at $F^{c}$.

  * (iii)

Suppose the DPM with Beta kernel has the prior specified as Equation (5). Let
$F^{c}$ (and the density $f^{c}$) be supported on $[a_{1},a_{2}]$ and assume
that $f^{c}(x)$ is continuous density on $[a_{1},a_{2}]$. Then, the posterior
$p(F|\mathbf{X}_{m})$ from DPM with Beta kernel is weakly consistent at
$F^{c}$.

The proof of Theorem 1 is given in Appendix B. This theorem indicates that the
posterior from DPM with Gamma, Gaussian and Beta kernels can consistently
estimate the true input distributions under some conditions, including the
existence of moments and entropy of $F^{c}$ (or $f^{c}$), and the boundedness
and continuity of $f^{c}$. No assumptions on the analytic forms of $F^{c}$ and
$f^{c}$ are required for the posterior consistency. In Appendix B, we give
examples of posterior consistency on commonly used distributions in
simulations, such as normal, logistic, Student’s $t$, Cauchy, uniform,
triangular, power function, Beta, truncated normal, log-normal, Gamma with
shape parameter greater than 2, Weibull with shape parameter greater than 3,
log-logistic with shape parameter greater than 2, Pearson Type V and Type VI,
Johnson $S_{B}$, Johnson $S_{L}$, Johnson $S_{U}$, by applying Theorem 1,
though one of the main motivations for using DPM is to flexibilly model the
underlying distribution that does not belong to any known or commonly used
distribution families.

Theorem 1 gives consistency results for DPM with Gaussian, Beta and Gamma
density kernels. These three density kernels are examples that can cover a
wide range of distributions. Just like any other statistical models, they
cannot cover all situations. Some distributions, such as standard (unshifted)
versions of the gamma, Pearson VI, Pearson V, and log-logistic with shape
parameter less than or equal to 2, and Weibull distributions with shape
parameter less than or equal to 3, do not satisfy the finite moment condition
or bounded density requirement in Theorem 1 (i). This does not necessarily
mean that the proposed nonparametric approach cannot consistently estimate
such true distributions. Theorem 1 provides only sufficient conditions on
which posterior consistence holds, while our empirical study indicates that
the performance of DPM with Gamma kernel is robust to the conditions in Part
(i). Theoretically, by choosing different kernel functions, the DPM can lead
to consistent posterior on the distributions that are not covered by Theorem
1. For example, if the true density is completely monotone on $[0,\infty)$,
such as Gamma, Pearson Type V, log-logistic, and Weibull distributions with
shape parameters less than or equal to 2 or 3, the DPM priors with exponential
density or scaled uniform density kernels have consistent posterior on them,
as shown by Theorems 16 and 17 in Wu and Ghosal (2008).

Part (iii) of Theorem 1 shows the asymptotic consistency of
$p(F|\mathbf{X}_{m})$ for DPM with Beta kernel. Existing Bayesian asymptotic
results in the literature mainly focus on slightly different versions of Beta
mixtures, such as the finite mixtures of Bernstein polynomials (Petrone and
Wasserman 2002, Wu and Ghosal 2008), or the finite Beta mixtures in Rousseau
(2010). Ghosal et al. (2008) contains partial results on the classes of
distributions that can be expressed with an infinite mixture of Betas. We
presented a general result for posterior consistency of DPM with Beta kernel
and its proof is given in Appendix B. Our empirical study also demonstrates
the flexibility and adaptiveness of DPM with Beta kernel for fitting the
continuous distributions with known bounded support.

Theorem 1 directly applies to the situation that the support of the input
distribution is known. As pointed out by Law (2015), such assumption could be
a limitation for some cases. These distributions include, but not limited to,
the shifted (scaled) versions of the distributions on $[0,\infty)$ and
$[0,1]$. To deal with bounded or semi-bounded supported underlying input
distributions with unknown boundaries, we mildly extend the DPM prior as
follows: Let $\xi$ denote the boundary value(s) of the support. For the
bounded support case, we take $\xi=(a_{1},a_{2})$, while for the half-bounded
case, we take $\xi=a_{0}$ and let the support be $[a_{0},\infty)$ without loss
of generality. We assign the prior $\pi$ on $\xi$ and the complete
nonparametric prior on the input distribution is given by:

$X_{i}|\boldsymbol{\psi}_{i}\sim
h(\cdot|\boldsymbol{\psi}_{i}),\quad\boldsymbol{\psi}_{i}\mid G\sim G,\quad
G\mid\xi\sim DP(\alpha,G_{0})\times 1\\!\mathrm{l}(\xi),\quad\xi\sim\pi.$ (11)

The mixing distribution $G$ here is slightly different from the one defined in
Model (1) and it assigns point mass $1$ on the boundary $\xi$ for given $\xi$,
where $\xi$ is assigned a prior $\pi$ as an index parameter.

With the prior defined as (11), we have corollaries as follows. Please refer
to Appendix B Remark 3 for the detailed discussion and proofs.

######  Corollary 1.

Suppose that the true density $f^{c}$ is continuous, and has bounded support
with unknown boundary $\xi=(a_{1},a_{2})$. Let the DPM prior with location-
scale transformed beta kernel as described above. Assume that the prior on the
index parameter (boundary values) satisfies that for any $\delta>0$,
$\pi([a_{1}-\delta,a_{1}]\times[a_{2},a_{2}+\delta])>0$. Then the posterior
$p(F|\mathbf{X}_{m})$ from DPM with transformed Beta kernel is weakly
consistent at $F^{c}$.

######  Corollary 2.

Suppose the true density $f^{c}$ satisfies all the conditions in Theorem 1
Part (i), except that the support is $[a_{0},\infty)$ with $a_{0}$ unknown.
Let the DPM prior with shifted Gamma density kernel as described in (69).
Assume the prior on the index parameter (boundary value) satisfies
$\pi([a_{0}-\delta,a_{0}])>0$.Then the posterior $p(F|\mathbf{X}_{m})$ from
DPM with the shifted Gamma kernel is weakly consistent at $F^{c}$.

###  3.4 Asymptotic Properties of the CrI

In this section, we study the asymptotic properties of the empirical CrI
constructed from our framework in Section 3.2. In many situations, it could be
hard or expensive to collect more real-world data when we make decisions.
Therefore, in Theorem 2 part (i), we show that given finite real-world data
$\mathbf{X}_{m}$, as the simulation budget increases, the interval constructed
by our approach,
$\left[\bar{Y}_{\left(\lceil(\alpha^{*}/2)B\rceil\right)},\bar{Y}_{\left(\lceil(1-\alpha^{*}/2)B\rceil\right)}\right]$,
converges to the $(1-\alpha^{*})100\%$ percentile CrI induced by the input
uncertainty with the true mean response surface $\mu(\cdot)$ known,
$\left[q_{W}(\alpha^{*}/2|\mathbf{X}_{m}),q_{W}(1-\alpha^{*}/2|\mathbf{X}_{m})\right]$.
In Theorem 2 part (ii), we show that as the amount of real-world data and the
simulation budget go to infinity, the constructed CrI
$\left[\bar{Y}_{\left(\lceil(\alpha^{*}/2)B\rceil\right)},\bar{Y}_{\left(\lceil(1-\alpha^{*}/2)B\rceil\right)}\right]$
shrinks to the true mean response $\mu^{c}$.

The convergence between two credible intervals in Theorem 2 is measured by the
Hausdorff distance, denoted by $d_{H}(\cdot,\cdot)$, which is widely used for
measuring the distance between two sets. It has a simplified expression when
$A_{1}$ and $A_{2}$ are both closed intervals: If $A_{1}=[a_{1},b_{1}]$ and
$A_{2}=[a_{2},b_{2}]$, then
$d_{H}(A_{1},A_{2})=\max(|a_{1}-a_{2}|,|b_{1}-b_{2}|)$. In this case, the
convergence under Hausdoff distance is the same as the point-wise convergence
for the two endpoints of CrIs.

For two generic distributions (and measures) $F_{1}$ and $F_{2}$ on $\Re$ with
the Borel sigma algebra $\mathcal{B}(\Re)$, their Lévy-Prokhorov (L-P)
distance (Billingsley 1999) is defined by
$d_{LP}(F_{1},F_{2})\equiv\inf\\{\eta>0~{}|~{}F_{1}(A)\leq
F_{2}(A^{\eta})+\eta\text{ and }F_{2}(A)\leq F_{1}(A^{\eta})+\eta,\text{ for
all }A\in\mathcal{B}(\Re)\\}$, where $A^{\eta}\equiv\\{a\in\Re~{}|~{}\exists
b\in A,|a-b|<\eta\\}$. The L-P distance, denoted by $d_{LP}$, is a metric
under which the convergence is equivalent to the weak convergence of measures
on $\Re$.

######  Theorem 2.

Let $n_{\min}=\min\\{n_{1},n_{2},\ldots,n_{B}\\}\geq 1$. Suppose that the
following conditions hold:

  1. (1)

The posterior distribution function $F_{W}(\cdot|\mathbf{X}_{m})$ is
continuous with a positive density on its support; Furthermore,
$\int_{\Re}w^{2}dF_{W}(w|\mathbf{X}_{m})<\infty$ holds almost surely for all
$m$;

  2. (2)

The simulation model satisfies
$\text{E}(\bar{Y}_{b}|\widetilde{F}^{(b)})=\mu(\widetilde{F}^{(b)})$ and
$\text{Var}(\bar{Y}_{b}|\widetilde{F}^{(b)})=\sigma^{2}_{e}(\widetilde{F}^{(b)})/n_{b}$
for $b=1,\ldots,B$. For almost surely all $\widetilde{F}\sim
p(F|\mathbf{X}_{m})$, there exists a finite constant $C_{\sigma}>0$ with
$\sigma^{2}_{e}(\widetilde{F})\leq C^{2}_{\sigma}$;

  3. (3)

For any $\epsilon>0$, there exists a finite $\delta>0$ such that
$|\mu(F)-\mu(F^{c})|<\epsilon$ if $d_{LP}(F,F^{c})<\delta$;

  4. (4)

The posterior distribution $p(F|\mathbf{X}_{m})$ is weakly consistent at
$F^{c}$.

Then,

  1. (i)

If Conditions (1) and (2) hold, then for some absolute constant $C_{1}>0$ that
does not depend on $\mathbf{X}_{m}$ and $\mu(\cdot)$, the CrI in Equation (8)
satisfies

$\displaystyle\text{E}\left[\int_{0}^{1}d_{H}\left(\left[\bar{Y}_{\left(\lceil(\alpha^{*}/2)B\rceil\right)},\bar{Y}_{\left(\lceil(1-\alpha^{*}/2)B\rceil\right)}\right],\left[q_{W}\left(\alpha^{*}/2|\mathbf{X}_{m}\right),q_{W}\left(1-\alpha^{*}/2|\mathbf{X}_{m}\right)\right]\right)d\alpha^{*}~{}\Big{|}~{}\mathbf{X}_{m}\right]$
$\displaystyle\leq\frac{2C_{\sigma}}{\sqrt{n_{\min}}}+\frac{2C_{1}\int_{\Re}w^{2}dF_{W}(w|\mathbf{X}_{m})}{\sqrt{B}}.$
(12)

Furthermore, for any given $\epsilon>0$, $\delta>0$, and fixed
$\alpha^{*}\in(0,1)$, there exist some integers $B_{1}$ and $n_{\min,1}$ that
only depends on $\epsilon,\delta,\alpha^{*},C_{\sigma},\mathbf{X}_{m}$ and the
function $\mu(\cdot)$, such that for all $B>B_{1}$ and $n_{\min}>n_{\min,1}$,

$\displaystyle\text{P}\Big{(}d_{H}\Big{(}\big{[}\bar{Y}_{\left(\lceil(\alpha^{*}/2)B\rceil\right)},\bar{Y}_{\left(\lceil(1-\alpha^{*}/2)B\rceil\right)}\big{]},\big{[}q_{W}\left(\alpha^{*}/2|\mathbf{X}_{m}\right),q_{W}\left(1-\alpha^{*}/2|\mathbf{X}_{m}\right)\big{]}\Big{)}<\delta~{}\Big{|}~{}\mathbf{X}_{m}\Big{)}$
$\displaystyle>1-\epsilon.$ (13)

  2. (ii)

If Conditions (1)–(4) hold, then for any given
$\epsilon>0,\delta>0,\gamma\in(0,1)$, $\eta\in(0,2\min\\{\gamma,1-\gamma\\})$,
there exist sufficiently large integers $B_{2}$ and $n_{\min,2}$ that only
depend on $\epsilon,\delta,\gamma,\eta,C_{\sigma},\mathbf{X}_{m}$ and the
function $\mu(\cdot)$, and a sufficiently large integer $M_{0}$ that depends
on $\epsilon,\delta,\eta$, such that for all $B>B_{2}$, $n_{\min}>n_{\min,2}$,
and $m>M_{0}$,

$\displaystyle\text{P}_{F^{c}}\left[\text{P}\left(\left|\bar{Y}_{(\lceil\gamma
B\rceil)}-\mu(F^{c})\right|<\delta~{}\Big{|}~{}\mathbf{X}_{m}\right)>1-\eta\right]>1-\epsilon.$
(14)

Condition (1) is a mild condition that implies that the posterior distribution
function of the system response $\mu(\widetilde{F})$ with $\widetilde{F}\sim
p(F|\mathbf{X}_{m})$ is continuous and strictly increasing in its support,
with finite second moment. Condition (2) requires that the simulation errors
have a bounded variance. Condition (3) is about the continuity of the system
response $\mu(F)$ with respect to $F$ around $F^{c}$ in terms of the L-P
distance that can be used to characterize the weak posterior consistency. The
similar continuity assumption is commonly used in the literature on input
uncertainty and the Gaussian process metamodel when the parametric family of
input model is known; see for example Ankenman et al. (2010), Barton et al.
(2014) and Xie et al. (2014). Condition (3) generalizes it to the
nonparametric situations. Condition (4) is a direct consequence from Theorem 1
which only provides the asymptotic consistency for input models with support
on $\Re^{+}$ and $\Re$.

Given finite real-world data $\mathbf{X}_{m}$, Part (i) of Theorem 2 shows
that as the simulation budget goes to infinity with
$n_{\min},B\rightarrow\infty$, the empirical CrI obtained by our framework in
Equation (8) converges to the true underlying CrI
$[q_{W}(\alpha^{*}/2|\mathbf{X}_{m}),q_{W}(1-\alpha^{*}/2|\mathbf{X}_{m})]$.
This convergence happens in the integrated sense, in which we take an average
of the Hausdorff distance over the significance level $\alpha^{*}\in(0,1)$.
For finite $B$ and $n_{\min}$, the upper bound in ((i)) further provides a
detailed breakdown of the approximation error from the simulation estimation
uncertainty. The first error term in ((i)) comes from the finite replications
($n_{\min}$) allocated to the posterior samples of input model quantifying the
input uncertainty. The second error term in ((i)) comes from using finite
($B$) posterior samples. The convergence of CrI in ((i)) is stated in the
Bayesian setup conditional on the real-world data $\mathbf{X}_{m}$ and it does
not require the sample size $m\rightarrow\infty$. Therefore, the bound in Part
(i) is non-asymptotic in $m$ and only asymptotic in the simulation budget
$(n_{\min},B)$. The relation ((i)) further expresses this relation using the
convergence in posterior probability for each fixed $\alpha^{*}\in(0,1)$.
Notice that Part (i) only requires Conditions (1) and (2).

Part (ii) shows that the convergence of quantile as the amount of real-world
data $m$ increases to infinity, for which we have used the continuity of
$\mu(\cdot)$ at $F^{c}$ from Condition (3) and the weak consistency of the
posterior distribution $p(F|\mathbf{X}_{m})$ in Condition (4). The detailed
proof of Theorem 2 is provided in Appendix C.

###  3.5 Variance Decomposition

Following the procedure in Section 3.2, we obtain samples of system mean
response, $\bar{Y}_{b}$ with $b=1,2,\ldots,B$, quantifying the estimation
uncertainty of $\mu^{c}$. We derive the variance decomposition in Theorem
3(i). The overall uncertainty can be written as the sum of two variance
components, $\sigma^{2}_{T}=\sigma^{2}_{I}+\sigma^{2}_{S}$, quantifying the
relative contributions from input and simulation uncertainties, which can
guide how to improve the system performance estimation when the overall
uncertainty is too large. If $\sigma_{S}^{2}$ is dominant, more simulation
resource can be invested to reduce the estimation uncertainty of $\mu^{c}$.

We further study the asymptotic property of both variance components
$\sigma^{2}_{I}$ and $\sigma_{S}^{2}$ in Theorem 3(ii), which depends on
$\mathbf{X}_{m}$. The component $\sigma_{I}^{2}$ measuring the impact from
input uncertainty decreases as the amount of real-world data increases. For
the input model with support on $\Re^{+}$, $\Re$ or $[a_{1},a_{2}]$ satisfying
the conditions in Theorem 1, as $m\rightarrow\infty$, the impact of input
uncertainty disappears or $\sigma_{I}^{2}$ converges to zero in probability
$\sigma_{I}^{2}\stackrel{{\scriptstyle p}}{{\rightarrow}}0$, since the
posterior of input model $p(F|\mathbf{X}_{m})$ converges to $F^{c}$ and
$\mu(\cdot)$ is bounded and continuous around $F^{c}$ in terms of L-P
distance. As $n_{\min}\rightarrow\infty$, the variance component
$\sigma_{S}^{2}$ measuring the impact from the simulation uncertainty
disappears $\sigma_{S}^{2}\stackrel{{\scriptstyle p}}{{\rightarrow}}0$, if
$\sigma^{2}_{e}(\cdot)$ is bounded. The detailed derivation of Theorem 3 is
provided in Appendix D.

######  Theorem 3.

At any $\widetilde{F}^{(b)}$ with $b=1,2,\ldots,B$, let
$\mu_{b}=\mu(\widetilde{F}^{(b)})$ and
$\sigma_{b}^{2}=\sigma^{2}_{e}(\widetilde{F}^{(b)})$.

  * (i)

Given $\mathbf{X}_{m}$, the total variance of $\bar{Y}(\widetilde{F})$ can be
decomposed as

$\displaystyle\mbox{Var}(\bar{Y}(\widetilde{F})|\mathbf{X}_{m})=\mbox{E}_{\widetilde{F}^{(b)}}\left[\left.\frac{\sigma_{b}^{2}}{n_{b}}\right|\mathbf{X}_{m}\right]+\mbox{Var}_{\widetilde{F}^{(b)}}[\mu_{b}|\mathbf{X}_{m}]$
(15)

On the right side of equation,
$\sigma_{S}^{2}\equiv\mbox{E}_{\widetilde{F}^{(b)}}\left[\left.\frac{\sigma_{b}^{2}}{n_{b}}\right|\mathbf{X}_{m}\right]$
and
$\sigma_{I}^{2}\equiv\mbox{Var}_{\widetilde{F}^{(b)}}[\mu_{b}|\mathbf{X}_{m}]$
measure the impacts from simulation and input uncertainties. Since the sample
mean and variance $\bar{Y}_{b}$ and $S^{2}_{b}$ are the consistent estimators
for $\mu_{b}$ and $\sigma^{2}_{b}$, we estimate $\sigma_{S}^{2}$ with
$\widehat{\sigma}_{S}^{2}=\frac{1}{B}\sum_{b=1}^{B}\frac{S_{b}^{2}}{n_{b}}$
and estimate $\sigma_{I}^{2}$ with
$\widehat{\sigma}^{2}_{I}=\frac{1}{B}\sum_{b=1}^{B}(\bar{Y}_{b}-\bar{\bar{Y}})^{2}$,
where $\bar{\bar{Y}}=\frac{1}{B}\sum_{b=1}^{B}\bar{Y}_{b}$.

  * (ii)

Suppose that Conditions (2)–(4) in Theorem 2 hold. For almost surely all
$\widetilde{F}\sim p(F|\mathbf{X}_{m})$, there exists a finite constant
$C_{\mu}>0$ such that $|\mu(\widetilde{F})|\leq C_{\mu}$. Then, as $m$ and
$n_{\min}$ go to infinity, the variance components $\sigma_{I}^{2}$ and
$\sigma_{S}^{2}$ converge to zero in probability: (a)
$\sigma_{I}^{2}\stackrel{{\scriptstyle p}}{{\rightarrow}}0$ as
$m\rightarrow\infty$; and (b) $\sigma_{S}^{2}\stackrel{{\scriptstyle
p}}{{\rightarrow}}0$ as $n_{\min}\rightarrow\infty$.

##  4 Empirical Study

We first study the finite-sample performance of nonparametric input models by
using simulated data in Section 4.1 and real RM demand data collected from the
bio-pharmaceutical manufacturing in Appendix F. Results demonstrate that DPM
with appropriate kernel can capture the important properties in real-world
data and it has better and more robust finite-sample performance than existing
approaches, including finite mixture, empirical distribution, KDE and
parametric distributions selected by using the Anderson-Darling (AD) and
Kolmogorov-Smirnov (KS) tests. Since some test examples in Section 4.1 violate
the conditions in Theorem 1, results also indicate that the performance of DPM
is robust to the conditions required for input asymptotic consistency. Then,
we use an $M/G/1$ queue to study the performance of our Bayesian nonparametric
framework in Sections 4.2. Results show that our approach has good and robust
performance when there is no strong prior information on the input model and
the mean response surface. As the amount of real-world data and the simulation
budget increase, the empirical CrI
$\left[\bar{Y}_{\left(\lceil(\alpha^{*}/2)B\rceil\right)},\bar{Y}_{\left(\lceil(1-\alpha^{*}/2)B\rceil\right)}\right]$
shrinks closer to $\mu^{c}$. Further, the ratio
$\widehat{\sigma}_{I}/\widehat{\sigma}_{S}$ provides a good indicator of the
relative contributions from both input and simulation uncertainty. We also use
an $M/M/1$ queue to illustrate the performance robustness of our framework;
see results in Appendix G.

###  4.1 Input Density Estimation

In the empirical study, a Gamma prior is used for the dispersion parameter
$\alpha$. Escobar and West (1995) recommend to choose $\alpha$ around $1$. We
use the prior, $\alpha\sim\mbox{Gamma}(a,b)$ with $a=1$ and $b=1$, which puts
a fair degree of support at values around $\alpha=1$. Our sensitivity study in
Appendix E indicates that the input model performance is not sensitive to the
values of hyper-parameters $\boldsymbol{\theta}_{\alpha}$. As for the hyper-
parameters $\boldsymbol{\theta}_{G}$ for the base distribution $G_{0}$, we use
the noninformative prior. We set $\theta=0.01$, $r=2$ and $s=2$ for DPM with
Gamma kernel density, and set $\mu_{0}=0$, $v_{0}=1.5$, $m_{0}=0.01$ and
$\sigma_{0}=1$ for DPM with Gaussian kernel density. For DPM with Beta kernel
density, we set $\lambda_{0}=1,\lambda_{1}=\lambda_{2}=0.01$.

We study the performance of our nonparametric Bayesian input models by using
simulated data generated from the test examples listed in Table 1. Example 1
is an Pareto distribution with shape $1.1$ and the support on $\Re^{+}$, which
violates Condition (i) in Theorem 1. Examples 2 and 3 are shifted Gamma and
shifted Weibull with shape less than $1$ and the support on $\Re^{+}$, which
have unknown lower endpoints of their support, and violate Condition (i) in
Theorem 1. Example 4 is log-logistic with shape less than $1$ and the support
on $\Re^{+}$, which violates Condition (i) in Theorem 1. Example 5 is a
mixture distribution of log-normal with the support on $\Re^{+}$. Example 6 is
a mixture distribution of Gumbel with the support on $\Re$. Both log-normal
and Gumbel mixtures have heavy tails. Example 7 is a mixture of Beta
distributions, which has the support on $[0,1]$.

Table 1: Test examples to study the input distribution estimation

Example 1 | Pareto | Pareto(shape = 1.1, scale = 1)
---|---|---
Example 2 | Shifted Gamma | Gamma(0.5,1) with shift = 1
Example 3 | Shifted Weibull | Weibull(shape = 0.5, scale = 1) with shift = 1
Example 4 | Log-logistic | Log-logistic(shape = 0.5, scale = 1)
Example 5 | Log-normal (L) | 0.3L(0,0.1)+0.4L(1,0.1)+0.3L(2,0.1)
Example 6 | Gumbel (Gum) | 0.3Gum(1.5,0.1)+0.4Gum(2.5,0.3)+0.3Gum(5,0.5)
Example 7 | Beta (Be) | 0.3Be(10,90)+0.4Be(20,60)+0.3Be(10,10)

For Bayesian approaches, there exist various model selection criteria,
including Bayes Factor (Kass and Raftery, 1995), Posterior predictive density
(Gelman et al., 2004), and Deviance Information Criteria (Spiegelhalter et
al., 2002). However, they are not suitable here since we consider both
frequentist and Bayesian candidates. As the KS and AD test statistics are
commonly used to study the goodness of fit in the simulation community, we use
the KS and AD criteria to study the fitting performance obtained by various
approaches. Specifically, since the underlying true input model $F^{c}$ for
examples listed in Table 1 are known, we replace the hypothesized distribution
in these test statistics with $F^{c}$ to obtain corresponding distance
measures. The KS distance records the largest vertical distance between
$F^{c}(\cdot)$ and the distribution estimated by $m$ real-world data, denoted
by $\hat{F}_{m}(\cdot)$, which could be obtained by different approaches,
including DPM with various kernel densities, empirical distribution, KDE and
parametric approaches,
$D_{m}\equiv\underset{x\in\Re}{\sup}(|F^{c}(x)-\hat{F}_{m}(x)|).$ The KS
distance puts equal weight to all $x\in\Re$. Since it is typically more
challenging to estimate the tail behavior compared to the central part, the AD
distance places more weight on the tails of $F^{c}$, $A_{m}^{2}\equiv
m\int_{-\infty}^{\infty}|F^{c}(x)-\hat{F}_{m}(x)|^{2}w(x)dF^{c}(x)$, where the
weight function is $w(x)=1/\left(F^{c}(x)(1-F^{c}(x))\right)$. Thus, the AD
distance can detect the discrepancies at the tails better.

Table 2 records the statistical behaviors of KS and AD distances ($D_{m}$ and
$A_{m}$) obtained by DPM with Gamma, Gaussian and Beta kernel densities,
finite mixture (Cheng and Currie, 2003), empirical distribution, KDE, and
parametric distributions selected based on KS and AD criteria when
$m=50,100,500$. All results are based on $N=1000$ macro-replications. In the
$i$-th macro-replication, we first draw $m$ samples, denoted by
$\mathbf{X}_{m}^{(i)}$, from $F^{c}$ listed in Table 1 to mimic the procedure
collecting $m$ “real-world data”. Then, various approaches are used to fit the
real-world data, and calculate the KS and AD distances for the fitted
distributions. In the table, “parametric (AD)” and “parametric (KS)” refer to
the parametric distributions selected based on the AD and KS statistics by
using @Risk. KDE is obtained by using the R function, kde, and the bandwidth
is selected to minimize the mean integrated squared error (Sheather and Jones,
1991). For empirical distribution, KDE and parametric approaches, we find the
fitted distributions and then record the KS and AD distances for these fitted
distributions.

Differing from these frequentist approaches that provide the point estimates
of input distribution, DPM and finite mixture are Bayesian approaches.
According to Gelman et al. (2004), the posterior predictive distribution,
defined by $f(X|\mathbf{X}_{m})=\int f(X|F)dP(F|\mathbf{X}_{m})$, is
recommended for assessing the fit of input model to the real-world data. Thus,
the posterior predictive distribution is used to calculate the KS and AD
distances. We use the Gibbs samplers described in Section 3.2 to generate 100
samples of input models with the warmup equal to 500 and save the sample for
each 10 draws. Then, we aggregate these samples of input models to obtain the
posterior predictive distribution.

Table 2: KS and AD distances obtained from DPM with Gamma, Gaussian and Beta
kernel densities, the empirical distribution, KDE and parametric distributions
selected based on KS and AD tests.

$m=50$ | Example 1 | Example 2 | Example 3 | Example 4 | Example 5 | Example 6 | Example 7
---|---|---|---|---|---|---|---
DPM with Gamma | $D_{m}$ | 0.089$\pm$0.001 | 0.116$\pm$0.001 | 0.123$\pm$0.002 | 0.072$\pm$0.001 | 0.050$\pm$0.001 | NA | 0.075$\pm$0.001
$A_{m}$ | 12.079$\pm$0.153 | 15.649$\pm$0.210 | 13.172$\pm$0.231 | 11.253$\pm$0.144 | 7.093$\pm$0.098 | NA | 8.566$\pm$0.115
DPM with Gaussian | $D_{m}$ | 0.135$\pm$0.002 | 0.133$\pm$0.002 | 0.146$\pm$0.002 | 0.098$\pm$0.001 | 0.062$\pm$0.001 | 0.071$\pm$0.001 | 0.084$\pm$0.001
$A_{m}$ | 16.894$\pm$0.250 | 17.759$\pm$0.288 | 18.658$\pm$0.319 | 14.599$\pm$0.221 | 8.861$\pm$0.105 | 5.940$\pm$0.074 | 9.864$\pm$0.129
DPM with Beta | $D_{m}$ | NA | NA | NA | NA | NA | NA | 0.068$\pm$0.001
$A_{m}$ | NA | NA | NA | NA | NA | NA | 8.608$\pm$0.096
Finite Mixture | $D_{m}$ | 0.113$\pm$0.001 | 0.230$\pm$0.003 | 0.206$\pm$0.002 | 0.142$\pm$0.002 | 0.143$\pm$0.002 | 0.094$\pm$0.001 | 0.082$\pm$0.001
$A_{m}$ | 14.324$\pm$0.218 | 21.658$\pm$0.293 | 25.023$\pm$0.324 | 18.549$\pm$0.287 | 16.127$\pm$0.251 | 9.387$\pm$0.133 | 10.528$\pm$0.146
Empirical Distribution | $D_{m}$ | 0.104$\pm$0.001 | 0.122$\pm$0.002 | 0.120$\pm$0.002 | 0.115$\pm$0.002 | 0.082$\pm$0.001 | 0.081$\pm$0.001 | 0.084$\pm$0.001
$A_{m}$ | 12.224$\pm$0.163 | 16.181$\pm$0.188 | 13.769$\pm$0.169 | 13.934$\pm$0.174 | 9.734$\pm$0.118 | 6.592$\pm$0.087 | 9.809$\pm$0.120
KDE | $D_{m}$ | 0.161$\pm$0.002 | 0.174$\pm$0.002 | 0.207$\pm$0.002 | 0.085$\pm$0.001 | 0.141$\pm$0.002 | 0.130$\pm$0.002 | 0.084$\pm$0.001
$A_{m}$ | 19.086$\pm$0.266 | 20.654$\pm$0.291 | 23.996$\pm$0.325 | 13.404$\pm$0.190 | 18.951$\pm$0.283 | 11.245$\pm$0.148 | 10.408$\pm$0.142
Parametric (KS) | $D_{m}$ | 0.153$\pm$0.002 | 0.180$\pm$0.002 | 0.191$\pm$0.002 | 0.076$\pm$0.001 | 0.147$\pm$0.002 | 0.128$\pm$0.002 | 0.110$\pm$0.001
$A_{m}$ | 22.045$\pm$0.305 | 25.541$\pm$0.328 | 26.410$\pm$0.346 | 12.679$\pm$0.170 | 19.653$\pm$0.281 | 12.273$\pm$0.168 | 11.506$\pm$0.157
Parametric (AD) | $D_{m}$ | 0.154$\pm$0.002 | 0.184$\pm$0.002 | 0.191$\pm$0.002 | 0.078$\pm$0.001 | 0.148$\pm$0.002 | 0.131$\pm$0.002 | 0.113$\pm$0.001
$A_{m}$ | 21.922$\pm$0.292 | 25.528$\pm$0.303 | 26.207$\pm$0.352 | 12.100$\pm$0.174 | 19.744$\pm$0.283 | 11.924$\pm$0.159 | 11.233$\pm$0.152
$m=100$ | Example 1 | Example 2 | Example 3 | Example 4 | Example 5 | Example 6 | Example 7
DPM with Gamma | $D_{m}$ | 0.058$\pm$0.001 | 0.088$\pm$0.001 | 0.084$\pm$0.001 | 0.051$\pm$0.001 | 0.056$\pm$0.001 | NA | 0.053$\pm$0.001
$A_{m}$ | 8.095$\pm$0.114 | 9.093$\pm$0.142 | 10.516$\pm$0.166 | 7.160$\pm$0.109 | 6.121$\pm$0.088 | NA | 6.208$\pm$0.090
DPM with Gaussian | $D_{m}$ | 0.091$\pm$0.001 | 0.100$\pm$0.001 | 0.117$\pm$0.002 | 0.069$\pm$0.001 | 0.064$\pm$0.001 | 0.062$\pm$0.001 | 0.056$\pm$0.001
$A_{m}$ | 13.174$\pm$0.187 | 12.678$\pm$0.182 | 14.224$\pm$0.202 | 10.372$\pm$0.153 | 8.133$\pm$0.111 | 4.390$\pm$0.068 | 7.334$\pm$0.105
DPM with Beta | $D_{m}$ | NA | NA | NA | NA | NA | NA | 0.049$\pm$0.001
$A_{m}$ | NA | NA | NA | NA | NA | NA | 6.416$\pm$0.092
Finite Mixture | $D_{m}$ | 0.096$\pm$0.001 | 0.176$\pm$0.002 | 0.185$\pm$0.002 | 0.064$\pm$0.001 | 0.094$\pm$0.001 | 0.066$\pm$0.001 | 0.061$\pm$0.001
$A_{m}$ | 15.550$\pm$0.207 | 17.582$\pm$0.237 | 18.188$\pm$0.255 | 11.288$\pm$0.167 | 16.583$\pm$0.214 | 5.329$\pm$0.080 | 7.802$\pm$0.102
Empirical Distribution | $D_{m}$ | 0.071$\pm$0.001 | 0.086$\pm$0.001 | 0.087$\pm$0.001 | 0.083$\pm$0.001 | 0.064$\pm$0.001 | 0.068$\pm$0.001 | 0.059$\pm$0.001
$A_{m}$ | 8.319$\pm$0.115 | 9.621$\pm$0.134 | 9.687$\pm$0.139 | 8.860$\pm$0.128 | 6.874$\pm$0.096 | 5.443$\pm$0.075 | 7.230$\pm$0.098
KDE | $D_{m}$ | 0.130$\pm$0.002 | 0.176$\pm$0.002 | 0.211$\pm$0.001 | 0.062$\pm$0.001 | 0.132$\pm$0.002 | 0.078$\pm$0.001 | 0.060$\pm$0.001
$A_{m}$ | 15.431$\pm$0.204 | 18.544$\pm$0.239 | 22.895$\pm$0.315 | 10.623$\pm$0.161 | 18.037$\pm$0.243 | 6.386$\pm$0.094 | 8.095$\pm$0.112
Parametric (KS) | $D_{m}$ | 0.141$\pm$0.002 | 0.174$\pm$0.002 | 0.182$\pm$0.002 | 0.049$\pm$0.001 | 0.146$\pm$0.001 | 0.087$\pm$0.001 | 0.086$\pm$0.001
$A_{m}$ | 19.710$\pm$0.315 | 21.593$\pm$0.330 | 22.263$\pm$0.367 | 7.752$\pm$0.118 | 16.585$\pm$0.293 | 7.794$\pm$0.115 | 9.503$\pm$0.124
Parametric (AD) | $D_{m}$ | 0.143$\pm$0.002 | 0.175$\pm$0.002 | 0.184$\pm$0.002 | 0.050$\pm$0.001 | 0.146$\pm$0.002 | 0.089$\pm$0.001 | 0.086$\pm$0.001
$A_{m}$ | 19.448$\pm$0.322 | 21.081$\pm$0.335 | 22.928$\pm$0.354 | 7.680$\pm$0.114 | 16.330$\pm$0.287 | 7.542$\pm$0.109 | 9.337$\pm$0.133

In the $i$-th macro-replication, we obtain the KS and AD distances, denoted by
$D^{(i)}_{m}$ and $A^{(i)}_{m}$, with $i=1,2,\ldots,N$. Then, we record 95%
symmetric CIs for both KS and AD distances, denoted by $\bar{D}\pm
1.96S_{D}/\sqrt{N}$ and $\bar{A}\pm 1.96S_{A}/\sqrt{N}$, in Table 2, and
highlight the smallest values, where $\bar{D}=\sum_{i=1}^{N}D^{(i)}_{m}/N$,
$\bar{A}=\sum_{i=1}^{N}A^{(i)}_{m}/N$,
$S_{D}=\big{[}\sum_{i=1}^{N}(D^{(i)}_{m}-\bar{D})^{2}/(N-1)\big{]}^{1/2}$ and
$S_{A}=\big{[}\sum_{i=1}^{N}(A^{(i)}_{m}-\bar{A})^{2}/(N-1)\big{]}^{1/2}$. As
$m$ increases, the KS and AD distances obtained from all approaches decrease,
and the DPM with appropriate kernel density typically has the best
performance. Notice that DPM with Gamma and Beta kernel densities performs
better than DPM with Gaussian kernel, which is the main focus of study in both
statistics and machine learning communities. Further, DPM with Gamma kernel
fits different input models with support on $\Re^{+}$ well. Based on the
results of AD distance, DPM tends to provide better estimation on the tail
behavior compared with the finite mixture, empirical, KDE and parametric
distributions, especially when $m$ is not large.

To study the number of components identified by DPM and the finite mixture
(Cheng and Currie, 2003), we consider a test example with $F^{c}$ equal to
$0.3\mbox{Gum}(1,0.1)+0.3{L(2,0.1)}+0.4\mathcal{N}(4,0.5)$. It includes three
components from different parametric families, where the Gumbel and log-normal
components are asymmetric and the normal component is symmetric. Here, we use
the DPM with Gaussian kernel, and record the marginal distribution for the
number of active components,

$\bar{p}(K_{0}=k)=\int p(K_{0}=k|\mathbf{X}_{m})dF^{c}(\mathbf{X}_{m}).$

To differ with the prior distribution, we use the notation $\bar{p}(K_{0})$
here. The probability $\bar{p}(K_{0}=k)$ is estimated based on $N=100$ macro-
replications and $B_{0}=100$ posterior samples of input model obtained in each
macro-replication with results shown in Table 3. In the $i$-th macro-
replication, we generate $\mathbf{X}_{m}^{(i)}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}F^{c}$ with $i=1,2,\ldots,N$. The marginal probability
$\bar{p}(K_{0}=k)$ is estimated by using
$\frac{1}{B_{0}N}\sum_{i=1}^{N}\sum_{b=1}^{B_{0}}\delta(K_{0}^{(b)}=k|\mathbf{X}_{m}^{(i)})$,
where $k=1,2,\ldots,m$ and $\delta(\cdot|\mathbf{X}_{m}^{(i)})$ denotes an
indicator function conditional on $\mathbf{X}_{m}^{(i)}$ for $i=1,2,\ldots,N$.
The posterior samples of $K_{0}$ and input model can be obtained by following
the procedure in Section 3.2. We compare the posterior of $K_{0}$ obtained by
DPM with the finite Gaussian mixture using Maximum a posteriori importance
sampling (MAPIS) described in Cheng and Currie (2003). We record the estimated
marginal probability in Table 3 when $m=50,100,500$. The DPM provides a better
detection of the underlying number sources of uncertainty.

Table 3: The Estimated Marginal Distribution for the Number of Active
Components, $\bar{p}(K_{0}=k)$

$k$ | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | >=10
---|---|---|---|---|---|---|---|---|---|---
$m=50$ | DPM | 0.014 | 0.165 | 0.453 | 0.239 | 0.084 | 0.026 | 0.013 | 0.004 | 0.002 | 0
MAPIS | 0.047 | 0.079 | 0.109 | 0.102 | 0.108 | 0.106 | 0.109 | 0.116 | 0.110 | 0.115
$m=100$ | DPM | 0 | 0.138 | 0.503 | 0.255 | 0.072 | 0.023 | 0.008 | 0.001 | 0 | 0
MAPIS | 0.036 | 0.070 | 0.114 | 0.111 | 0.118 | 0.107 | 0.112 | 0.117 | 0.115 | 0.099
$m=500$ | DPM | 0 | 0.097 | 0.558 | 0.246 | 0.077 | 0.016 | 0.006 | 0 | 0 | 0
MAPIS | 0.032 | 0.050 | 0.123 | 0.122 | 0.106 | 0.113 | 0.123 | 0.115 | 0.108 | 0.109

In addition, given the data $\mathbf{X}_{m}$, Figures 2 and 2 give the
representative posterior samples of input model obtained by DPM and finite
mixture (Cheng and Currie, 2003) when $m=500$. The posterior density function
of $\widetilde{F}^{(b)}$ is given in Equation (9), where the integration can
be estimated by $\frac{1}{N_{G}}\sum_{i=1}^{N_{G}}h(x|\boldsymbol{\psi}_{i})$
with $\boldsymbol{\psi}_{i}\sim G_{0}$, and we use $N_{G}=1000$ here. In
Figures 2 and 2, the solid line represents the true density function, and the
dashed lines represent the posterior samples of input model. The DPM can
deliver more accurate density estimation.

Figure 1: Posterior samples of input density obtained by DPM

Figure 2: Posterior samples of input density obtained by finite mixture

###  4.2 An $M/G/1$ Queue

An $M/G/1$ queue is used to study the performance of our Bayesian
nonparametric framework, and we compare it with the direct bootstrap (Barton
and Schruben, 2001). Suppose that the arrival process is known with the
arrival rate equal to $\lambda$, and the true distribution of service time
$F^{c}$ is unknown. We are interested in the probability of each customer
having the time staying in the system greater than a threshold, denoted by
$\tau$.

The unknown distribution for service time is estimated by using $m=50,500$
observations drawn from $F^{c}$. Among the seven examples in Table 1, Example
1 has infinite variance and Example 6 has support on $\Re$. Thus, we consider
the remaining five examples as shown in Table 4 as the underlying input
distribution. For each case, the true mean system response $\mu^{c}$ is
estimated by a side experiment with the runlength equal to $10^{6}$ customers.

Table 4: The settings for five $M/G/1$ test examples

Example | Distribution of Service Time | Inter-arrival Time | Threshold | $\mu^{c}$
---|---|---|---|---
Log-normal | $0.3L(0,0.1)+0.4L(1,0.1)+0.3L(2,0.1)$ | $\mbox{Exp}(\lambda=0.2)$ | $\tau=25$ | $0.0843\pm 0.0002$
Log-logistic | $\mbox{Log-logistic}(\mbox{shape}=0.5,\mbox{scale}=1)$ | $\mbox{Exp}(\lambda=0.1)$ | $\tau=20$ | $0.1123\pm 0.0003$
Shifted-Gamma |  $\mbox{Gamma}(0.5,1)$ with $\mbox{shift}=1$ | $\mbox{Exp}(\lambda=0.5)$ | $\tau=8$ | $0.1303\pm 0.0003$
Shifted-Weibull |  $\mbox{Weibull}(0.5,1)$ with $\mbox{shift}=1$ | $\mbox{Exp}(\lambda=0.25)$ | $\tau=40$ | $0.1198\pm 0.0003$
Beta | $0.3\mbox{Be}(10,90)+0.4\mbox{Be}(20,60)+0.3\mbox{Be}(10,10)$ | $\mbox{Exp}(\lambda=3)$ | $\tau=3$ | $0.0837\pm 0.0002$

For the DPM, we use the Gamma kernel. By following the sampling procedure
described in Section 3.2, we construct the 95% percentile empirical CrI,
$\left[\bar{Y}_{\left(\lceil(\alpha^{*}/2)B\rceil\right)},\bar{Y}_{\left(\lceil(1-\alpha^{*}/2)B\rceil\right)}\right]$,
accounting for both input and simulation uncertainties. We use
$\bar{\bar{Y}}=\frac{1}{B}\sum_{b=1}^{B}\bar{Y}_{b}$ as the point estimator
for the system mean response. For the direct bootstrap, we first generate $B$
sets of bootstrapped input data, and build empirical distributions,
$\widehat{F}^{(1)},\widehat{F}^{(2)},\ldots,\widehat{F}^{(B)}$, and at each
$\widehat{F}^{(b)}$ with $b=1,2,\ldots,B$, record the output sample mean
$\bar{Y}_{b}$. A percentile CI is constructed according to Barton and Schruben
(2001), and similarly, $\bar{\bar{Y}}=\frac{1}{B}\sum_{b=1}^{B}\bar{Y}_{b}$ is
used as the point estimator. We set $B=1000$, and assign equal replications,
$n=100,1000$, to each posterior or bootstrapped sample of input model. Each
simulation run starts with the empty system with both warmup and runlength
equal to $1000$ customers.

To compare the performance of our approach with direct bootstrap, we first
report the mean and standard deviation (SD) of the deviation of the point
estimator from $\mu^{c}$, defined by $\mbox{Err}=|\bar{\bar{Y}}-\mu^{c}|$, and
then record the mean and SD of the CrI and CI width as shown in Table 5, where
“DPM” and “EMP” represent our approach and direct bootstrap, respectively. The
results are based on $1000$ macro replications. The proposed framework
provides a smaller Err and a shorter CrI width. Since the frequentist and
Bayesian approaches are based on different philosophies, we do not record the
CrI and CI coverage. We also report the ratio
${\widehat{\sigma}_{I}^{2}}/{\widehat{\sigma}_{S}^{2}}$ from our approach
according to Section 3.5, which provide insights on the contributions of input
and simulation estimation uncertainties.

Table 5: Probability having the time staying in system greater than $\tau$
obtained by using our framework (DPM) and direct bootstrap (EMP) (inside the
brackets is the standard deviation)

Log-normal | Err (DPM) | CrI Width (DPM) | Err (EMP) | CI Width (EMP)
---|---|---|---|---
$m=500,n=1000$ | 0.020 (0.002) | 0.068 (0.008) | 0.026 (0.003) | 0.095 (0.011)
$m=500,n=100$ | 0.024 (0.003) | 0.073 (0.009) | 0.029 (0.003) | 0.102 (0.013)
$m=50,n=1000$ | 0.051 (0.005) | 0.206 (0.024) | 0.061 (0.006) | 0.258 (0.027)
$m=50,n=100$ | 0.055 (0.006) | 0.223 (0.023) | 0.069 (0.006) | 0.271 (0.030)
Log-logistic | Err (DPM) | CrI Width (DPM) | Err (EMP) | CI Width (EMP)
$m=500,n=1000$ | 0.038 (0.004) | 0.074 (0.009) | 0.044 (0.005) | 0.098 (0.011)
$m=500,n=100$ | 0.045 (0.005) | 0.080 (0.010) | 0.052 (0.006) | 0.105 (0.012)
$m=50,n=1000$ | 0.103 (0.011) | 0.186 (0.023) | 0.122(0.013) | 0.254 (0.031)
$m=50,n=100$ | 0.108 (0.012) | 0.198 (0.026) | 0.128 (0.014) | 0.263 (0.032)
Shifted-Gamma | Err (DPM) | CrI Width (DPM) | Err (EMP) | CI Width (EMP)
$m=500,n=1000$ | 0.043 (0.005) | 0.087 (0.011) | 0.044 (0.005) | 0.096 (0.012)
$m=500,n=100$ | 0.048 (0.005) | 0.095 (0.012) | 0.047 (0.005) | 0.101 (0.012)
$m=50,n=1000$ | 0.099 (0.011) | 0.282 (0.032) | 0.105 (0.011) | 0.296 (0.034)
$m=50,n=100$ | 0.104 (0.011) | 0.290 (0.034) | 0.109 (0.012) | 0.307 (0.038)
Shifted-Weibull | Err (DPM) | CrI Width (DPM) | Err (EMP) | CI Width (EMP)
$m=500,n=1000$ | 0.064 (0.007) | 0.125 (0.018) | 0.065 (0.007) | 0.133 (0.020)
$m=500,n=100$ | 0.069 (0.008) | 0.131 (0.019) | 0.071 (0.008) | 0.139 (0.022)
$m=50,n=1000$ | 0.184 (0.022) | 0.409 (0.045) | 0.220 (0.026) | 0.509 (0.055)
$m=50,n=100$ | 0.191 (0.023) | 0.417 (0.047) | 0.228 (0.027) | 0.522 (0.058)
Beta | Err (DPM) | CrI Width (DPM) | Err (EMP) | CI Width (EMP)
$m=500,n=1000$ | 0.015 (0.002) | 0.041 (0.006) | 0.019 (0.003) | 0.045 (0.006)
$m=500,n=100$ | 0.018 (0.002) | 0.044 (0.006) | 0.022 (0.003) | 0.049 (0.006)
$m=50,n=1000$ | 0.029 (0.004) | 0.100 (0.012) | 0.037 (0.005) | 0.121 (0.016)
$m=50,n=100$ | 0.034 (0.005) | 0.105 (0.013) | 0.041 (0.006) | 0.128 (0.016)

Table 6: Ratio of Input and Simulation Uncertainties (inside the brackets is
the standard deviation)

${\widehat{\sigma}_{I}^{2}}/{\widehat{\sigma}_{S}^{2}}$ | Log-normal | Log-logistic | Shifted-Gamma | Shifted-Weibull | Beta
---|---|---|---|---|---
$m=500,n=1000$ | 4.957 (1.238) | 8.792 (3.525) | 6.183 (1.736) | 8.034 (2.454) | 5.870 (1.416)
$m=500,n=100$ | 0.635 (0.142) | 1.455 (0.247) | 0.924 (0.280) | 1.252 (0.314) | 0.908 (0.241)
$m=50,n=1000$ | 8.492 (3.043) | 22.560 (8.784) | 11.273 (4.028) | 16.813 (6.255) | 9.405 (3.719)
$m=50,n=100$ | 1.280 (0.573) | 4.642 (1.839) | 1.765 (0.844) | 3.508 (1.493) | 2.072 (0.835)

##  5 Conclusions

Without strong prior information on the true input models and the system mean
response surface, in this paper, a Bayesian nonparametric framework is
proposed to quantify the overall uncertainty of system mean performance
estimation. The DPM can capture the important properties in the real-world
data, including heterogeneity, multi-modality, skewness and tails. The
posteriors of flexible input models can automatically account for both model
selection and parameters value uncertainty. Then, direct simulation is used to
propagate the input uncertainty to the outputs with the simulation uncertainty
quantified by the sampling distribution of system mean responses. Therefore,
given the real-world input data, our framework leads to a sampling procedure
that can deliver a conditional distribution of the system mean response and
provide a percentile empirical CrI accounting for both input and simulation
uncertainties. A variance decomposition is further developed to quantify the
relative contributions from both sources of uncertainty. Our approach is
supported with a rigorous asymptotic study. Given a finite amount of real-
world data, as the simulation budget increases, our CrI converges to the CrI
accounting for input uncertainty with the true mean response surface known. As
both real-world data and simulation budget go to infinity, our empirical CrI
converges to the true system response.

The empirical study demonstrates obvious advantages of Bayesian nonparametric
DPM for input density estimation compared to existing approaches, including
empirical distribution, KDE and parametric approaches. The simulation results
indicate that our framework is robust to the conditions required for the
asymptotic study. It demonstrates better performance than the nonparametric
direct bootstrap. The ratio $\sigma_{I}/\sigma_{S}$ provides a good measure of
the relative contributions from input and simulation uncertainties.

##  Appendix A Gibbs Samplers for DPM with Gamma, Gaussian and Beta Kernels

For DPM with Gamma, Gaussian and Beta kernels, we provide the posterior
inference and sampling for the indicator variables $\mathbf{c}$ and component
parameters $\boldsymbol{\psi}^{\star}_{j}$ for $j=1,2,\ldots,K_{0}$ used in
Steps 1 and 2 of the Gibbs samplers presented in Section 3.2. We describe the
main results to support the Gibbs sampling in Section A.1. Then, in Section
A.2, we provide the detailed derivation of the results used in the sampling
procedure.

###  A.1 Gibbs Sampling for $\mathbf{c}$ and $\Psi^{\star}$

####  A.1.1 DPM with Gamma Kernel

Here, we present a posterior sampler for the DPM with Gamma kernel. Given the
base distribution $G_{0}$ in Equation (3), we first generate samples of latent
variables $\mathbf{c}$ for Step 1 of the Gibbs sampler in Section 3.2.
According to Equation (6), the conditional posterior probabilities of $c_{i}$
in DPM with Gamma kernel is

$p(c_{i}=j|\mathbf{c}^{-i},\boldsymbol{\psi}^{\star}_{j},\alpha,X_{i})=\left\\{\begin{array}[]{ll}b_{0}\frac{m_{j}^{-i}}{m+\alpha-1}X_{i}^{V_{j}-1}e^{-\frac{V_{j}}{u_{j}}X_{i}}\mbox{
if }\exists c_{q}=j\mbox{ for all }q\neq i\\\
b_{0}\frac{\alpha}{m+\alpha-1}\int
X_{i}^{V-1}e^{-\frac{V}{u}X_{i}}dG_{0}(V,u)\mbox{ otherwise
}\end{array}\right.$ (16)

where $b_{0}$ denotes the normalizing constant. When $X_{i}$ comes from a new
component, the conditional posterior for $c_{i}$ in Equation (16) is not
analytically tractable and a sampling approach is used to generate samples of
$\mathbf{c}$ by following Algorithm 4 in Neal (2000).

Next we generate samples of the parameters
$\boldsymbol{\psi}^{\star}_{j}=\left(V_{j},u_{j}\right)^{\top}$ for Step 2 of
the Gibbs sampler. By the Bayes’ rule, $p(V_{j}|u_{j},\mathbf{X}^{j})\propto
p(V_{j})f(\mathbf{X}^{j}|V_{j},u_{j})$ and
$p(u_{j}|V_{j},\mathbf{X}^{j})\propto p(u_{j})f(\mathbf{X}^{j}|V_{j},u_{j})$,
the conditional posteriors of $V_{j}$ and $u_{j}$ are given by

$\displaystyle V_{j}|u_{j},\mathbf{X}^{j}$ $\displaystyle\propto$
$\displaystyle\frac{V_{j}^{m_{j}V_{j}}}{\Gamma(V_{j})^{m_{j}}}\exp\left[-V_{j}\left(\theta+\frac{\sum_{k=1}^{m_{j}}{X_{k}^{j}}}{u_{j}}+m_{j}\log(u_{j})-\sum_{k=1}^{m_{j}}{\log(X_{k}^{j})}\right)\right]$
(17) $\displaystyle u_{j}|V_{j},\mathbf{X}^{j}$ $\displaystyle\sim$
$\displaystyle\mbox{Inv-
Gamma}\left(r+m_{j}V_{j},s+V_{j}\sum_{k=1}^{m_{j}}X_{k}^{j}\right)$

where $X_{k}^{j}$ are the $k$th observation associated to the $j$-th component
and $m_{j}$ is the size of $\mathbf{X}^{j}$. The detailed derivation for these
posteriors can be found in Appendix A.2.1.

The conditional posterior $p(V_{j}|u_{j},\mathbf{X}^{j})$ in Equation (17) is
not a standard distribution. A Metropolis-Hasting (M-H) nested Gibbs sampler
is developed to generate samples of $V_{j}$ from the conditional posterior.
Specifically, denote the sample from the previous iteration in the nested M-H
sampling by $V_{j}^{0}$. We first generate a candidate sample
$\widetilde{V}_{j}$ from a proposal distribution, denoted by
$g(\cdot,V_{j}^{0})$, and accept it with probability

$\min\left\\{1,\frac{p(\widetilde{V}_{j}|u_{j},\mathbf{X}^{j})g(V_{j}^{0},\widetilde{V}_{j})}{p(V_{j}^{0}|u_{j},\mathbf{X}^{j})g(\widetilde{V}_{j},V_{j}^{0})}\right\\},$

where $p(V_{j}^{0}|u_{j},\mathbf{X}^{j})$ and
$p(\widetilde{V}_{j}|u_{j},\mathbf{X}^{j})$ are the conditional posteriors
from Equation (17). Otherwise, retain the value of $V_{j}^{0}$. The proposal
distribution $g(\cdot,V_{j}^{0})$ is chosen to be
$\mbox{Gamma}(d,d/V_{j}^{0})$ with mean located at $V_{j}^{0}$. This proposal
distribution is determined by using the Stirling approximation so that it can
capture the tail of the conditional posterior $p(V_{j}|u_{j},\mathbf{X}^{j})$
well. The detailed derivation can be found in Appendix A.2.1. To make the
proposal distribution relatively flat, we recommend that the value of $d$ is
set to be small, e.g., $d=2$ used in our empirical study.

####  A.1.2 DPM with Gaussian Kernel

Given the base distribution $G_{0}$ in Equation (4), we first generate samples
of the latent variables $\mathbf{c}$ for Step 1 of the Gibbs sampler. If
$c_{i}$ is associated with an existing $j$th component, then

$p(c_{i}=j|\mathbf{c}^{-i},\boldsymbol{\psi}^{\star}_{j},\alpha,X_{i})=b_{0}\frac{m_{j}^{-i}}{m+\alpha-1}\frac{1}{{\sqrt{2\pi}\sigma_{j}}}e^{{{-\left({X_{i}-u_{j}}\right)^{2}}\mathord{\left/{\vphantom{{-\left({x-u_{j}}\right)^{2}}{2\sigma_{j}^{2}}}}\right.\kern-1.2pt}{2\sigma_{j}^{2}}}}.$

If $c_{i}$ is associated with a new component, then

$p(c_{i}=j|\mathbf{c}^{-i},\boldsymbol{\psi}^{\star}_{j},\alpha,X_{i})=b_{0}\frac{\alpha}{m+\alpha-1}\frac{\left(v_{0}/2\right)^{v_{0}/2}}{\Gamma(v_{0}/2)}\sigma_{0}^{v_{0}}\sqrt{\frac{m_{0}}{2\pi(m_{0}+1)}}\frac{\Gamma(A)}{B^{A}}$

where $A=(v_{0}+1)/2$,
$B=[v_{0}\sigma_{0}^{2}+m_{0}(X_{i}-u_{0})^{2}/(m_{0}+1)]/2$ and $b_{0}$ is
the normalizing constant. The detailed derivation for this conditional
posterior can be found in Appendix A.2.2.

Next we generate samples of the parameters
$\boldsymbol{\psi}^{\star}_{j}=(u_{j},\sigma_{j}^{2})^{\top}$ for Step 2 of
the Gibbs sampler. The conditional posteriors for $u_{j}$ and $\sigma_{j}$ are
derived by following Chapter 3 in Gelman et al. (2004)

$\displaystyle u_{j}|\sigma_{j},\mathbf{X}^{j}$ $\displaystyle\sim$
$\displaystyle\mathcal{N}\left(\frac{m_{0}}{m_{0}+m_{j}}u_{0}+\frac{m_{j}}{m_{0}+m_{j}}\bar{X}^{j},\frac{\sigma_{0j}^{2}}{m_{0}+m_{j}}\right),$
$\displaystyle\sigma_{j}^{2}/\sigma_{0}^{2}\Big{|}\mathbf{X}^{j}$
$\displaystyle\sim$ $\displaystyle\mbox{Inv-
Gamma}\left(\frac{v_{0}+m_{j}}{2},\frac{1}{2}\right),$

where

$\sigma_{0j}^{2}=\frac{v_{0}\sigma_{0}^{2}+\sum_{k=1}^{m_{j}}(X^{j}_{k}-\bar{X}^{j})^{2}+\frac{m_{0}m_{j}(\bar{X}^{j}-u_{0})^{2}}{m_{0}+m_{j}}}{v_{0}+m_{j}}\mbox{
with }\bar{X}^{j}=\frac{1}{m_{j}}\sum_{k=1}^{m_{j}}X_{k}^{j}.$

####  A.1.3 DPM with Beta Kernel

Here we develop a posterior sampler for DPM with the Beta kernel density to
fit the input models with compact supports. We assume that
$X_{i}|c_{i}=j,\omega_{j},\beta_{j}\sim\mbox{Beta}(\omega_{j},\beta_{j})$ and
denote the parameters for the $j$th component by
$\boldsymbol{\psi}^{\star}_{j}=\left(\omega_{j},\beta_{j}\right)^{\top}$.
Equation (5) provides the base function $G_{0}(\omega,\beta)$. The derivation
for this prior can be founded in Appendix A.2.3.

We first generate samples of the latent variable $c_{i}$ for Step 1 of the
Gibbs sampler. According to Equation (6), the conditional posterior
probabilities of $c_{i}$ in DPM with Beta kernel is

$p(c_{i}=j|\mathbf{c}^{-i},\boldsymbol{\psi}^{\star}_{j},\alpha,X_{i})=\left\\{\begin{array}[]{ll}b_{0}\frac{m_{j}^{-i}}{m+\alpha-1}X_{i}^{\omega_{j}-1}(1-X_{i})^{\beta_{j}-1}\mbox{
if }\exists c_{q}=j\mbox{ for all }q\neq i\\\
b_{0}\frac{\alpha}{m+\alpha-1}\int
X_{i}^{\omega-1}(1-X_{i})^{\beta-1}dG_{0}(\omega,\beta)\mbox{ otherwise
}\end{array}\right.$

where $b_{0}$ denotes the normalizing constant. Since the conditional
posterior for $c_{i}$ associated with a new component does not have a closed
form, we use the sampling approach by following Algorithm 4 in Neal (2000) to
generate samples of $c_{i}$.

Next we generate samples of the parameters
$\boldsymbol{\psi}^{\star}_{j}=(\omega_{j},\beta_{j})^{\top}$ for Step 2 of
the Gibbs sampler. By applying the Bayes’ rule,
$p(\omega_{j}|\beta_{j},\mathbf{X}^{j})\propto
p(\omega_{j})p(\mathbf{X}^{j}|\omega_{j},\beta_{j})$ and
$p(\beta_{j}|\omega_{j},\mathbf{X}^{j})\propto
p(\beta_{j})p(\mathbf{X}^{j}|\omega_{j},\beta_{j})$, the conditional
posteriors of component parameters $\omega_{j}$ and $\beta_{j}$ are given by

$\displaystyle\omega_{j}|\beta_{j},\mathbf{X}^{j}$ $\displaystyle\propto$
$\displaystyle\exp\left\\{\left[-\lambda_{1}+\sum_{k=1}^{m_{j}}{\log(X_{k}^{j})}\right]\omega_{j}-(\lambda_{0}+m_{j})\log\left[\frac{\Gamma(\omega_{j})}{\Gamma(\omega_{j}+\beta_{j})}\right]\right\\},$
(18) $\displaystyle\beta_{j}|\omega_{j},\mathbf{X}^{j}$ $\displaystyle\propto$
$\displaystyle\exp\left\\{\left[-\lambda_{2}+\sum_{k=1}^{m_{j}}{\log(1-X_{k}^{j})}\right]\beta_{j}-(\lambda_{0}+m_{j})\log\left[\frac{\Gamma(\beta_{j})}{\Gamma(\omega_{j}+\beta_{j})}\right]\right\\}.$
(19)

The detailed derivation for these posteriors can be found Appendix A.2.3.

Since the conditional posteriors in Equations (18) and (19) are not standard
distributions, we again develop an M-H nested Gibbs sampler to generate
samples for $\omega_{j}$ and $\beta_{j}$. Denote the samples from the previous
iteration in the M-H sampling by $\omega_{j}^{0}$ and $\beta_{j}^{0}$. By
using the Stirling approximation, we choose $\mbox{Gamma}(d,d/a)$ with
relatively small $d$ and mean $a$ equal to $\omega_{j}^{0}$ or $\beta_{j}^{0}$
as the proposal distribution; See the detailed derivation in Appendix A.2.3.
Denote the proposal density by $g(\cdot,a)$. Specifically, for $\omega_{j}$,
we randomly sample a candidate $\widetilde{\omega}_{j}$ from the proposal
distribution $\mbox{Gamma}(d,d/\omega_{j}^{0})$, and accept
$\widetilde{\omega}_{j}$ with probability

$\min\left\\{1,\frac{p(\widetilde{\omega}_{j}|\beta_{j}^{0},\mathbf{X}^{j})g(\omega_{j}^{0},\widetilde{\omega}_{j})}{p(\omega_{j}^{0}|\beta_{j}^{0},\mathbf{X}^{j})g(\widetilde{\omega}_{j},\omega_{j}^{0})}\right\\},$

where $p(\widetilde{\omega}_{j}|\beta_{j}^{0},\mathbf{X}^{j})$ and
$p(\omega_{j}^{0}|\beta_{j}^{0},\mathbf{X}^{j})$ are the conditional posterior
in Equation (18). Otherwise, retain the value of $\omega_{j}^{0}$. Similarly,
for $\beta_{j}$, we randomly sample a candidate $\widetilde{\beta}_{j}$ from
the proposal distribution $\mbox{Gamma}(d,d/\beta_{j}^{0})$, and accept
$\widetilde{\beta}_{j}$ with probability

$\min\left\\{1,\frac{p(\widetilde{\beta}_{j}|\omega_{j}^{0},\mathbf{X}^{j})g(\beta_{j}^{0},\widetilde{\beta}_{j})}{p(\beta_{j}^{0}|\omega_{j}^{0},\mathbf{X}^{j})g(\widetilde{\beta}_{j},\beta_{j}^{0})}\right\\},$

where $p(\widetilde{\beta}_{j}|\omega_{j}^{0},\mathbf{X}^{j})$ and
$p(\beta_{j}^{0}|\omega_{j}^{0},\mathbf{X}^{j})$ are the conditional
posteriors in Equation (19). Otherwise, retain the value of $\beta_{j}^{0}$.
In our empirical study, we set $d=2$ when we sample both $\omega_{j}$ and
$\beta_{j}$.

###  A.2 Derivation of the Results Used in the Gibbs Sampling

In this section, we provide the detailed derivation of priors, proposal
distributions, and conditional posteriors used in the Gibbs samplers for DPM
with Gamma, Gaussian and Beta kernel densities in Section A.1.

####  A.2.1 Conditional Posteriors of DPM with Gamma Kernel

We derive the conditional posteriors of parameters
$\boldsymbol{\psi}^{\star}_{j}=(V_{j},u_{j})$ with $j=1,2,\ldots,K_{0}$ for
DPM with Gamma kernel. Given the priors $V_{j}\sim\mbox{exp}(\theta)$,
$u_{j}\sim\mbox{Inv-Gamma}(r,s)$ and the likelihood
$X_{i}|c_{i}=j,\boldsymbol{\psi}^{\star}_{j}\sim\mbox{Gamma}(V_{j},V_{j}/u_{j})$,
by the Bayes’ rule, we have the conditional posterior for $V_{j}$

$\displaystyle p(V_{j}|\mathbf{X}^{j},u_{j})\propto
p(V_{j})\prod_{k=1}^{m_{j}}p\left(X^{j}_{k}|V_{j},u_{j}\right)$ (20)
$\displaystyle\propto$ $\displaystyle e^{-\theta
V_{j}}\prod_{k=1}^{m_{j}}\frac{(V_{j}/u_{j})^{V_{j}}}{\Gamma(V_{j})}(X_{k}^{j})^{V_{j}-1}e^{-(V_{j}/u_{j})X_{k}^{j}}$
$\displaystyle\propto$
$\displaystyle\frac{V_{j}^{m_{j}V_{j}}}{\Gamma(V_{j})^{m_{j}}}\exp\left\\{-V_{j}\left[\theta+\frac{\sum_{k=1}^{m_{j}}{X_{k}^{j}}}{u_{j}}+m_{j}\log(u_{j})-\sum_{k=1}^{m_{j}}{\log\left(X_{k}^{j}\right)}\right]\right\\}.$

Since the conditional posterior of $V_{j}$ in Equation (20) is not a standard
distribution, we develop an M-H sampling algorithm to generate samples of
$V_{j}$. We first find an appropriate proposal distribution for the M-H
sampling. To get a fair degree of probability drawing samples from the tail
part of the conditional posterior $p(V_{j}|\mathbf{X}^{j},u_{j})$, the
Stirling approximation, $n!\approx\sqrt{2\pi n}(n/\mbox{e})^{n}$ for large n,
is used to find an appropriate family for the proposal distribution. Since
$\Gamma(n)=(n-1)!$,

$\displaystyle
p\left(V_{j}|\mathbf{X}^{j},u_{j}\right)\propto\frac{V_{j}^{m_{j}V_{j}}}{\Gamma(V_{j})^{m_{j}}}\exp\left\\{-V_{j}\left[\theta+\frac{\sum_{k=1}^{m_{j}}{X_{k}^{j}}}{u_{j}}+m_{j}\log(u_{j})-\sum_{k=1}^{m_{j}}{\log(X_{k}^{j})}\right]\right\\}$
$\displaystyle\approx$
$\displaystyle\left[\frac{V_{j}^{V_{j}}}{\sqrt{2\pi(V_{j}-1)}{\left(\frac{V_{j}-1}{e}\right)^{V_{j}-1}}}\right]^{m_{j}}e^{-V_{j}B},\mbox{
if }V_{j}\mbox{ is large}$ $\displaystyle\approx$
$\displaystyle\left[\frac{V_{j}^{V_{j}}}{(V_{j}-1)^{V_{j}}}\sqrt{\frac{V_{j}-1}{2\pi}}e^{V_{j}-1}\right]^{m_{j}}e^{-V_{j}B}$
$\displaystyle\approx$
$\displaystyle\left(\sqrt{\frac{V_{j}-1}{2\pi}}e^{V_{j}}\right)^{m_{j}}e^{-V_{j}B}\approx\left(\frac{1}{2\pi}\right)^{m_{j}/2}(V_{j})^{m_{j}/2}e^{-V_{j}(B-m_{j})}$

where
$B=\theta+\sum_{k=1}^{m_{j}}{X_{k}^{j}}/u_{j}+m_{j}\log(u_{j})-\sum_{k=1}^{m_{j}}{\log(X_{k}^{j})}$.
This approximation holds when $V_{j}$ is large and it returns a Gamma kernel
function. Thus, we choose the proposal distribution to be
$\mbox{Gamma}(d,d/V_{j}^{0})$ with mean $V_{j}^{0}$ denoting the sample
obtained from the previous M-H iteration. To have a non-negligible probability
to draw samples far from $V_{j}^{0}$, the value of $d$ is recommended to be
small, e.g., $d=2$ used in our empirical study.

Next we derive the conditional posterior for parameter $u_{j}$. By applying
the Bayes’ rule, we have

$\displaystyle p(u_{j}|\mathbf{X}^{j},V_{j})\propto
p(u_{j})\prod_{k=1}^{m_{j}}p(X^{j}_{k}|V_{j},u_{j})$ $\displaystyle\propto$
$\displaystyle
u_{j}^{-(r+1)}e^{-s/u_{j}}\prod_{i=1}^{m_{j}}\frac{(V_{j}/u_{j})^{V_{j}}}{\Gamma(V_{j})}(X_{k}^{j})^{V_{j}-1}e^{-(V_{j}/u_{j})X_{k}^{j}}$
$\displaystyle\propto$ $\displaystyle
u_{j}^{-(r+1+m_{j}V_{j})}\exp\left[-\frac{s+V_{j}\sum_{k=1}^{m_{j}}X_{k}^{j}}{u_{j}}\right]$
$\displaystyle\sim$ $\displaystyle\mbox{Inv-
Gamma}\left(r+m_{j}V_{j},s+V_{j}\sum_{k=1}^{m_{j}}X_{k}^{j}\right).$

####  A.2.2 Conditional Posteriors of DPM with Gaussian Kernel

For DPM with Gaussian kernel, we choose a conditional conjugate joint prior
distribution for the component parameters
$\boldsymbol{\psi}^{\star}_{j}=(u_{j},\sigma^{2}_{j})$ with
$j=1,2,\ldots,K_{0}$,

$u_{j}|\sigma_{j}^{2}\sim\mathcal{N}(u_{0},\sigma_{j}^{2}/m_{0})\mbox{ and
}\sigma_{j}^{2}/\sigma_{0}^{2}\sim\mbox{Inv-
Gamma}\left(\frac{v_{0}}{2},\frac{1}{2}\right)$

which determines the base function $G_{0}(u,\sigma^{2})$ with hyper-parameters
$\boldsymbol{\theta}_{G}=(u_{0},m_{0},v_{0},\sigma_{0})$.

Here, we derive the conditional posteriors of the latent variables
$\mathbf{c}$. For $i=1,2,\ldots,m$, if $X_{i}$ is associated to an existing
component, by applying the Bayes’ rule,

$p(c_{i}=j|\mathbf{c}^{-i},\boldsymbol{\psi}^{\star}_{j},\alpha,X_{i})=b_{0}\cdot
p(c_{i}=j|\alpha,\mathbf{c}^{-i})p(X_{i}|c_{i}=j,\boldsymbol{\psi}^{\star}_{j})=b_{0}\frac{m_{j}^{-i}}{m+\alpha-1}\frac{1}{{\sqrt{2\pi}\sigma_{j}}}e^{-(X_{i}-u_{j})^{2}/2\sigma_{j}^{2}}.$

If $X_{i}$ is associated to a new component,

$\displaystyle
p(c_{i}=j|\mathbf{c}^{-i},\boldsymbol{\psi}^{\star}_{j},\alpha,X_{i})=b_{0}\cdot
p(c_{i}=j|\alpha,\mathbf{c}^{-i})p(X_{i}|c_{i}=j,\boldsymbol{\psi}^{\star}_{j})$
$\displaystyle=$ $\displaystyle
b_{0}\frac{\alpha}{m+\alpha-1}\int_{0}^{\infty}\int_{-\infty}^{\infty}p(X_{i}|u_{j},\sigma_{j}^{2})p(u_{j}|\sigma_{j}^{2})p(\sigma_{j}^{2})du_{j}d\sigma_{j}^{2}$
$\displaystyle=$ $\displaystyle
b_{0}\frac{\alpha}{m+\alpha-1}\int_{0}^{\infty}\int_{-\infty}^{\infty}(2\pi\sigma_{j}^{2})^{-1/2}e^{-\frac{(X_{i}-u_{j})^{2}}{2\sigma_{j}^{2}}}\times\left(\frac{2\pi\sigma_{j}^{2}}{m_{0}}\right)^{-1/2}\exp\left[-\frac{m_{0}(u_{j}-u_{0})^{2}}{2\sigma_{j}^{2}}\right]$
$\displaystyle\times\frac{(v_{0}/2)^{v_{0}/2}}{\Gamma(v_{0}/2)}\sigma_{0}^{v_{0}}(\sigma_{j}^{2})^{-(v_{0}/2+1)}e^{\left(-\frac{v_{0}\sigma_{0}^{2}}{2\sigma_{j}^{2}}\right)}du_{j}d\sigma_{j}^{2}$
$\displaystyle=$ $\displaystyle
b_{0}\frac{\alpha}{m+\alpha-1}\frac{(v_{0}/2)^{v_{0}/2}}{\Gamma(v_{0}/2)}\sigma_{0}^{v_{0}}\sqrt{m_{0}/2\pi}\int_{0}^{\infty}\int_{-\infty}^{\infty}(\sigma_{j}^{2})^{-(\frac{v_{0}+3}{2})}(2\pi\sigma_{j}^{2})^{-1/2}$
$\displaystyle\times\exp\left[-\left(\frac{(m_{0}+1)(u_{j}-\frac{X_{i}+m_{0}u_{0}}{m_{0}+1})^{2}+\frac{m_{0}(X_{i}-u_{0})^{2}}{m_{0}+1}+v_{0}\sigma_{0}^{2}}{2\sigma_{j}^{2}}\right)\right]du_{j}d\sigma_{j}^{2}$
$\displaystyle=$ $\displaystyle
b_{0}\frac{\alpha}{m+\alpha-1}\frac{(v_{0}/2)^{v_{0}/2}}{\Gamma(v_{0}/2)}\sigma_{0}^{v_{0}}\sqrt{\frac{m_{0}}{2\pi(m_{0}+1)}}\int_{0}^{\infty}\int_{-\infty}^{\infty}\left(\frac{2\pi\sigma_{j}^{2}}{m_{0}+1}\right)^{-1/2}$
$\displaystyle\times\exp\left[-\left(\frac{\left(u_{j}-\frac{X_{i}+m_{0}u_{0}}{m_{0}+1}\right)^{2}}{2\sigma_{j}^{2}/(m_{0}+1)}\right)\right]du_{j}\exp\left[-\left(\frac{\frac{m_{0}(X_{i}-\mu_{0})^{2}}{m_{0}+1}+v_{0}\sigma_{0}^{2}}{2\sigma_{j}^{2}}\right)\right](\sigma_{j}^{2})^{-(\frac{v_{0}+1}{2}+1)}d\sigma_{j}^{2}$
$\displaystyle=$ $\displaystyle
b_{0}\frac{\alpha}{m+\alpha-1}\frac{(v_{0}/2)^{v_{0}/2}}{\Gamma(v_{0}/2)}\sigma_{0}^{v_{0}}\sqrt{\frac{m_{0}}{2\pi(m_{0}+1)}}\int_{0}^{\infty}\exp\left[-\left(\frac{\frac{m_{0}(X_{i}-u_{0})^{2}}{m_{0}+1}+v_{0}\sigma_{0}^{2}}{2\sigma_{j}^{2}}\right)\right](\sigma_{j}^{2})^{-(\frac{v_{0}+1}{2}+1)}d\sigma_{j}^{2}$
$\displaystyle=$ $\displaystyle
b_{0}\frac{\alpha}{m+\alpha-1}\frac{(v_{0}/2)^{v_{0}/2}}{\Gamma(v_{0}/2)}\sigma_{0}^{v_{0}}\sqrt{\frac{m_{0}}{2\pi(m_{0}+1)}}\frac{\Gamma(A)}{B^{A}}$

where $b_{0}$ is a normalization constant, $A=\frac{v_{0}+1}{2}$ and
$B=\left[v_{0}\sigma_{0}^{2}+\frac{m_{0}(X_{i}-u_{0})^{2}}{m_{0}+1}\right]/2$.

####  A.2.3 Conditional Posteriors of DPM with Beta Kernel

In this section, we first find a conjugate joint prior and then derive the
conditional posteriors of parameters
$\boldsymbol{\psi}^{\star}_{j}=(\omega_{j},\beta_{j})$ with
$j=1,2,\ldots,K_{0}$ for DPM with Beta kernel density. The likelihood is
$X_{i}\mid
c_{i}=j,\boldsymbol{\psi}^{\star}_{j}\sim\mbox{Beta}(\omega_{j},\beta_{j})$.
Since Beta distribution belongs to the exponential family, we rewrite the Beta
density into the general form

$p(x|\omega_{j},\beta_{j})=\frac{\Gamma(\omega_{j}+\beta_{j})}{\Gamma(\omega_{j})\Gamma(\beta_{j})}x^{\omega_{j}-1}(1-x)^{\beta_{j}-1}=\frac{\Gamma(\omega_{j}+\beta_{j})}{\Gamma(\omega_{j})\Gamma(\beta_{j})}e^{(\omega_{j}-1)\log(x)+(\beta_{j}-1)\log(1-x)}.$

Thus, we choose a conjugate joint prior for $(\omega_{j},\beta_{j})$ with the
hyper-parameters
$\boldsymbol{\theta}_{G}=(\lambda_{0},\lambda_{1},\lambda_{2})$ (Chick, 2001)

$\omega_{j},\beta_{j}|\lambda_{0},\lambda_{1},\lambda_{2}\propto\exp\left\\{-\lambda_{1}\omega_{j}-\lambda_{2}\beta_{j}-\lambda_{0}\log\left[\frac{\Gamma(\omega_{j})\Gamma(\beta_{j})}{\Gamma(\omega_{j}+\beta_{j})}\right]\right\\}.$

Then, we derive the conditional posteriors for parameters
$(\omega_{j},\beta_{j})$ used in the Gibbs sampler in Appendix A.1.3. By
applying the Bayes’ rule, the conditional posterior for $\omega_{j}$ is

$\displaystyle p(\omega_{j}|\beta_{j},\mathbf{X}^{j})\propto
p(\omega_{j}|\beta_{j})p(\mathbf{X}^{j}|\omega_{j},\beta_{j})$ (21)
$\displaystyle\propto$
$\displaystyle\exp\left\\{-\lambda_{1}\omega_{j}-\lambda_{0}\log\left[\frac{\Gamma(\omega_{j})}{\Gamma(\omega_{j}+\beta_{j})}\right]\right\\}\prod_{k=1}^{m_{j}}\frac{\Gamma(\omega_{j}+\beta_{j})}{\Gamma(\omega_{j})}(X_{k}^{j})^{\omega_{j}-1}$
$\displaystyle\propto$
$\displaystyle\exp\left\\{\left(-\lambda_{1}+\sum_{k=1}^{m_{j}}{\log(X_{k}^{j})}\right)\omega_{j}-(\lambda_{0}+m_{j})\log\left[\frac{\Gamma(\omega_{j})}{\Gamma(\omega_{j}+\beta_{j})}\right]\right\\}.$

Since the conditional posterior for $\omega_{j}$ in Equation (21) is not a
standard distribution, we develop an M-H sampling algorithm to draw samples of
$\omega_{j}$ by following the similar procedure used in DPM with Gamma kernel
density. The Stirling approximation is used to find an appropriate proposal
distribution family. As $\omega_{j}$ is large, the conditional posterior
distribution can be approximated by

$\displaystyle p(\omega_{j}|\beta_{j},\mathbf{X}^{j})\propto
e^{\left(-\lambda_{1}+\sum_{k=1}^{m_{j}}{\log(X_{k}^{j})}\right)\omega_{j}-(\lambda_{0}+m_{j})\log\left[\frac{\Gamma(\omega_{j})}{\Gamma(\omega_{j}+\beta_{j})}\right]}$
$\displaystyle\approx$ $\displaystyle
e^{\left(-\lambda_{1}+\sum_{k=1}^{m_{j}}{\log(X_{k}^{j})}\right)\omega_{j}}\left[\frac{(\omega_{j}+\beta_{j}-1)!}{(\omega_{j}-1)!}\right]^{\lambda_{0}+m_{j}}$
$\displaystyle\approx$ $\displaystyle
e^{-\left(\lambda_{1}-\sum_{k=1}^{m_{j}}{\log(X_{k}^{j})}\right)\omega_{j}}\left(\omega_{j}^{\beta_{j}}\right)^{\lambda_{0}+m_{j}},\mbox{
if }\omega_{j}\mbox{ is large}$ $\displaystyle\sim$
$\displaystyle\mbox{Gamma}\left(\beta_{j}(\lambda_{0}+m_{j})+1,\lambda_{1}-\sum_{k=1}^{m_{j}}{\log(X_{k}^{j})}\right).$

Thus, $\mbox{Gamma}(d,d/\omega_{j}^{0})$ with small $d$, e.g., $d=2$ used in
the empirical study, is used as the proposal distribution, where
$\omega_{j}^{0}$ denotes the sample obtained from the previous M-H sampling
iteration.

Next, by applying the Bayes’ rule, we derive the conditional posterior for
$\beta_{j}$

$\displaystyle p(\beta_{j}|\omega_{j},\mathbf{X}^{j})\propto
p(\beta_{j}|\omega_{j})p(\mathbf{X}^{j}|\omega_{j},\beta_{j})$ (22)
$\displaystyle\propto$
$\displaystyle\exp\left\\{-\lambda_{2}\beta_{j}-\lambda_{0}\log\left[\frac{\Gamma(\beta_{j})}{\Gamma(\omega_{j}+\beta_{j})}\right]\right\\}\prod_{k=1}^{m_{j}}\frac{\Gamma(\omega_{j}+\beta_{j})}{\Gamma(\beta_{j})}\left(1-X_{k}^{j}\right)^{\beta_{j}-1}$
$\displaystyle\propto$
$\displaystyle\exp\left\\{\left(-\lambda_{2}+\sum_{k=1}^{m_{j}}{\log(1-X_{k}^{j})}\right)\beta_{j}-(\lambda_{0}+m_{j})\log\left[\frac{\Gamma(\beta_{j})}{\Gamma(\omega_{j}+\beta_{j})}\right]\right\\}.$

Notice that Equations (21) and (22) have the similar form, and they do not
belong to any standard distribution. Thus, an M-H sampling approach is
developed to generate samples for $\beta_{j}$. An appropriate proposal
distribution family is found by applying the Stirling approximation,

$\displaystyle p(\beta_{j}|\omega_{j},\mathbf{X}^{j})\propto
e^{\left(-\lambda_{2}+\sum_{k=1}^{m_{j}}{\log(1-X_{k}^{j})}\right)\beta_{j}-(\lambda_{0}+m_{j})\log\left[\frac{\Gamma(\beta_{j})}{\Gamma(\omega_{j}+\beta_{j})}\right]}$
$\displaystyle\approx$ $\displaystyle
e^{\left(-\lambda_{2}+\sum_{k=1}^{m_{j}}{\log(1-X_{k}^{j})}\right)\beta_{j}}\left[\frac{(\omega_{j}+\beta_{j}-1)!}{(\beta_{j}-1)!}\right]^{\lambda_{0}+m_{j}}$
$\displaystyle\approx$ $\displaystyle
e^{-\left(\lambda_{2}-\sum_{k=1}^{m_{j}}{\log(1-X_{k}^{j})}\right)\beta_{j}}\left(\beta_{j}^{\omega_{j}}\right)^{\lambda_{0}+m_{j}},\mbox{
if }\beta_{j}\mbox{ is large}$ $\displaystyle\sim$
$\displaystyle\mbox{Gamma}\left(\omega_{j}(\lambda_{0}+m_{j})+1,\lambda_{2}-\sum_{k=1}^{m_{j}}{\log(1-X_{k}^{j})}\right).$

In the M-H sampling, $\mbox{Gamma}(d,d/\beta_{j}^{0})$ with small $d$ is used
as the proposal distribution, where $\beta_{j}^{0}$ denotes the sample
obtained from the previous iteration.

##  Appendix B Posterior Consistency of DPM Nonparametric Input Models

We first introduce a series of basic definitions and theorems related to
posterior consistency from Bayesian nonparametrics theory. We refer the
readers to the textbooks of Ghosh and Ramamoorthi (2003) and Ghosal and van
der Vaart (2017) for technical details.

######  Definition 1.

(Ghosal et al. 1999b) Let $\mathscr{F}$ be the set of all densities on $\Re$
with respect to the Lebesgue measure on $\Re$. Let $f^{c}\in\mathscr{F}$
denote the true probability density and $P_{f^{c}}$ be its associated
probability measure. A weak neighborhood $U$ of $f^{c}$ is a set containing a
set of the form

$V=\left\\{f\in\mathscr{F}:~{}\left|\int\phi_{i}(x)f(x)dx-\int\phi_{i}(x)f^{c}(x)dx\right|<\epsilon,~{}i=1,\ldots,k\right\\},$

where $\phi_{i}$’s are bounded continuous functions on $\Re$ and $k$ is a
positive integer.

The weak neighborhood is defined on the space of probability measures
topologized by weak convergence (convergence in distribution). Refer to
Billingsley (1999) Chapter 1 Section 2 for more details on weak convergence.

######  Definition 2.

(Ghosal et al. 1999b) Let $\mathbf{X}_{m}=\\{X_{1},\ldots,X_{m}\\}$ be an
i.i.d. sample from $F^{c}$ (with density $f^{c}$). The posterior distribution
$p(\cdot\mid\mathbf{X}_{m})$ is said to be weekly consistent at $F^{c}$ or
$f^{c}$, if with $P_{f^{c}}$-probability 1,

$p(U\mid\mathbf{X}_{m})\to 1,\quad\text{as }m\to\infty,$ (23)

for all weak neighborhoods $U$ of $f^{c}$.

Then we define the concept of Kullback-Leibler (K-L) support:

######  Definition 3.

(Ghosal et al. 1999b) Let $p$ be a prior distribution over the space
$\mathscr{F}$, the set of all densities on $\Re$ with respect to the Lebesgue
measure on $\Re$. A density $f^{c}$ is said to be in the Kullback-Leibler
(K-L) support of the prior $p$ (denoted by $f^{c}\in KL(p)$), if for all
$\epsilon>0$, $p(\mathcal{K}_{\epsilon}(f^{c}))>0$, where
$\mathcal{K}_{\epsilon}(f^{c})=\\{g\in\mathscr{F}:\int
f^{c}(x)\log\frac{f^{c}(x)}{g(x)}dx<\epsilon\\}$ is the K-L neighborhood of
$f^{c}$.

Following this definition, we cite Theorem 4.4.2 in Ghosh and Ramamoorthi
(2003), which is essentially derived from the Schwartz theorem (Schwartz
1965):

######  Theorem 4.

(Ghosh and Ramamoorthi 2003 Theorem 4.4.2) Let $p$ be a prior distribution
over the space $\mathscr{F}$, the set of all densities on $\Re$ with respect
to the Lebesgue measure on $\Re$. If $f^{c}$ is in the Kullback-Leibler (K-L)
support of $p$, then the posterior is weakly consistent at $f^{c}$.

Theorem 4 shows that $f^{c}$ being in the K-L support of the prior $p$ implies
the posterior weak consistency. Therefore, to prove Theorem 1, it is
sufficient to show that $f^{c}$ is in the K-L support of those priors under
the conditions of Theorem 1.

In this paper, we define the support of a probability measure using the
following standard definition on page 23 of Billingsley (1999):

######  Definition 4.

(Billingsley 1999) If $\mathcal{F}$ is a $\sigma$-field in $\Omega$ and $P$ is
a probability measure on $\mathcal{F}$, the triple $(\Omega,\mathcal{F},P)$ is
called a probability measure space, or simply a probability space. A support
of $P$ is any $\mathcal{F}$-measurable set $A$ for which $P(A)=1$, denoted by
$A=supp(P)$.

To make our proof self-contained, we cite the original theorems and lemmas
given in Wu and Ghosal (2008) that are used in the proof of Theorem 1. In the
citation, we changed the notation system in the original paper to the one used
in this paper for easy understanding.

######  Theorem 5.

(Theorem 1 of Wu and Ghosal 2008) Let $f^{c}$ be the true density. Let
$h(x;\psi,\phi)$ be a kernel density, where $\psi$ is the mixing parameter and
$\phi$ is the hyper-parameter which lies in the parameter space $\Phi$.
Consider the mixture distribution $f_{G,\phi}=\int h(x;\psi,\phi)dG(\psi)$,
where $G$ lies in the space of all mixing distributions $\mathscr{G}$. Let
$\mu$ and $\Pi$ be priors for the hyper parameter and the mixing distribution,
and $p$ be the prior on $\mathscr{F}$ induced by $\mu\times\Pi$. If for any
$\epsilon>0$, there exists a mixing distribution $G_{\epsilon}$, a hyper-
parameter $\phi_{\epsilon}$, a set $A\subset\Phi$ with $\mu(A)>0$, and a set
$\mathscr{W}\subset\mathscr{G}$ with $\Pi(\mathscr{W})>0$, such that

  * A1.

$\int f^{c}\log\frac{f^{c}}{f_{G_{\epsilon},\phi_{\epsilon}}}<\epsilon$,

  * A2.

$\int
f^{c}\log\frac{f_{G_{\epsilon},\phi_{\epsilon}}}{f_{G_{\epsilon},\phi}}<\epsilon$
for every $\phi\in A$, and

  * A3.

$\int f^{c}\log\frac{f_{G_{\epsilon},\phi}}{f_{G,\phi}}<\epsilon$ for every
$G\in\mathscr{W}$, $\phi\in A$,

then $f^{c}\in KL(p)$.

This is Theorem 1 in Wu and Ghosal (2008), which was constructed for Bayesian
nonparametric kernel mixture models that are more general than the model (1).
The kernel function for DPM considered in the cited theorem contains two of
parameters $\psi$ and $\phi$, while the prior being focused in this paper as
presented in (1) has only one parameter $\psi$. As same as presented in (1),
$\phi$ in Theorem 5 is mixed over by the mixing distribution $G$, which is
further given a prior $DP(\alpha,G_{0})$. The additional parameter $\phi$ in
Theorem 5 is known as the hyper-parameter and is given a separate prior $\mu$
directly. We use $\Psi$ to denote the sample space of $\psi$, and
$\mathscr{M}(\Psi)$ to denote the space of mixing distributions on $\Psi$.
Notice that $\mathscr{W}\subset\mathscr{M}(\Psi)$.

Prior (1) does not involve the hyper-parameter $\phi$, and hence it is
equivalent to the prior induced by $\Pi$ only, instead of by $\mu\times\Pi$.
Specifically, $\Pi$ is the Dirichlet Process (DP) with parameters $\alpha$ and
base measure $G_{0}$ on $\mathscr{G}$, the space of mixing distribution $G$;
and $p$ is the induced Dirichlet Process Mixtrure (DPM) on $\mathscr{F}$, the
space of density functions. Therefore, to prove Theorem 1, we only need to
verify Conditions A1 and A3 in Theorem 5.

Applying Theorem 5 to prove the consistency is to verify that Conditions A1,
A2 and A3 are satisfied for a given prior. Usually, Condition A1 is directly
verified by construction, while A2 and A3 are verified through the following
two lemmas.

######  Lemma 1.

(Lemma 2 of Wu and Ghosal 2008) Let $f^{c}$, $\Pi$, $\mu$ and $p$ be the same
as in Theorem 5. If for any $\epsilon>0$, there exist a mixing distribution
$G_{\epsilon}$, a set $D\supseteq supp(G_{\epsilon})$, and $\phi_{\epsilon}\in
supp(\mu)$ such that Condition A1 holds and the kernel density function $h$
satisfies

  * A4.

for any given $x$ and $\psi$, the map $\phi\mapsto h(x;\psi,\phi)$ is
continuous on the interior of $supp(\mu)$;

  * A5.

$\int_{\mathfrak{X}}f^{c}(x)\left\\{\left|\log\frac{\sup_{\psi\in
D}h(x;\psi,\phi_{\epsilon})}{\inf_{\psi\in
D}h(x;\psi,\phi)}\right|+\left|\log\frac{\sup_{\psi\in
D}h(x;\psi,\phi)}{\inf_{\psi\in
D}h(x;\psi,\phi_{\epsilon})}\right|\right\\}dx<\infty$ for every $\phi\in
N(\phi_{\epsilon})$, where $N(\phi_{\epsilon})$ is an open neighborhood of
$\phi_{\epsilon}$;

  * A6.

for any given $x\in\mathfrak{X}$, $\psi\in D$ and $\phi\in
N(\phi_{\epsilon})$, there exists a function $k(x,\psi)$ such that
$k(x,\psi)\geq h(x;\psi,\phi)$, and $\int
k(x,\psi)dG_{\epsilon}(\psi)<\infty$;

then there exists a set $A\subset\Phi$ such that Condition A2 holds.

######  Lemma 2.

(Lemma 3 of Wu and Ghosal 2008) Let $f^{c}$, $\Pi$, $\mu$ and $p$ be the same
as in Theorem 5. If for any $\epsilon>0$, there exist a mixing distribution
$G_{\epsilon}\in supp(\Pi)$, a hyper-parameter $\phi_{\epsilon}\in supp(\mu)$,
and a set $A\in\Phi$ with $\mu(A)>0$, such that Conditions A1 and A2 hold and
for some $D\supseteq supp(G_{\epsilon})$, the kernel density function $h$ and
the prior $\Pi$ satisfy

  * A7.

for any $\phi\in A$,
$\int_{\mathfrak{X}}f^{c}(x)\log\frac{f_{G_{\epsilon},\phi}(x)}{\inf_{\psi\in
D}h(x;\psi,\phi)}dx<\infty$;

  * A8.

for any $\epsilon>0$, there exists a compact set $C\subset\mathfrak{X}$ with
the complement set denoted by $C^{c}$, such that

$\int_{C^{c}}f^{c}(x)\log\frac{f_{G_{\epsilon},\phi}(x)}{\inf_{\psi\in
D}h(x;\psi,\phi)}dx<\epsilon/4,$ (24)

and $P_{f^{c}}(C^{c})<\epsilon/(4\log 2)$, we have that $c:=\inf_{x\in
C}\inf_{\psi\in D}h(x;\psi,\phi)>0$;

  * A9.

for any given $\phi\in A$ and some compact $C\subset\mathfrak{X}$ as required
in A8, such that the family of maps $\\{\psi\mapsto h(x;\psi,\phi),x\in C\\}$
is uniformly equicontinuous on $D$;

then there exists a set $\mathscr{W}\subset\mathscr{M}(\Psi)$ such that
Condition A3 holds and $\Pi(\mathscr{W})>0$.

Notice that, Conditions A8 and A9 are a little different from their original
form in Wu and Ghosal (2008). The modification of Condition A9 follows Wu et
al. (2009), while the modification on A8 is justified by the fact that in
detailed proof of this lemma in Wu and Ghosal (2008). We only need the
existence of $C$ that satisfies (24) and $P_{f^{c}}(C^{c})<\epsilon/(4\log
2)$, where the existence is implied by Condition A7. Since the conditions of
the lemma are different from the original version in Wu and Ghosal (2008), we
include the proof of this lemma below, which reflexes the change in the
conditions.

Proof of Lemma 2:

For any $\phi\in A$, write

$\displaystyle\int_{\mathfrak{X}}f^{c}(x)\log\frac{f_{G_{\epsilon},\phi}(x)}{f_{G,\phi}(x)}dx$
$\displaystyle=$
$\displaystyle\int_{C^{c}}f^{c}(x)\log\frac{f_{G_{\epsilon},\phi}(x)}{f_{G,\phi}(x)}dx$
(25)
$\displaystyle+\int_{C}f^{c}(x)\log\frac{f_{G_{\epsilon},\phi}(x)}{f_{G,\phi}(x)}dx.$

Now, since $G_{\epsilon}(D)=1>\frac{1}{2}$, $\mathscr{V}=\\{G:\
G({D})>\frac{1}{2}\\}$ is an open neighborhood of $G_{\epsilon}$ by the
Portmanteau Theorem. For any $G\in\mathscr{V}$ and $\phi\in A$,

$\displaystyle\int_{C^{c}}f^{c}(x)\log\frac{f_{G_{\epsilon},\phi}(x)}{f_{G,\phi}(x)}dx$
$\displaystyle\leq\int_{C^{c}}f^{c}(x)\log\frac{f_{G_{\epsilon},\phi}(x)}{\int_{\psi\in
D}\inf_{\psi\in D}h(x;\psi,\phi)dG(\psi)}dx$
$\displaystyle\leq\int_{C^{c}}f^{c}(x)\log\frac{f_{G_{\epsilon},\phi}(x)}{\inf_{\psi\in
D}h(x;\psi,\phi)\int_{\psi\in D}dG(\psi)}dx$
$\displaystyle<{\int_{C^{c}}f^{c}(x)\log\frac{2f_{G_{\epsilon},\phi}(x)}{\inf_{\psi\in
D}h(x;\psi,\phi)}dx}$

The last term in the above inequality is due to the fact that $G(D)>1/2$ as
defined in the definition of $\mathscr{V}$. Now we have that

$\displaystyle\int_{C^{c}}f^{c}(x)\log\frac{2f_{G_{\epsilon},\phi}(x)}{\inf_{\psi\in
D}h(x;\psi,\phi)}dx$
$\displaystyle=\int_{C^{c}}f^{c}(x)\log\frac{f_{G_{\epsilon},\phi}(x)}{\inf_{\psi\in
D}h(x;\psi,\phi)}dx+(\log 2)P_{f^{c}}(C^{c});$

where $P_{f^{c}}(C^{c})=\int_{C^{c}}f^{c}(x)dx$, and $P_{f^{c}}$ denotes the
probability measure corresponding to $f^{c}$. By Condition A7, there exists
compact $C\subset\mathfrak{X}$, such that

$\int_{C^{c}}f^{c}(x)\log\frac{f_{G_{\epsilon},\phi}(x)}{\inf_{\psi\in
D}h(x;\psi,\phi)}dx<\epsilon/4.$ (26)

We can further ensure that $P_{f^{c}}(C^{c})<\epsilon/(4\log 2)$, so the bound
for $\int_{C^{c}}f^{c}\log\frac{f_{G_{\epsilon,\phi}}}{f_{G,\phi}}$ is less
than $\epsilon/2$. Now, if we can show that for the given $\epsilon>0$, there
exists a weak neighborhood $\mathscr{U}$ of $G_{\epsilon}$, such that
$\int_{C}f^{c}(x)\log\frac{f_{G_{\epsilon},\phi}(x)}{f_{G,\phi}(x)}dx<\epsilon/2$
for any $G\in\mathscr{U}$ and $\phi\in A$, then Lemma 2 is proved by letting
$\mathscr{W}=\mathscr{U}\cap\mathscr{V}$.

Observing that for any given $\phi\in A$, the family of maps $\\{\psi\mapsto
h(x;\psi,\phi):x\in C\\}$ is uniformly equicontinuous on $D\subset\Psi$. By
the Arzela-Ascoli theorem (see Royden (1988) [pp. 169]), for any $\delta>0$,
there exist $x_{1},x_{2},\ldots,x_{m}$, such that, for any $x\in C$,

$\sup_{\psi\in D}|h(x;\psi,\phi)-h(x_{i};\psi,\phi)|<c\delta.$ (27)

for some $i=1,2,\ldots,k$.

Let
$\mathscr{U}=\\{G:|\int_{D}h(x_{i};\psi,\phi)dG_{\epsilon}(\psi)-\int_{D}h(x;\psi,\phi)dG(\psi)|<c\delta,\
\ x\in C,i=1,2,\ldots,k\\}.$ Then $\mathscr{U}$ is an open weak neighborhoods
of $G_{\epsilon}$ since $G_{\epsilon}\in\mathrm{supp}(\Pi)$ and
$G_{\epsilon}(\partial D)=0$. Hence, given $\phi\in A$, for any
$G\in\mathscr{U}\cap\mathscr{V}$ and $x\in C$, by applying Condition A8,

$\displaystyle\left|\frac{\int_{\Psi}h(x;\psi,\phi)dG_{\epsilon}(\psi)}{\int_{\Psi}h(x;\psi,\phi)dG(\psi)}-1\right|$
$\displaystyle\leq$
$\displaystyle\left|\frac{\int_{D}h(x;\psi,\phi)dG_{\epsilon}(\psi)}{\int_{D}h(x;\psi,\phi)dG(\psi)}-1\right|$
$\displaystyle=$
$\displaystyle\left|\frac{\int_{D}h(x;\psi,\phi)dG_{\epsilon}(\psi)-\int_{D}h(x;\psi,\phi)dG(\psi)}{\int_{D}h(x;\psi,\phi)dG(\psi)}\right|$
$\displaystyle{\leq}$ $\displaystyle\frac{c\delta}{c/2}=2\delta,$

since
$|\int_{D}h(x;\psi,\phi)dG_{\epsilon}(\psi)-\int_{D}h(x;\psi,\phi)dG(\psi)|<c\delta$
for any $G\in\mathscr{U}$, and $\int_{D}h(x;\psi,\phi)dG(\psi)\leq\inf_{x\in
C}\inf_{\psi\in D}h(x;\psi,\phi)G(D)\geq c/2$. By choosing $\delta$ small
enough, we can ensure that the right hand side (RHS) of the last display is
less than $\epsilon/2$. Hence, for any given $\phi\in A$

$\displaystyle\int_{C}f^{c}(x)\log\frac{f_{G_{\epsilon},\phi}(x)}{f_{G,\phi}(x)}dx$
$\displaystyle\leq\sup_{x\in
C}\log\frac{f_{G_{\epsilon},\phi}(x)}{f_{G,\phi}(x)}\leq\sup_{x\in
C}\left|\frac{f_{G_{\epsilon},\phi}(x)}{f_{G,\phi}(x)}-1\right|$
$\displaystyle\leq\sup_{x\in
C}\left|\frac{\int_{\Psi}h(x;\psi,\phi)dG_{\epsilon}(\psi)}{\int_{\Psi}h(x;\psi,\phi)dG(\psi)}-1\right|<\epsilon/2$

for any $G\in\mathscr{U}\cap\mathscr{V}$. $\mbox{}\Box$

Next, we present several lemmas which will be used for proving Theorem 1 Part
(i), together with their proofs when necessary. Differing with Theorem 14 in
Wu and Ghosal (2008), the restriction of $f^{c}(0)>0$ is not necessary here.
In the proofs of the Lemmas and Theorem 1 (i), we will provide the detailed
explanations.

######  Lemma 3.

(Lemma 7 in Wu and Ghosal (2008)) Let $l$ be a positive integer and
$A_{l}=[a_{l},b_{l}]\subset\mathfrak{X}$, and let $h_{l}(x,t)$ be a sequence
of continuous functions for $x\in\mathfrak{X}$ and $t\in A_{l}$. Define
$f_{l}(x)=\int_{A_{l}}h_{l}(x,t)f(t)dt$, $m=1,2,\ldots$, where $f$ is bounded,
uniformly continuous and integrable on $\mathfrak{X}$. If $h_{l}$ satisfies
that for any $a_{l},b_{l}\in\mathfrak{X}$,

  * C1.

$\int_{A_{l}}h_{l}(x,t)dt\rightarrow 1$ as $l\rightarrow\infty$,

  * C2.

for each $\delta>0$, $\int_{|x-t|>\delta,t\in A_{l}}|h_{l}(x,t)|dt\rightarrow
0$ as $l\rightarrow\infty$,

  * C3.

$\int_{A_{l}}|h_{l}(x,t)|dt\leq M(x)<\infty$ for each $x\in\mathfrak{X}$,
$l=1,2\ldots$, where the bound $M(x)$ may depend on $x$ but not on $l$,

then $f_{l}(x)\rightarrow f(x)$ for each $x\in\mathfrak{X}$ as $l\to\infty$.

Proof of Lemma 3:

For any $\epsilon>0$, there is $\delta>0$ so small that that
$|f(t)-f(x)|<\epsilon$ for $|x-t|\leq\delta$. By Condition C1,

$f_{l}(x)-f(x)=\int_{A_{l}}[f(t)-f(x)]h_{l}(x,t)dt+o(1),$ (28)

where the last term goes to 0 as $l\rightarrow\infty$, for each
$x\in\mathfrak{X}$. To complete the proof, we will show that
$\int_{A_{l}}[f(t)-f(x)]h_{l}(x,t)dt\to 0$ as $l\to\infty$ for all
$x\in\mathfrak{X}$. By Condition C3, we have that

$\Big{|}\int_{|x-t|\leq\delta,\,t\in
A_{l}}[f(t)-f(x)]h_{l}(x,t)dt\Big{|}\leq\epsilon\int_{|x-t|\leq\delta,\,t\in
A_{l}}|h_{l}(x,t)|dt\leq\epsilon M(x).$ (29)

We also have that

$\int_{|x-t|>\delta,\,t\in A_{l}}[f(t)-f(x)]h_{l}(x,t)dt\rightarrow 0\mbox{ as
}l\rightarrow\infty,$ (30)

by Condition C2, and the condition that function $f$ is bounded. Combining
(28)-(30), we have $|f_{l}(x)-f(x)|\leq\epsilon M(x)+o(1)$, and hence the
result follows. $\Box$

In the following, we first fix some notation for proving the posterior
consistency of DPM of Gamma. Let $h_{l}(x;V)\equiv h(x;V,V/l)$ be the Gamma
density kernel with shape parameter $V$ and mean equals to $V/l$, where $V>0$
and $l>0$. Define

$h_{l}(x;V)=h(x;V,V/l)=\frac{l^{V}}{\Gamma(V)}x^{V-1}e^{-x/l}$ (31)

Let

$f_{l}(x)=t_{l}\int_{2}^{1+l^{2}}h_{l}(x;V)l^{-1}f^{c}((V-1)/l)dV$ (32)

where $t_{l}=(\int_{l^{-1}}^{l}f^{c}(s)ds)^{-1}$.

######  Lemma 4.

(Lemma 8 in Wu and Ghosal (2008)) Let $h_{l}(x;V)$ be defined as (31). If
Condition (d) in Theorem 1 is satisfied, then there exists a function
$0<C(x)<1$ such that for all sufficiently large $l>0$,

$C(x)\leq\left\\{\begin{array}[]{ll}\int_{l^{-1}\vee
x}^{x+\delta}h_{l}(x;lv+1)dv,&\ l^{-1}<x<1,\\\ \mbox{}&\mbox{}\\\
\int_{x-\delta}^{l\wedge x}h_{l}(x;lv+1)dv,&\ 1\leq x\leq
l+l^{-1},\end{array}\right.$ (33)

and $\int_{\mathfrak{X}}f^{c}(x)\log\frac{1}{C(x)}dx<\infty$.

Proof of Lemma 4:

For $l^{-1}<x<1$, applying Stirling’s inequality and noting that
$v<x+\delta<1+\delta$ (with $\delta>0$ chosen later) in the following
integral, it follows that

$\displaystyle\int_{{l}^{-1}\vee x}^{x+\delta}h_{l}(x;lv+1)dv$ (34)
$\displaystyle=$ $\displaystyle\int_{l^{-1}\vee
x}^{x+\delta}\frac{l^{lv+1}x^{lv}e^{-lx}}{\Gamma(lv+1)}dv$ $\displaystyle\geq$
$\displaystyle\int_{l^{-1}\vee
x}^{x+\delta}\frac{l^{lv+1}x^{lv}e^{-lx}}{\sqrt{2\pi}(lv+1)^{lv+1/2}\exp\\{-(lv+1)+(12x)^{-1}\\}}dv$
$\displaystyle=$
$\displaystyle\sqrt{\frac{l}{2\pi}}\exp(1-(12x)^{-1})\int_{l^{-1}\vee
x}^{x+\delta}\frac{x^{lv}e^{l(v-x)}}{(v+l^{-1})^{lv+{1}/{2}}}dv$
$\displaystyle\geq$
$\displaystyle\frac{\sqrt{l}}{\sqrt{2\pi(1+\delta+l^{-1})}}\exp(1-(12x)^{-1})\int_{l^{-1}\vee