11footnotetext: Min-ge Xie is a Distinguished Professor, Department of
Statistics, Rutgers, The State University of New Jersey, Piscataway, NJ 08854,
Email<EMAIL_ADDRESS>Peng Wang is an Associate Professor, Department
of Operations, Business Analytics and Information System, Email:
<EMAIL_ADDRESS>The research is supported in part by NSF grants
DMS1812048, DMS2015373 and DMS2027855.

# Repro Samples Method for Finite- and Large-Sample Inferences

Min-ge Xie and Peng Wang

###### Abstract

This article presents a novel, general, and effective simulation-inspired
approach, called repro samples method, to conduct statistical inference. The
approach studies the performance of artificial samples, referred to as repro
samples, obtained by mimicking the true observed sample to achieve uncertainty
quantification and construct confidence sets for parameters of interest with
guaranteed coverage rates. Both exact and asymptotic inferences are developed.
An attractive feature of the general framework developed is that it does not
rely on the large sample central limit theorem and is likelihood-free. As
such, it is thus effective for complicated inference problems which we can not
solve using the large sample central limit theorem. The proposed method is
applicable to a wide range of problems, including many open questions where
solutions were previously unavailable, for example, those involving discrete
or non-numerical parameters. To reduce the large computational cost of such
inference problems, we develop a unique matching scheme to obtain a data-
driven candidate set. Moreover, we show the advantages of the proposed
framework over the classical Neyman-Pearson framework. We demonstrate the
effectiveness of the proposed approach on various models throughout the paper
and provide a case study that addresses an open inference question on how to
quantify the uncertainty for the unknown number of components in a normal
mixture model. To evaluate the empirical performance of our repro samples
method, we conduct simulations and study real data examples with comparisons
to existing approaches. Although the development pertains to the settings
where the large sample central limit theorem does not apply, it also has
direct extensions to the cases where the central limit theorem does hold.

Key words: Artificial samples; Conformal or matching samples; Discrete and
continuous parameters; Parametric and nonparametric settings; Mixture model;
Nuisance parameters

##  1 Introduction

Statistical inference includes quantifying uncertainty inherited from sampling
based on an observed copy of data. The approach of repeatedly creating similar
artificial (Monte-Carlo) data to help assess the uncertainty has proved in the
literature to be an effective method. Examples include bootstrap, approximate
Bayesian computing (ABC), generalized fiducial inference (GFI) algorithm,
permutation, among others; see, e.g., Efron and Tibshirani (1993); Fishman
(1996); Robert (2016); Hannig et al. (2016), and references therein. Most of
these works, however, rely on the large sample central limit theorem (CLT) or
approximations to justify their validity of statistical inference. Their
applicability to some complicated inference problems, especially those to
which the large sample CLT does not apply, is limited. In this article, we
develop a novel and wide reaching artificial-sample-based inferential
framework, which we support with both finite and large sample theories. The
proposed framework has the advantage of providing valid inference without
relying on likelihood functions or large sample theories. It is especially
effective for difficult inference problems such as those involving discrete
and/or non-numerical parameters (e.g., the unknown number of clusters, model
A, B, C, etc.) among others. Although the repro samples method pertains to
settings in which the large sample CLT does not apply, it also has direct
extensions to the cases where the large sample theorem does hold.

Suppose that the sample data $\bm{Y}=(Y_{1},\ldots,Y_{n})^{\top}\in{\cal
Y}\subset\mathbb{R}^{n}$ are generated from a model
$\bm{Y}|{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits}_{0}\sim
M_{{\bf\theta}_{0}}(\cdot),$ where
${\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits}_{0}\in\Theta$ is the
true value of the model parameters that can be either numerical or
nonnumerical or a mixture of both. We assume (and assume only) that we know
how to simulate data from $M_{\bf\theta}$, given
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$. Often,
$\bm{Y}|{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits}_{0}\sim
M_{{\bf\theta}_{0}}(\cdot)$ can be re-expressed in a structure equation form
(also referred to as an “algorithmic model”):

$\bm{Y}=G({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits}_{0},\bm{U}),$
(1)

where $G(\cdot,\cdot)$ is a known mapping from $\Theta\times{\cal
U}\mapsto{\cal Y}$ and $\bm{U}=(U_{1},\ldots U_{r})^{\top}\in{\cal
U}\subset\mathbb{R}^{r}$, $r\geq n$, is a random vector whose distribution is
known or from which it can be simulated. The reason that we adopt the
formulation of (1) is that this model specification is very general. For
example, it includes the traditional statistical model specification using a
density function as its special case: Let $Y\sim f_{\theta_{0}}(y)$ be a
sample from a density function $f_{\theta}(y)$ with true parameter value
$\theta_{0}$. We then express $Y$ in the form of (1) as $Y=G(\theta_{0},U)$
where $G(\theta,U)=F_{\theta}^{-1}(U)$, $U\sim U(0,1)$ and $F_{\theta}(y)$ is
the cumulative distribution function of $f_{\theta}(y)$. Additional examples,
including more complicated ones, will be provided throughout the paper.

The repro samples method uses a simple yet fundamental idea: study the
performance of artificial samples that are obtained by mimicking the sampling
mechanism of the observed data; the artificial samples then help to quantify
the uncertainty in estimation of the associated model and parameters.
Specifically, let

$\bm{y}_{obs}=G({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits}_{0},\bm{u}^{rel})$

be the realized observations, where $\bm{u}^{rel}$ is the corresponding
unobserved (thus unknown) realization of $\bm{U}$. For any
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits\in\Theta$ and a copy
of $\bm{u}^{*}\sim\bm{U}$, we can obtain a copy of artificial sample
$\bm{y}^{*}=G({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits},\bm{u}^{*})$,
which we refer to in this paper as repro samples. Note that, when
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits=\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}$
and $\bm{u}^{*}$ matches $\bm{u}^{rel}$, the corresponding repro sample
$\bm{y}^{*}=G(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits,\bm{u}^{*})$
will match $\bm{y}_{obs}$. However, both
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}$ and
$\bm{u}^{rel}$ are unknown. Since $\bm{u}^{rel}$ and $\bm{u}^{*}$ are both
realizations from the same $\bm{U}$ distribution, we use this fact, together
with a so-called nuclear mapping function on ${\cal
U}\times\Theta\to\mathbb{R}^{q}$, $q\leq n$, to quantify how likely a
$\bm{u}^{*}$ is a potential match for $\bm{u}^{rel}$. For the collection of
$\bm{u}^{*}$’s that are deemed as “likely” matches, we look for the
corresponding $\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$ that
produces a matching repro sample
$\bm{y}^{*}=G(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits,\bm{u}^{*})$
with $\bm{y}^{*}=\bm{y}_{obs}$. We use these
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$ to help recover
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}$. We formally
formulate the procedure and show that a subset of $\Theta$ constructed by this
procedure has a desired frequency coverage rate and thus is a confidence set
for $\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}$.

In Section 2, we formally describe the repro samples framework. For a given
potential value of $\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$,
we construct a fixed level-$\alpha$ Borel set of $\bm{u}^{*}$ through a
nuclear mapping function, and analyze whether any $\bm{u}^{*}$ in the Borel
set produces repro samples
$\bm{y}^{*}=G({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits},{\bm{u}}^{*})$
that match $\bm{y}_{obs}$. By controlling $\bm{u}^{*}$ through the fixed Borel
set and seeking for the matching, we obtain a theoretically guaranteed
level-$\alpha$ confidence set for
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}$. The use of
nuclear mapping functions is new and general, and it greatly increases the
flexibility and utility of the repro samples method. If the distribution of
the nuclear mapping has an explicit form, so does the confidence set. When the
the confidence set cannot be explicitly expressed, we provide a generic
algorithm. Moreover, the repro samples method can also be used to obtain
p-values to conduct hypothesis testing for a wide range of problems. In
addition, we further extended the framework beyond the parametric model (1),
where the data generating model
$G(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits,\bm{U})$ is not
completely known. Although the repro samples method shares many insights with
some Bayesian and fiducial procedures, its development is purely under the
frequentist paradigm. It sidesteps the need to specify a likelihood or a loss
function, and also the need to use large sample theories or the Dumpster-
Shafer calculus. It has a wide range of applications, including problems
involving any type of parameter spaces, exact and asymptotic inferences,
parametric and nonparametric models, among others.

Section 3 provides additional developments we use to facilitate our
understanding and broaden the application of the repro samples method. In
Section 3.1, we explore the relationship of our repro samples method with the
classical Neyman-Pearson hypothesis testing procedure by considering a special
case in which we choose a test statistic as our nuclear mapping function. In
general, a nuclear mapping function does not need to be a test statistics.
Even in special cases in which the nuclear mapping is defined by a typical
test statistics, we show that a confidence set obtained by the repro samples
method is either the same or smaller than the one obtained by inverting the
classical Neyman-Pearson hypothesis testing procedure.

We also obtain a set of optimality results. Furthermore, we develop in Section
3.2 a general technique to handle nuisance parameters through profiling. The
method maintains the desired coverage rates for confidence sets constructed
for target parameters, even in finite sample settings. Finally in Section 3.3,
we improve the computational efficiency by investigating the use of data-
dependent candidate sets to greatly reduce the search region in the parameter
spaces. Using a candidate set is particularly effective when the target
parameter of the inference is discrete, since we can take advantage of a many-
to-one mapping that is uniquely inherent in a repro samples procedure. We
provide a Monte-Carlo approach to obtain the candidate set and supporting
theories to ensure the coverage rate of the final confidence set.

To demonstrate the utility and performance of the repro samples method, we use
a normal mixture model as a case study example in Section 4. Accounting for
uncertainty in an estimation of the unknown number of components in a mixture
model, say $\tau_{0}$, is a notoriously difficult and unsolved inference
problem (e.g., Roeder, 1994; Chen et al., 2012), partly because the parameter
space of $\tau_{0}$ is discrete and the large sample CLT does not apply on any
point estimator $\hat{\tau}$ of $\tau_{0}$. Here, we construct a finite-sample
confidence set for the number of components $\tau_{0}$ in the normal mixture
model. We demonstrate that the repro samples method can effectively solve this
open inference problem while handling a large and varying number of nuisance
parameters. We perform simulation and real data analyses, and provide
comparisons to existing frequentist and Bayesian approaches. The analyses is
revealing. First, the regular criterion-based point estimators (even
consistent estimators such as the BIC estimators) are often biased for smaller
$\tau$ and they rarely recover the true $\tau_{0}$. Second, existing
(frequentist) likelihood-based procedures do not provide any upper bounds, and
they have low power and work only asymptotically. Third, Bayesian credible
sets obtained using the posteriors of $\tau$ perform poorly in terms of
covering $\tau_{0}$ in repeated runs, even when a uniform flat prior of $\tau$
(over a range including the true value $\tau_{0}$) is used. In addition,
credible sets for $\tau_{0}$ are very sensitive to the prior choices of both
the target parameter $\tau$ and also (surprisingly) the nuisance parameters.
By contrast, the proposed repro sample confidence sets always recover
$\tau_{0}$ with the desired coverage rate. They also provide a full picture of
the sampling uncertainty pertaining to the estimation of $\tau_{0}$.

The repro samples method stems from several existing artificial-sample-based
inference procedures across Bayesian, frequentist and fiducial (BFF)
paradigms, particularly, approximate Bayesian computing (ABC), generalized
fiducial inference (GFI) and bootstrap methods, in which artificial samples
are used to help quantify inference uncertainty. It also extends a key idea of
the inferential model (IM) (Martin and Liu, 2015) that the inference
uncertainty is quantified through the auxiliary $\bm{U}$. We view the repro
samples approach as an improvement over the current existing inference methods
and provide comparisons throughout the paper when the topic arises, At the
outset, we emphasize that the repro samples method side-steps many
restrictions encountered in the aforementioned BFF approaches. First,
bootstrap, GFI and ABC methods rely on large sample CLT or some version of
Bernstein von Mises-type theorem (thus CLT) to justify their (frequentist)
frequency performance, therefore does not guarantee Small sample frequency
performance. Second, both ABC and GFI methods are approximation methods that
require a predetermined threshold $\epsilon$ to judge whether an artificial
sample is sufficiently close to the observed sample. Their supporting theories
need $\epsilon\to 0$ at a fast rate as the sample size $n\to\infty$. However,
when $\epsilon\to 0$ at a fast rate, the sample rejection rate in their
algorithms tends to $1$. Thus, how to choose an appropriate $\epsilon$ to
balance the procedures’ computational efficiency and inference validity is an
unsolved problem (Li and Fearnhead, 2018). Finally, IM method is designed for
epistemic probabilistic inference, and its complexity and reliance on
Dempster-Shafer calculus make it hard to implement under complex settings. The
IM method is not easily adaptable for asymptotic inference, and the
restriction of using an ancillary to quantify uncertainty also limits its use
in applications. The repro samples method, developed fully following the
frequentist setting, addresses or sidesteps all these issues.

The remainder of this article is arranged as follows. Section 2 develops the
general framework of the repro samples method with supporting theories and
associated algorithms. Sections 3 considers several important aspects that
further support the general development, including connection to the Neyman-
Pearson procedure, handling of nuisance parameters, and a unique approach of
getting a data-driven candidate set to improve computational efficiency.
Section 4 is a case study on a normal mixture model that addresses a long
standing open problem in statistical inference. Section 5 contains further
remarks and discussions.

##  2 A general inference framework by repro samples

Let $(\Omega,{\cal F},\mathbb{P})$ be a probability space and $\bm{U}\in{\cal
U}$ is a measurable $r\times 1$ random vector on $\Omega$ with
$\\{\bm{U}\leq\bm{u}\\}=\\{\omega\in\Omega:\bm{U}(\omega)\leq\bm{u}\\}$.
Suppose $\bm{Y}$ is the sample data from (1) and the $m\times 1$ vector
$\bm{Z}={\bf m}(\bm{Y})$, $m\leq n$, is a function (or a summary) of the
sample data $\bm{Y}$. A slightly more general model than (1) is to assume that
$\bm{Z}$ is directly generated from

$\bm{Z}=G_{z}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits}_{0},\bm{U}),$
(2)

where $G_{z}(\cdot,\cdot)={\bf m}(G(\cdot,\cdot))$ is a given mapping from
$\Theta\times{\cal U}\to{\cal Z}\subset R^{m}$. In the special case that ${\bf
m}(\cdot)$ is the identity mapping function, then $G_{z}=G$, $\bm{Z}=\bm{Y}$
and model (2) reduces to model (1); thus (1) is a special case of (2).

Let $\bm{z}_{obs}$ be the observed $\bm{Z}$ and $\bm{u}^{rel}$ be the
(unknown) realization of $\bm{U}$ associated with $\bm{z}_{obs}$. By (2), we
have

$\bm{z}_{obs}=G_{z}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits}_{0},\bm{u}^{rel}).$

For any given $\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits^{*}$ and
a simulated $\bm{u}^{*}\sim\bm{U}$, we can thereafter construct an artificial
(‘fake’) copy of $\bm{z}_{obs}$:

$\bm{z}^{*}=G_{z}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits}^{*},\bm{u}^{*}).$

We call $\bm{z}^{*}$ a repro sample of $\bm{z}_{obs}$ and
$({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits}^{*},\bm{u}^{*})$ a
repro copy of
$({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits}_{0},\bm{u}^{rel})$.

The goal of Section 2.1-2.3 is to present how to use repro samples to
construct a confidence set, based on the observed $\bm{z}_{obs}$, for the
unknown parameter $\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}$
with a guaranteed coverage rate. In Section 2.1 we consider the case of giving
the algorithmic model (2) where $G_{z}(\cdot,\cdot)$ is completely known.
Section 2.2 contains a generic algorithm when an explicit mathematical
expression of the confidence set is not available. Section 2.3 extends the
developments to a more general case in which the algorithmic model
$G_{z}(\cdot,\cdot)$ is not completely given.

###  2.1 A general formulation of a basic repro samples method

Under the algorithmic model (2) with a fixed
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}$, the sampling
uncertainty of $\bm{Z}$ is solely determined by the uncertainty of $\bm{U}$,
whose distribution is free of the unknown model parameters
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}$. Let us first
focus on the easier task of quantifying the uncertainty in $\bm{U}$.
Specifically, let $B_{\alpha}\subset{\cal U}$ be a level-$\alpha$ Borel set
such that $\mathbb{P}\\{\bm{U}\in B_{\alpha}\\}\geq\alpha.$ Although
$\bm{u}^{rel}$ is unknown, we have more than $\alpha$-level confidence that
$\bm{u}^{rel}\in B_{\alpha}$. If $\bm{u}^{*}$ is confined within the same set
$B_{\alpha}$ (i.e., requiring $\bm{u}^{*}\in B_{\alpha}$) and we collect any
potential $\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$ values that
can create a repro sample
$\bm{z}^{*}=G_{z}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits,\bm{u}^{*})$
that matches $\bm{z}^{*}=\bm{z}_{obs}$, then it leads to a subset on $\Theta$:

$\Gamma_{\alpha}(\bm{z}_{obs})=\big{\\{}\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits:\exists\,\bm{u}^{*}\in
B_{\alpha}\mbox{ s.t.
}\bm{z}_{obs}=G_{z}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits},\bm{u}^{*})\big{\\}}\subset\bm{\Theta},$
(3)

where s.t. is short for such that. In another words, for a potential value
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$, if there exists a
$\bm{u}^{*}\in B_{\alpha}$ such that the repro sample
${\bm{z}}^{*}=G({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits},{\bm{u}}^{*})$
matches $\bm{z}_{obs}$ with $\bm{z}^{*}=\bm{z}_{obs}$, then we keep this
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$ in the set. Since
$\bm{z}_{obs}=G_{z}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits}_{0},\bm{u}^{rel})$,
if $\bm{u}^{rel}\in B_{\alpha}$ then
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}\in\Gamma_{\alpha}(\bm{z}_{obs})$.
Similarly, under the model
$\bm{Z}=G_{z}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits}_{0},\bm{U})$,
we have $\left\\{\bm{U}\in
B_{\alpha}\right\\}\subseteq\left\\{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}\in\Gamma_{\alpha}(\bm{Z})\right\\}$.
Thus, by construction,

$\mathbb{P}\left\\{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}\in\Gamma_{\alpha}(\bm{Z})\right\\}\geq\mathbb{P}\big{\\{}\bm{U}\in
B_{\alpha}\big{\\}}\geq\alpha.$

The set $\Gamma_{\alpha}(\bm{z}_{obs})$ constructed in (3) is a $100\alpha\%$
confidence set for
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}$. Here,
$B_{\alpha}$ is any single fixed set. As long as $\mathbb{P}(\bm{U}\in
B_{\alpha})\geq\alpha$, the validity statement holds.

To expand the scope of this development, the requirement $\mathbb{P}(\bm{U}\in
B_{\alpha})\geq\alpha$ is replaced by

$\mathbb{P}\left\\{T(\bm{U},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\in
B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\right\\}\geq\alpha,$
(4)

for any $\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$, where
$T(\cdot,\cdot)$ is a mapping function from ${\cal U}\times\Theta\to$ ${\cal
T}\subseteq\mathbb{R}^{q}$, for some $q\leq n$. Then, the constructed
confidence set in (3) becomes

$\Gamma_{\alpha}(\bm{z}_{obs})=\big{\\{}\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits:\exists\,\bm{u}^{*}\in{\cal
U}\mbox{ s.t.
}\bm{z}_{obs}=G_{z}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits},\bm{u}^{*}),\,T(\bm{u}^{*},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\in
B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\big{\\}}\subset\Theta.$
(5)

The function $T(\cdot,\cdot)$, which we refer to as a nuclear mapping, adds
much flexibility to our repro-sampling method, as we will show throughout the
paper. For the moment, we assume that the nuclear mapping $T$ is given;
further discussions are provided in later sections.

From now on the set $\Gamma_{\alpha}(\bm{z}_{obs})$ refers to the general
expression of (5) instead of the special case in (3), unless specified
otherwise. In the construction of $\Gamma_{\alpha}(\bm{z}_{obs})$ in (5), for
any potential value $\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$,
we impose two constraints: (a) The repro error $\bm{u}^{*}$ is confined by the
restriction
$T(\bm{u}^{*},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\in
B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$; (b) The
repro sample
$\bm{z}^{*}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)=G_{z}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits},\bm{u}^{*})$
generated using $\bm{\theta}$ matches the observed sample, i.e.,
$\bm{z}_{obs}=\bm{z}^{*}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$.
The quantification of inference uncertainty is achieved only through
constraint (a) requiring
$T(\bm{u}^{*},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\in
B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ for the
given $\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$, and it does
not involve the unknown true
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}$ that generates
$\bm{z}_{obs}$. The input of the observed sample $\bm{z}_{obs}$ is discernible
by matching the artificial repro sample
$\bm{z}^{*}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)=G_{z}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits},\bm{u}^{*})$
with $\bm{z}_{obs}$ in constraint (b).

As an illustration of the method, we consider below a simple binomial example:

######  Example 1 (Binomial sample).

Assume $Y\sim Binomial(r,\theta_{0})$, $0<\theta_{0}<1$, and we observe
$y_{obs}$. We can re-express the model as $Y=\sum_{i=1}^{r}{\bf
1}{(U_{i}\leq\theta_{0})}$ for $U_{1},\ldots,U_{r}\sim U(0,1)$. It follows
that $y_{obs}=\sum_{i=1}^{r}{\bf 1}{(u_{i}^{rel}\leq\theta_{0})},$ where
$\bm{u}^{rel}=(u_{1}^{rel},\ldots,u_{r}^{rel})$ is the realization of
$\bm{U}=(U_{1},\ldots,U_{r})$ corresponding to $y_{obs}$.

For a given vector $\bm{u}=(u_{1},\ldots,u_{r})$, we consider the nuclear
mapping function $T(\bm{u},\theta)=\sum_{i=1}^{r}{\bf 1}{(u_{i}\leq\theta)}$,
and we have $T(\bm{U},\theta)=\sum_{i=1}^{r}{\bf 1}{(U_{i}\leq\theta)}\sim
Binomial(r,\theta)$, for each given $\theta\in\Theta=(0,1)$. Let
$B_{\alpha}(\theta)=\\{\bm{u}\big{|}a_{L}(\theta)\leq T(\bm{u},\theta)\leq
a_{U}(\theta)\\}$, where

$\left(a_{L}(\theta),a_{U}(\theta)\right)=\underset{\left\\{(i,j):\sum^{j}_{k=i}{r\choose
k}\theta^{k}(1-\theta)^{(r-k)}\geq\alpha\right\\}}{\arg\min}|j-i|,$ (6)

is the pair of $(i,j)$ in $\\{(i,j):\sum^{j}_{k=i}{r\choose
k}\theta^{k}(1-\theta)^{(r-k)}\geq\alpha\\}$ that makes the shortest interval
$\left[a_{L}(\theta),a_{U}(\theta)\right]$. By a direct calculation, it
follows that
$P\left\\{T(\bm{U},\theta)\in[a_{L}(\theta),a_{U}(\theta)]\right\\}\geq\alpha$,
for any $\theta\in(0,1)$. The confidence set in (5) is then

$\displaystyle\Gamma_{\alpha}(y_{obs})$ $\displaystyle=$
$\displaystyle\bigg{\\{}\theta\big{|}\exists\,\bm{u}^{*}\in{\cal
U}\,\hbox{s.t.}\,\,y_{obs}=\sum_{i=1}^{r}{\bf
1}{(u_{i}^{*}\leq\theta)},\sum_{i=1}^{r}{\bf
1}{(u_{i}^{*}\leq\theta)}\in[a_{L}(\theta),a_{U}(\theta)]\bigg{\\}}$
$\displaystyle=$ $\displaystyle\left\\{\theta\big{|}y_{obs}=\sum_{i=1}^{r}{\bf
1}{(u_{i}^{*}\leq\theta)},y_{obs}\in[a_{L}(\theta),a_{U}(\theta)],\exists\,\bm{u}^{*}\in{\cal
U}\right\\}$ $\displaystyle=$
$\displaystyle\left\\{\theta\big{|}y_{obs}=\sum_{i=1}^{r}{\bf
1}{(u_{i}^{*}\leq\theta)},\exists\,\bm{u}^{*}\in{\cal
U}\right\\}\cap\left\\{\theta\big{|}y_{obs}\in[a_{L}(\theta),a_{U}(\theta)]\right\\}$
$\displaystyle=$ $\displaystyle\left\\{\theta\big{|}a_{L}(\theta)\leq
y_{obs}\leq a_{U}(\theta)\right\\}.$

The last equation holds, because for any given $\theta\in\Theta=(0,1)$ there
always exists at least a $\bm{u}^{*}\in{\cal U}$ such that
$y_{obs}=\sum_{i=1}^{r}{\bf 1}{(u_{i}^{*}\leq\theta)}$, i.e.,
$\\{\theta\big{|}y_{obs}=\sum_{i=1}^{r}{\bf
1}{(u_{i}^{*}\leq\theta)},\exists\,\bm{u}^{*}\in{\cal U}\\}=\Theta$.

Table 1 provides a numerical study comparing the empirical performance of
$\Gamma_{.95}(y_{obs})$ (Repro) against 95% confidence intervals obtained
using the traditional Wald method and the fiducial approach (GFI) of Hannig
(2009), both of which are asymptotic methods that can ensure coverage only
when $n$ is large. The rerpo sampling method is an exact method that improves
the performance of existing methods, especially when $n\theta_{0}<5$ (a rule
of thumb on asymptotic approximation of binomial data; cf., Casella and Berger
(1990, p.106)).

|  | $n=20,\theta_{0}=0.1$ |  | $n=20,\theta_{0}=0.4$ |  | $n=20,\theta_{0}=0.8$
---|---|---|---|---|---|---
|  | Coverage | Width |  | Coverage | Width |  | Coverage | Width
Repro |  | 0.949(0.007) | 0.281(0.059) |  | 0.963(0.006) | 0.408(0.026) |  | 0.959(0.006) | 0.342(0.045)
Wald |  | 0.877(0.010) | 0.236(0.052) |  | 0.927(0.008) | 0.418(0.028) |  | 0.915(0.009) | 0.332(0.071)
GFI |  | 0.988(0.003) | 0.293(0.065) |  | 0.963(0.006) | 0.438(0.022) |  | 0.973(0.004) | 0.365(0.056)

Table 1: Comparison of $95\%$ confidence intervals by the repro, Wald and GFI
methods in Binomial$(n,\theta_{0})$ data; repetitions $=1000$; standard errors
are enclosed in brackets.

The next theorem states that $\Gamma_{\alpha}(\bm{z}_{obs})$ in (5) is a
level-$\alpha$ confidence set.

######  Theorem 1.

Assume model (2) holds. If inequality (4) holds exactly, then the following
inequality holds exactly

$\mathbb{P}\left\\{{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits}_{0}\in\Gamma_{\alpha}(\bm{Z})\right\\}\geq\alpha\,\,\,\hbox{for}\,\,0<\alpha<1.$
(7)

Furthermore, if the inequality (4) holds approximately with
$\mathbb{P}\big{\\{}T(\bm{U},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\in
B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\big{\\}}\geq\alpha\\{1+o(\delta^{{}^{\prime}})\\}$,
then (7) holds approximately with
$\mathbb{P}\left\\{{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits}_{0}\in\Gamma_{\alpha}(\bm{Z})\right\\}\geq\alpha\\{1+o(\delta^{{}^{\prime}})\\}$,
for $0<\alpha<1$, where $\delta^{{}^{\prime}}>0$ is a small value that may
depend on sample size $n$.

We provide a proof of the theorem is in the Appendix. The
$\delta^{{}^{\prime}}$ in Theorem 1 may or may not depend on sample size $n$.
In examples involving large sample approximations, $\delta^{{}^{\prime}}$ is
often a function of $n$ with $\delta^{{}^{\prime}}\to 0$ as $n\to 0$. However,
there are also examples in which $\delta^{{}^{\prime}}$ does not involve $n$.
For example, suppose that $Y|\theta=\lambda\sim\text{Poisson}(\lambda)$. Then,
$U=(Y-\lambda)/\sqrt{\lambda}\to N(0,1)$, when $\lambda$ is large. So, if we
take $T(U,\lambda)=U$, $\mathbb{P}\\{T(U,\lambda)\in B\\}=\int_{{t}\in
B}\phi({t})d{t}\\{1+o(\lambda^{-1})\\}$, for any Borel set $B$, thus
$\delta^{{}^{\prime}}=o(\lambda^{-1})$. Here, $\phi(t)$ is the density
function of $N(0,1)$.

Moreover, the repro samples method can be used to do hypothesis testing. From
Theorem 1, we have the following corollary. A proof is given in Appendix.

######  Corollary 1.

Suppose we are interested in testing a hypothesis $H_{0}:\theta\in\Theta_{0}$
vs $H_{1}:\theta\not\in\Theta_{0}$. In that case, we can define a p-value
$p(\bm{z}_{obs})$ as

$\displaystyle
p(\bm{z}_{obs})=1-\inf_{\theta\in\Theta_{0}}\left[\inf\left\\{\alpha^{\prime}:\theta\in\Gamma_{\alpha^{\prime}}(\bm{z}_{obs})\right\\}\right],$

where $\Gamma_{\alpha^{\prime}}(\bm{z}_{obs})$ is the level-$\alpha^{\prime}$
repro samples confidence set defined in (5). Rejecting $H_{0}$ when
$p(\bm{z}_{obs})\leq\gamma$ leads to a size $\gamma$ test, for any
$0<\gamma<1$.

Although the repro samples method is developed following the frequentist
principles, thus making it a completely frequentist approach, it is also
closely related to several existing approaches across the Bayesian, fiducial,
and frequentist (BFF) paradigms. In Remarks 2.1 - 2.3 below, we provide
discussions to connect and compare these approaches. See also Appendix II for
further and more elaborated discussions.

Remark 2.1 (comparison with Neyman-Pearson testing method). Corollary 1
suggests that we can conduct hypothesis testing using the repro samples
method. Section 3.1 further investigates the relationship between the repro
samples method and the classical Neyman-Pearson (N-P) hypothesis testing
procedure under a special case where a nuclear mapping $T(\bm{u},\theta)$ in
(5) is chosen through a test statistic, i.e., we set
$T(\bm{u},\theta)=\widetilde{T}(\bm{z}^{\prime},\theta)$ where
$\bm{z}^{\prime}=G_{z}(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$
is a copy of data generated by the given
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$ and
$\widetilde{T}(\bm{z}^{\prime},\theta)$ is a test statistic of a hypothesis
testing problem. In general, the repro samples approach is broader than the
classical N-P procedure. In particular, the nuclear mapping
$T(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ does not
need to be a test statistic since it is not necessarily a function of data
$\bm{z}^{\prime}.$ Allowing the nuclear mapping to be a function of $\bm{u}$
provides more flexibility for developing inference procedures than the
conventional testing approach. We illustrate this point using a simple example
where a single $\text{Bernoulli}(\theta)$ observation $y_{obs}$ is generated
by $Y=I(U<\theta_{0})$, $U\sim\text{Uniform}(0,1).$ To get a confidence set
for $\theta_{0}$, we can let $T(u,\theta)=u$, and define
$B_{\alpha}(\theta)=\begin{cases}[1-\alpha,1]&\theta>0.5\\\
[\frac{1-\alpha}{2},\frac{1+\alpha}{2}]&\theta=0.5\\\
[0,\alpha]&\theta<0.5\end{cases}.$ Then, following (5), we construct the repro
samples confidence set
$\Gamma_{\alpha}(y_{obs})=\begin{cases}[1-\alpha,1]&y_{obs}=1\\\
[0,\alpha]&y_{obs}=0\end{cases}.$ Apparently, $U$ is not a function of $Y$,
since it is impossible to solve for $U$ in the structure equation
$Y=I(U<\theta)$ when given $(Y,\theta)$. In this case, $T(U,\theta)=U$ is not
a test statistic. Section 3.1 provides further discussion and additional
examples comparing the repro samples methods and the classical N-P method. We
show that the repro samples method can always produce an either equal or
smaller confidence set than the traditional N-P method when the nuclear
mapping $T(\bm{u},\theta)$ is chosen through a test statistic. There are also
options in which the nuclear mapping function $T(\bm{u},\theta)$ is not a test
statistic.

Remark 2.2 (comparison with GFI and ABC methods) In modern statistical
practice, R.A. Fisher’s fiducial method is understood as an inversion method
that solves the structural equation
$\bm{z}_{obs}=G_{z}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits},\bm{u}^{*})$
for parameter ${\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits}$ with
any $\bm{u}^{*}\sim\bm{U}$ (cf., Brenner et al., 1983; Hannig et al., 2016;
Thornton and Xie, 2022). The matching equation
$\bm{z}_{obs}=G_{z}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits},\bm{u}^{*})$
in (5) plays a key role in the GFI development (Hannig et al., 2016). Hannig
et al. (2016) and Thornton and Xie (2022) have also explored the connection of
GFI to the Bayesian ABC method, noting that trying to match $\bm{z}_{obs}$
with
$\bm{z}^{*}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)=G_{z}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits},\bm{u}^{*})$
is a key step in both GFI and ABC procedures. However, as stated in Hannig et
al. (2016), an exact matching of
$\bm{z}_{obs}=\bm{z}^{*}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$
for any given $\bm{u}^{*}$ is difficult and sometime even impossible. Both GFI
and ABC methods adopts a tuning threshold value $\epsilon$ to judge whether an
artificial sample
$\bm{z}^{*}=G_{z}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits},\bm{u}^{*})$
is close to the observed $\bm{z}_{obs}$. Unfortunately, the use of
approximating threshold $\epsilon$ has piratical issues. For instance, as
shown in Li and Fearnhead (2018), the validity of the ABC method requires that
$\epsilon\to 0$ at a fast rate, but this fast rate in turn leads to
inefficient sampling with a high (often 100%) rejection rate. How to choose an
appropriate $\epsilon$ to balance inferential validity and computational
efficiency is still an unaddressed question in both ABC and GFI. Our repro
sample approach avoids this issue of using a threshold $\epsilon$ by directly
working with a set of
$\bm{u}^{*}\in\\{T(\bm{u}^{*},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\in
B_{\alpha}(\theta)\\}$ and multiple copies of
$\bm{z}^{*}=G_{z}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits},\bm{u}^{*})$
for each given $\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$. As
further explained in Appendix II the repro samples method in effect compares
$\bm{z}_{obs}$ with multiple copies of
$\bm{z}^{*}=G_{z}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits},\bm{u}^{*})$
at each given $\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$ and
thus enabling us to use a confidence (significant) level $\alpha$ to replace
the arbitrary $\epsilon$. In addition, the repro sample method can provide
finite sample inference whereas GFI does not. Finally, ABC is a Bayesian
procedure that not only requires additional prior assumption, but also
requires a sufficient summary statistic in its algorithm. The repro samples
approach compares favorably in that it does not have these constraints.

Remark 2.3 (comparison with IM method). Inferential model (IM) (Martin and
Liu, 2015) also works directly on $\bm{U}$ to quantify inference uncertainty,
but contains notable differences from the proposed repro samples method. First
and foremost, the IM method exists to provide an epistemic probabilistic
inference for $\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}$, a
task not considered under the current frequentest framework; cf., Martin and
Liu (2015). The repro samples on the other hand is a frquentist approach
developed fully within the frequentist framework. Also importantly, the IM
method uses random sets to quantify the uncertainty of the auxiliary term
$\bm{U}$. But the repro samples method quantifies inference uncertainty
through a nuclear mapping function
$T(\bm{U},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ that
involves a potential parameter
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$ and an associated
Borel set. Only one fixed level-$\alpha$ Borel set is needed. Although IM can
sometimes potentially provide some finer details in special cases (since
Dempster-Shefer calculus automatically includes both upper and lower
probabilities; i.e., plausible and belief functions), with additional
constraints, the application of IM is more limited than repro samples method,
especially under more complicated settings. Moreover, using a nuclear mapping
function involving $\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$
greatly extends the scope and flexibility of the repro samples method. As a
special case, we are able to directly work on auxiliary $\bm{U}$ by setting
$T(\bm{U},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)=\bm{U}$. To
the best of our knowledge, the framework of IM method cannot be easily
extended to incorporate varying
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$. Furthermore, IM
development requires Dempster-Shefer calculus and asymptotic inference is not
available under the development. When both apply, the IM and repro samples
methods provide either the same or comparable inference conclusions. Finally,
the repro samples development affords us the opportunity to explore the unique
aspect of matching to obtain a data-driven candidate set and greatly reduce
computational burden; details are provided in Section 3.3. On the other hand,
how to incorporate a data-driven candidate set with the IM method to reduce
computational cost is not apparent.

###  2.2 A general Monte-Carlo algorithm

In Example 1, the distribution of
$T({\bm{U}},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ has an
explicit form, enabling us to express the level-$\alpha$ confidence set
$\Gamma_{\alpha}(\bm{y}_{obs})$ in (5) explicitly. In cases when the
distribution of
$T({\bm{U}},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ is not
explicit, we can use a Monte-Carlo method to obtain a level-$\alpha$ Borel set
$B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ in (4)
and then construct the level-$\alpha$ confidence set
$\Gamma_{\alpha}(\bm{y}_{obs})$ in (5) as stated in Algorithm 1 below:

Algorithm 1 Monte-Carlo algorithm (general case)

1:Simulate many copies of ${\bm{u}^{s}}\sim\bm{U}$; the collection of
${\bm{u}^{s}}$ forms a set ${\cal V}$.

2:For each given (grid) value of
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits\in\Theta$,
(a) Compute
$\\{T({\bm{u}^{s}},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits),{\bm{u}^{s}}\in{\cal
V}\\}$, from which obtain an empirical distribution of
$T(\bm{U},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$. Based on
the empirical distribution, we get a level-$\alpha$ Borel set
$B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$.
(b) Check whether there exists a $\bm{u}^{*}\in{\cal U}$ such that
${\bm{z}_{obs}}=G_{z}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits},\bm{u}^{*})$
and $T(\bm{u}^{*},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\in
B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$. If both
of the above criteria are satisfied, keep the
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$.

3:Collect all kept $\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$ to
form the set $\Gamma_{\alpha}(\bm{y}_{obs})$.

This general algorithm, grid search algorithm, searches the space of $\Theta$
and is used only when an explicit expression of
$T(\bm{U},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$’s
distribution is not available. Our remarks on the implementation of the
algorithm are as follows.

Remark 2.4 To obtain the level-$\alpha$ Borel set
$B_{\alpha}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits})$ for a
given $\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$ in Step 2(a), w
often need to compute the (empirical) distribution of
$T(\bm{U},{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits})$ based on
the set of Monte-Carlo points ${{\cal
S}_{\theta}}=\\{T(\bm{u}^{s},{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits})$,
$\bm{u}^{s}\in{\cal V}\\}$. Then according to the empirical distribution, we
can always numerically calculate
$B_{\alpha}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits})$. If
$T(\bm{u},{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits})$ is a
scalar, we just directly compute its empirical distribution $F_{\cal V}(t)=$
$\sum_{\bm{u}^{s}\in{\cal
V}}I\big{(}T(\bm{u}^{s},{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits})\leq
t\big{)}\big{/}|{\cal V}|$ and use the (upper/lower) quantiles of $F_{\cal
V}(t)$ to obtain a level-$\alpha$ interval as the Borel set
$B_{\alpha}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits})$. When
$T(\bm{u},{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits})$ is a
vector, say in ${\cal T}\subseteq\mathbb{R}^{q}$, we can use data depth and
follow Liu et al. (2021, §3.2.2) to construct the Borel set
$B_{\alpha}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits})$.
Specifically, we let $D_{{\cal S}_{\theta}}(\bm{t})$ be the empirical depth
function for any $\bm{t}\in\cal T$ that is computed based on the Monte-Carlo
points in ${\cal
S}_{\theta}=\\{T(\bm{u},{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits})$,
$\bm{u}\in{\cal V}\\}$ and $F_{{\cal V}|D}(t)=$ $\sum_{\bm{u}^{s}\in{\cal
V}}I\big{(}D_{{\cal
S}_{\theta}}(T(\bm{u}^{s},{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits}))\leq
t\big{)}\big{/}|{\cal V}|$ be the empirical CDF of $D_{{\cal
S}_{\theta}}(T(\bm{u}^{s},{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits})).$
Then, following Lemma 1 (a) of Liu et al. (2021), the level-$\alpha$
(empirical) central region
$B_{\alpha}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits})=\\{\bm{t}:F_{{\cal
V}|D}(D_{{\cal S}_{\theta}}(\bm{t}))\geq 1-\alpha\\}$ is a level-$\alpha$
Borel set on ${\cal T}$. See also Example 2 (B) of Section 2.3 for a numerical
example, where the empirical depth function $D_{{\cal S}_{\theta}}(\cdot)$ is
computed by R package ddalpha (Pokotylo et al., 2019).

Remark 2.5 The algorithm needs to calculate the empirical distribution of
$T(\bm{U},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ for each
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$. To improve the
computing efficiency (at the cost of potentially losing the flexibility of the
repro samples method), we can sometimes choose a nuclear mapping function
$T(\bm{U},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ whose
distribution of is free (or approximately free) of
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$. In this case,
$B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ in (4) is
free (or approximately free) of
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$ with
$B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\equiv
B_{\alpha}$. A benefit of using this type of nuclear mapping function is that
we can modify Algorithm 1 to save computing time: in step 2(a), we only need
to obtain $B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$
for a single value of $\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$
instead of for each grid value of
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits\in\Theta.$

Remark 2.6 Later in Section 3.3, we will discuss the technique of using a
data-driven candidate set to significantly reduce the search space. From the
outset, we would like to comment that the search algorithm works in most cases
that one commonly encounter. In many statistical analysis, for instance the
commonly used likelihood-based inference, the regular conditions on the
parameter space include that $\Theta$ is a continuous space and the test
statistic is also continuous (and often differentiable) in its inner space
$\Theta^{o}$. Similarly, if the nuclear mapping function
$T(y_{obs},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ is
continuous in $\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$ and
$\Theta$ is a continuous space, one can in principle use a grid search method
to carry out an analysis in this situation. What is new in the repro samples
development is that $\Theta$ can either be a discrete space or a set of non-
numerical subjects. In this case, even if $\Theta$ is (countable) infinite, we
provide in Section 3.3 a way to take advantage of a many-to-one mapping
structure inherited in the repro samples procedure to create a finite
candidate set and carry out an efficient grid search method. See Section 3.3
for further details.

###  2.3 A further extension to more general cases

In this section, we further extend the framework developed above to settings
where the matching equation (2) of the algorithmic model is not available.
This includes some nonparametric inference problems that are not covered by
the existing ABC or GFI methods. Consider the example of making inference for
a population quantile of a completely unknown distribution $F$, say,
$\theta_{0}=F^{-1}(\zeta)$, the $\zeta$-th quantile parameter of $F$. Since
$I(Y<\theta_{0})\sim\text{Bernoulli}(\zeta)$ for $Y\sim F$, we have equation
$I(Y<\theta_{0})=U$, where $U\sim\text{Bernoulli}(\zeta)$. Now suppose we
observe data $\bm{y}_{obs}=(y_{1}\ldots,y_{n})^{\top}$, for a fixed $n$, then
the equation becomes

$\sum_{i=1}^{n}I(Y_{i}-\theta_{0}<0)=\sum_{i=1}^{n}U_{i}$ (8)

where $Y_{i}\overset{iid}{\sim}F$ and
$U_{i}\overset{iid}{\sim}\text{Bernoulli}(\zeta)$. The corresponding realized
version is $\sum_{i=1}^{n}I(y_{i}-\theta_{0}<0)=\sum_{i=1}^{n}u_{i}^{rel}$.
Equation (8) can not be expressed in the form of either (1) or (2). In fact,
one cannot generate a repro sample $\bm{y}^{*}$ based on (8) that is directly
comparable to the observed sample $\bm{y}_{obs}$. However, we can still use a
repo sampling approach to provide an exact inference for such a nonparametrc
inference problem.

In particular, we consider the following generalized data generating equation:

$g(\bm{Z},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0},\bm{U})=\bm{0},$
(9)

where $g$ is a given mapping function from ${\cal Z}\times\Theta\times{\cal
U}\to\mathbb{R}^{s}$ and $\bm{Z}={\bf m}(\bm{Y})$ is a $m\times 1$ vector as
defined at the start of Section 2. The corresponding realization version is
$g(\bm{z}_{obs},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0},\bm{u}^{rel})=\bm{0}$.
Equation (8) is a special case of (9) with
$g(\bm{Y},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0},\bm{U})=\sum_{i=1}^{n}I(Y_{i}-\theta_{0}<0)-\sum_{i=1}^{n}U_{i}$.

The level-$\alpha$ confidence set in (5) is modified to be

$\displaystyle\Gamma_{\alpha}(\bm{z}_{obs})=\big{\\{}\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits:\exists\,\bm{u}^{*}\in{\cal
U}\mbox{ s.t.
}g(\bm{z}_{obs},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits,\bm{u}^{*})=0,T(\bm{u}^{*},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\in
B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\big{\\}}\subset\Theta.$
(10)

We have the following theorem for the set $\Gamma_{\alpha}(\bm{z}_{obs})$
defined in (10). The proof is similar to that of Theorem 1 and is briefly
described in Appendix I.

######  Theorem 2.

Assume the generalized model equation (9) holds. If the inequality (4) holds
exactly, then for $\Gamma_{\alpha}(\bm{z}_{obs})$ defined in (10) the
following inequality holds exactly,

$\mathbb{P}\left\\{{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits}_{0}\in\Gamma_{\alpha}(\bm{Z})\right\\}\geq\alpha\quad\hbox{for
$0<\alpha<1$.}$ (11)

Furthermore, if (4) holds approximately with
$\mathbb{P}\left\\{T(\bm{U},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\in
B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\right\\}\geq\alpha\\{1+o(\delta^{{}^{\prime}})\\}$,
then (11) holds approximately with
$\mathbb{P}\left\\{{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits}_{0}\in\Gamma_{\alpha}(\bm{Z})\right\\}\geq\alpha\\{1+o(\delta^{{}^{\prime}})\\}$,
for $0<\alpha<1$, where $\delta^{{}^{\prime}}>0$ is a small value that may
depend on sample size $n$.

As discussed in Sections 2.2, when the distribution of
$T(\bm{u}^{*},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ cannot
be explicitly expressed, we may use a Monte-Carlo algorithm to carry out the
construction of $\Gamma_{\alpha}(\bm{y}_{obs})$. Specifically, we may still
use Algorithm 1 to get $\Gamma_{\alpha}(\bm{y}_{obs})$ but replace the
matching equation
${\bm{z}_{obs}}=G_{z}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits},\bm{u}^{*})$
in the algorithm with
$g(\bm{z}_{obs},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits,\bm{u}^{*})=0$.

The example below provides an illustration of this extended framework. Part
(A) is on nonparametric quantile inference and part (B) is on censored
quantile regression where making a finite-sample (semi-parametric) inference
for model parameters is new and previously unavailable in the literature. Both
parts demonstrate that the repro samples method performs well and is compared
favorably to the corresponding large sample bootstrap methods.

######  Example 2 (Nonparametric and semi-parametric inference of quantiles).

(A) Assume $\bm{y}_{obs}=(y_{1},\ldots,y_{n})^{\top}$ are from an unknown
distribution $F$. We are interested in making inference about the population
quantile $\theta_{0}=F^{-1}(\zeta)$, for a given $0<\zeta<1$. Write
$g(\bm{Y},\theta,\bm{U})=\sum_{i=1}^{n}I(Y_{i}\leq\theta)-\sum_{i=1}^{n}U_{i}$,
where $Y_{i}\sim F$ and $U_{i}\sim\text{Bernoulli}(\zeta)$. We define
$T(\bm{U},\theta)=T(\bm{U})=\sum_{i=1}^{n}U_{i}$. It follows that
$T(\bm{U})\sim\text{Binomial}(n,\zeta)$ and
$B_{\alpha}=[a_{L}(\zeta),a_{U}(\zeta)]$, where $a_{L}(\cdot)$ and
$a_{U}(\cdot)$ are defined in (6). By (10), a finite sample level-$\alpha$
confidence set is
$\Gamma_{\alpha}(\bm{y}_{obs})=\big{\\{}\theta:\sum_{i=1}^{n}I(y_{i}\leq\theta)-\sum_{i=1}^{n}u_{i}=0,\sum_{i=1}^{n}u_{i}\in[a_{L}(\zeta),a_{U}(\zeta)]\big{\\}}=\big{\\{}\theta:a_{L}(\zeta)\leq\sum_{i=1}^{n}I(y_{i}\leq\theta)\leq
a_{U}(\zeta))\big{\\}}=\big{[}y_{(a_{L}(\zeta))},y_{(a_{U}(\zeta)+1)}\big{)},$
where $y_{(k)}$ is the $k$th order statistic of the sample
$\\{y_{1},\ldots,y_{n}\\}$, $y_{(0)}$ and $y_{(n+1)}$ are the infimum and
supremum of support of $F.$

Figure 1 provides a numerical study to compare the empirical performance of
$\Gamma_{.95}(y_{obs})$ against 95% intervals obtained by the conventional
bootstrap method under two types of $F$. The repro samples method works well
even for $\zeta$’s close to $0$ or $1$, whereas the bootstrap method has
coverage issues under these settings.

(B) Consider a censored regression model (Powell, 1986), $Y_{i}=\max\\{0,$
${\bm{x}}^{\top}_{i}\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}+\epsilon_{i}\\}$,
$i=1,\ldots,n,$ where $\epsilon_{i}$ are independent error terms from an
unknown distribution $F$ with medium $0$. Following Bilias et al. (2000), the
estimating equation for performing a quantile regression is
$\sum_{i=1}^{n}\bm{x}_{i}\big{\\{}I(Y_{i}-{\bm{x}}_{i}^{T}\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits\leq
0)-\zeta\big{\\}}I(\bm{x}^{T}_{i}\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits>0)=0$,
where $\zeta$ is the quantile level and the parameter of interest is
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}(\zeta)=(F^{-1}(\zeta),0,\ldots,0)^{\top}+\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}$.
Since
$I(Y_{i}-{\bm{x}}_{i}^{T}\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}(\zeta)\leq
0)\sim\text{Bernoulli}(\zeta)$, equation (9) can be expressed as
$g(\bm{Y},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits,\bm{U})=\sum_{i=1}^{n}{\bm{x}}_{i}I{({\bm{x}}_{i}^{T}\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits>0)}\big{\\{}I(Y_{i}-{\bm{x}}_{i}^{T}\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits\leq
0)-I(U_{i}\leq\zeta)\big{\\}},$ where $U_{i}\sim U(0,1).$ We consider a
nuclear mapping function
$T(\bm{U},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)=\big{\\{}\sum_{i=1}^{n}\bm{x}_{i}\bm{x}_{i}^{T}I({\bm{x}}_{i}^{T}\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits>0)\big{\\}}^{-1}\sum_{i=1}^{n}{\bm{x}}_{i}I{({\bm{x}}_{i}^{T}\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits>0)I(U_{i}\leq\zeta)}$.
Since the distribution of
$T(\bm{U},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ does not
have an explicit expression, we use Algorithm 1 to obtain a level-$\alpha$
confidence set for
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}(\zeta)$. Here,
$T(\bm{U},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ is a vector
(of the same length of
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$) and
$B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ is
obtained using the data-depth approach described in Remark 2.4. Table 2
provides a numerical study to compare $\Gamma_{\alpha}(\bm{y}_{obs})$ in (10)
with the confidence sets obtained using the enhanced bootstrap method of
Bilias et al. (2000). Again, the repro samples method provides level-$95\%$
confidence sets with finite sample coverage guarantee and it works well for
this difficult inference problem even at a small $n=75$, while the asymptotic
bootstrap method has coverage issues under the same settings for $\zeta=0.5$.
For $\zeta=0.9,$ both methods provide good coverage rates, but the bootstrap
confidence sets are much larger.

Quantile | Error Distribution | Repro Samples | Bootstrap | Relative Volume
---|---|---|---|---
$\zeta=0.50$ | (a) Normal | 0.945(0.007) | 0.885(0.010) | 0.779(0.028)
(b) Mixture | 0.927(0.008) | 0.891(0.010) | 0.784(0.040)
(c) Heteroscedastic | 0.943(0.007) | 0.893(0.010) | 0.762(0.008)
$\zeta=0.90$ | (a) Normal | 0.945(0.007) | 0.970(0.005) | 17.518(2.723)
(b) Mixture | 0.949(0.007) | 0.977(0.005) | 214.962(9.595)
(c) Heteroscedastic | 0.951(0.007) | 0.980(0.004) | 346.202(11.403)

Table 2: Coverage rates and relative volumes of repro samples method versus
the bootstrap method by Bilias et al. (2000). The study is based on model
$y_{i}=\max\\{0.5+x_{i1}+x_{i2}+\epsilon_{i},0\\},$ where
$x_{i1}\sim\text{Bernoulli}(0.5)$ and $x_{i2}\sim N(0,1),$ and $n=75.$ The
$\epsilon_{i}$ are generated from: (a) normal distribution $\epsilon_{i}\sim
N(0,1),$ (b) mixture normal $\epsilon_{i}\sim
w_{i}N(0,1)+(1-w_{i})N(0,2),w_{i}\sim\text{Bernoulli}(0.75)$ and (c)
heterocedastic distribution $\epsilon_{i}\sim(1+0.15x_{i2})N(0,1).$ The
relative volume of the two three dimensional confidence sets $C_{repro}$ and
$C_{boot}$ is defined as $\frac{\text{the proportion of points in $C_{repro}$
that are also in $C_{boot}$}}{\text{the proportion of points in $C_{boot}$
that are also in $C_{repro}$}}$. Each setting is repeated 1000 times.

(a) Cauchy Distribution

(b) Negative Binomial Distribution

Figure 1: Coverage rates and widths of 95% confidence intervals for
$\theta_{0}=F^{-1}(\zeta)$ of an (unknown) distribution $F$ for a range of
percentile $\zeta$ values. Bars and dots colored in blue are produced by the
repro samples method and those in black are by the conventional bootstrap
method. In (a) $F$ is $\text{Cauchy}(0,1)$ and in (b) $F$ is
$\text{NB}(2,0.1)$. Both cases have sample size $n=60.$ The simulations are
replicated for $1000$ times. The error bars display the mean $\pm 2$ standard
errors.

##  3 Test statistic, nuisance parameters and candidate set

For any nuclear mapping function $T$, as long as we have a set
$B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ such that
(4) holds, the set $\Gamma_{\alpha}(\bm{z}_{obs})$ constructed accordingly is
a level-$\alpha$ confidence set for
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$. In other words, it
is possible that one might have multiple choices of $T$ that all lead to valid
(and perhaps different) confidence sets. Therefore, the role of $T$ under the
repro samples framework is similar to that of a test statistic under the
classical (Neyman-Pearson) hypothesis testing framework. Accordingly, a good
choice for $T$ is problem specific and we do not have a unique solution to all
situations. In this section, we investigate the choice of the nuclear mapping
function $T$ in two specific yet common settings. Specifically, in Section
3.1, we consider the special case of using a test statistic, including an
(approximate) pivot quantity such as a likelihood ratio test statistic, as a
nuclear mapping function. We provide comparisons to Neyman-Pearson testing
methods. In Section 3.2, we consider the situation of existing nuisance
parameters and provide a method to define a nuclear mapping function while
minimizing the impact of the nuisance parameters.

To improve the efficiency of computing Algorithms 1, we also discuss the use
of a pre-screening candidate set in Section 3.3. We pay special attention to
discrete target parameters, a unique and also important development under the
repro samples framework.

###  3.1 Define a nuclear mapping through a test statistic

In a classical testing problem
$H_{0}:\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}=\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$
versus
$H_{1}:\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}\not=\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$,
one often constructs a test statistic, say
$\widetilde{T}(\bm{z}_{obs},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$,
and assume that we know the distribution of
$\widetilde{T}(\bm{Z}^{\prime},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$
for
$\bm{Z}^{\prime}=G_{z}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits,\bm{U})$
generated under
$H_{0}:\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}=\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$.
Under the Neyman-Pearson framework, we can derive a level-$\alpha$ acceptance
region
$A_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)=\big{\\{}\bm{z}_{obs}\big{|}\widetilde{T}(\bm{z}_{obs},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\in
B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\big{\\}}$,
where the set
$B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ satisfies
$\mathbb{P}\big{\\{}\widetilde{T}(\bm{Z}^{\prime},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\in
B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\big{\\}}\geq\alpha$.
By the property of test duality, a level-$\alpha$ confidence set is

$\widetilde{\Gamma}_{\alpha}(\bm{z}_{obs})=\left\\{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits:\bm{z}_{obs}\in
A_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\right\\}=\left\\{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits:\widetilde{T}(\bm{z}_{obs},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\in
B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\right\\}.$

Suppose the nuclear mapping function
$T(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ is defined
through this test statistic:

$T(\bm{U},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)=\widetilde{T}(\bm{Z}^{\prime},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)=\widetilde{T}(G_{z}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits},\bm{U}),\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits),$
(12)

where
$\bm{Z}^{\prime}=G_{z}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits,\bm{U})$.
Since
$\mathbb{P}\left\\{T(\bm{U},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\in
B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\right\\}=\mathbb{P}\big{\\{}\widetilde{T}(\bm{Z}^{\prime},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\in
B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\big{\\}}\geq\alpha$
for the same $B_{\alpha}(\cdot)$, a level-$\alpha$ confidence set of (5) by
the repro samples method is

$\displaystyle\Gamma_{\alpha}(\bm{z}_{obs})=\left\\{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits:\exists\bm{u}^{*}\in{\cal
U}\mbox{ s.t.
}\bm{z}_{obs}=G_{z}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits},\bm{u}^{*}),\widetilde{T}(G_{z}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits},\bm{u}^{*}),\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\in
B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\right\\}.$

We have the following lemma suggesting
$\Gamma_{\alpha}(\bm{z}_{obs})\subseteq\widetilde{\Gamma}_{\alpha}(\bm{z}_{obs})$.
See Appendix I for a proof.

######  Lemma 1.

If the nuclear mapping function is defined through a test statistic as in
(12), then we have
$\Gamma_{\alpha}(\bm{z}_{obs})\subseteq\widetilde{\Gamma}_{\alpha}(\bm{z}_{obs})$.
The two sets are equal
$\Gamma_{\alpha}(\bm{z}_{obs})=\widetilde{\Gamma}_{\alpha}(\bm{z}_{obs})$,
when
$\widetilde{\Gamma}_{\alpha}(\bm{z}_{obs})\subseteq\big{\\{}\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits:\bm{z}_{obs}=G_{z}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits},\bm{u}^{*}),\exists\,\bm{u}^{*}\in{\cal
U}\big{\\}}$.

It is possible that $\Gamma_{\alpha}(\bm{z}_{obs})$ is strictly smaller than
$\widetilde{\Gamma}_{\alpha}(\bm{z}_{obs})$ with
$\Gamma_{\alpha}(\bm{z}_{obs})\subsetneq\widetilde{\Gamma}_{\alpha}(\bm{z}_{obs})$.
The following is an example in which $\Gamma_{\alpha}(\bm{z}_{obs})$ is
strictly smaller than $\widetilde{\Gamma}_{\alpha}(\bm{z}_{obs})$.

######  Example 3.

Suppose $\bm{y}_{obs}=(y_{1},\ldots,y_{n})^{\top}$ is a set of sample from the
model $Y_{i}=\theta_{0}+U_{i}$, $U_{i}\sim U(-1,1)$, $i=1,\ldots,n$. Here,
$\Theta=(-\infty,+\infty)$, ${\cal Y}=(-\infty,\infty)$. The point estimator
$\bar{Y}=\frac{1}{n}\sum_{i=1}^{n}Y_{i}$ is unbiased (and also
$\sqrt{n}$-consistent). Let’s use the test statistic
$\widetilde{T}(\bm{y}_{obs},\theta)=\bar{y}-\theta$. Since
$n\\{(\bar{Y}-\theta)+1\\}/2=\sum_{i=1}^{n}\frac{U_{i}+1}{2}$ follows a Irwin-
Hall distribution when $\theta$ is the true value $\theta_{0}$, by the
classical testing method a level-$95\%$ confidence interval is
$\widetilde{\Gamma}_{.95}(\bm{y}_{obs})=(\bar{y}-\frac{2}{n}q_{.975}+1,\bar{y}+\frac{2}{n}q_{.975}-1)$.
Here, $q_{.975}$ is the $97.5\%$ quantile of the Irwin-Hall distribution. If
we use the repro samples method, our confidence set is then
$\Gamma_{.95}(\bm{y}_{obs})=\big{\\{}\theta:y_{i}=\theta+u_{i}^{*},i=1,\ldots,n;\,\bar{y}\in(\theta-\frac{2}{n}q_{.975}+1,\theta+\frac{2}{n}q_{.975}-1);\exists\bm{u}^{*}\in(-1,1)^{n}\big{\\}}=\left\\{\theta:\theta\in\cap_{i=1}^{n}(y_{i}-1,y_{i}+1)\right\\}\cap\widetilde{\Gamma}_{.95}(\bm{y}_{obs})$.
Since the constraint $\\{\theta\in\cap_{i=1}^{n}(y_{i}-1,y_{i}+1)\\}$ is non-
trivial, the strictly smaller statement
$\Gamma(\bm{y}_{obs})\subsetneq\widetilde{\Gamma}(\bm{y}_{obs})$ holds for
non-trivially many realizations of $\bm{y}_{obs}$.

As an illustration, consider the case $\theta=0$ and $n=3$, in which
$q_{.975}=2.4687$. Suppose we have a sample realization
$\bm{y}_{obs}=(-0.430,0.049,0.371)^{\top}$. By a direct calculation,
$\widetilde{\Gamma}_{\alpha}(\bm{z}_{obs})=(-.6458,.6458)$ and
$\Gamma_{\alpha}(\bm{y}_{obs})=\big{\\{}\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits:y_{obs,i}=\theta+u_{i}^{*},i=1,2,3;\,\exists(u_{1}^{*},u_{2}^{*},u_{3}^{*})\in(-1,1)\times(-1,1)\times(-1,1)\big{\\}}\cap\widetilde{\Gamma}_{\alpha}(\bm{y}_{obs})=\big{\\{}(-1.430,.570)\cap(-.951,1.049)\cap(-.629,1.371)\big{\\}}\cap(-.6458,.6458)=(-.629,.570)$.
Clearly,
$\Gamma_{\alpha}(\bm{z}_{obs})\subsetneq\widetilde{\Gamma}_{\alpha}(\bm{z}_{obs})$.
We repeated the experiment $1,000$ times. The coverage rates of
$\widetilde{\Gamma}_{\alpha}(\bm{y}_{obs})$ and
$\Gamma_{\alpha}(\bm{y}_{obs})$ are both $95.1\%$. Among the $1,000$
repetitions, $219$ times
$\Gamma_{\alpha}(\bm{y}_{obs})=\widetilde{\Gamma}_{\alpha}(\bm{y}_{obs})$ and
$781$ times
$\Gamma_{\alpha}(\bm{y}_{obs})\subsetneq\widetilde{\Gamma}_{\alpha}(\bm{y}_{obs})$.
The average lengths of $\Gamma_{\alpha}(\bm{y}_{obs})$ and
$\widetilde{\Gamma}_{\alpha}(\bm{y}_{obs})$ are $0.915$ and $1.292$,
respectively.

As mentioned in Remark 2.2, a nuclear mapping function does not need to be a
test statistic. However, if it is defined through a test statistics by (12),
there are some immediate implications from Lemma 1. First, the confidence set
by the repro samples method is either smaller or the same as the one obtained
using the test statistic by the classical testing approach. So the set
$\Gamma(\bm{z}_{obs})$ by the repro samples method is more desirable. Second,
if the test statistic is optimal (in the sense that it leads to a powerful
test), then the confidence set constructed by the corresponding nuclear
mapping function is also optimal. In particular, we have the following
corollary. Here, a level-$\alpha$ uniformly most accurate (UMA) confidence set
(also called Neyman shortest) refers to a level-$\alpha$ confidence set that
minimizes the probability of false coverage (i.e., probability of covering a
wrong parameter value) over a class of level-$\alpha$ confidence sets; cf.,
Casella and Berger (1990, §9.3.2). See Appendix I for a proof of the
corollary.

######  Corollary 2.

(a) If the test statistic
$\widetilde{T}(\bm{z}_{obs},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$
corresponds to the uniformly most powerful test and
$\widetilde{\Gamma}(\bm{z}_{obs})$ is a level-$\alpha$ UMA confidence set,
then the set $\Gamma(\bm{z}_{obs})$ by the corresponding repro sample method
is also a level-$\alpha$ UMA confidence set.
(b) If the test statistic
$\widetilde{T}(\bm{z}_{obs},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$
corresponds to the uniformly most powerful unbiased test and
$\widetilde{\Gamma}(\bm{z}_{obs})$ is a level-$\alpha$ UMA unbiased confidence
set, then the confidence set $\Gamma(\bm{z}_{obs})$ by the corresponding repro
sample method is also a level-$\alpha$ UMA unbiased confidence set.

Finally, we can define the nuclear mapping through a pivot or an approximately
pivot quantity, which is often used as a test statistic. The probability
distribution of a pivot quantity does not depend on the unknown parameters.
Hence it is convenient to define the nuclear mapping function through a pivot
or an approximately pivot quantity, i.e., set
$T(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)=\widetilde{T}(\bm{z}^{\prime},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$,
where
$\widetilde{T}(\bm{z}^{\prime},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$
is a pivot quantity whose distribution satisfies
$\mathbb{P}\\{\widetilde{T}(\bm{Z}^{\prime},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\in
B\\}=\int_{{\bf t}\in B}f({\bf t})d{\bf t}$ or an approximately pivot quantity
whose distribution satisfies
$\mathbb{P}\\{\widetilde{T}(\bm{Z}^{\prime},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\in
B\\}=\int_{{\bf t}\in B}f({\bf t})d{\bf t}\,\\{1+o(\delta^{{}^{\prime}})\\}.$
Here, $f({\bf t})$ is a density function that is free of the unknown parameter
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$ and
$B\subset\mathbb{R}^{d}$ is any Borel set, and $\delta^{{}^{\prime}}>0$ is a
small number as defined in Theorems 1 and 2. In this case, we have simplified
in (4) that the set
$B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\equiv
B_{\alpha}$ is free or approximately free of
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$.

Let’s consider the example of a likelihood ratio test (LRT) statistic, which
is often regarded as the most commonly used asymptotic pivot quantity; cf.,
Reid and Cox (2015). Example 4 below compares the repro sample method works
with the LRT test.

We both consider an example where the LRT is asymptotically chi-square
distributed, and revisit Example 3 where the regular asymptotic result does
not hold. It demonstrates the flexibility of the repro samples method and its
ability to achieve improved performance over the commonly used likelihood
inference especially in finite sample settings.

######  Example 4 (Likelihood inference).

(a) Assume $L(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits|\bm{Z})$
is the likelihood function
$L(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits|\bm{Z})$
corresponding to model (2). The LRT statistic for test
$H_{0}:\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}=\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$
vs
$H_{1}:\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}\not=\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$
is
$\lambda_{\theta}(\bm{z}_{obs})=\frac{L(\theta|\bm{z}_{obs})}{\sup_{\theta^{\prime}}L(\theta^{\prime}|\bm{z}_{obs})}$.
Suppose the nuclear mapping function is defined as
$T(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)=-2\log\\{\lambda_{\theta}(\bm{z})\\}=-2\log\\{\lambda_{\theta}(G_{z}(\bm{u},\theta))\\}$.
A level-$\alpha$ confidence set by the repro sample method is then
$\Gamma_{\alpha}(\bm{z}_{obs})=\big{\\{}\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits:\bm{z}_{obs}=G(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits,\bm{u}^{*}),-2\log\\{\lambda_{\theta}(\bm{z}_{obs})\\}\in
B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits),\exists\,\bm{u}^{*}\in{\cal
U}\big{\\}},$ where
$B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ is a set
that satisfies (4). This $\Gamma_{\alpha}(\bm{z}_{obs})$ is an exact
level-$\alpha$ confidence set, if
$B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ can be
obtained for each given
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$ or directly using
Algorithm 1. Frequently, as $n\to\infty$ and under some regularity conditions,
$-2\log\\{\lambda_{\theta}(\bm{Z})\\}$ is asymptotically
$\chi^{2}$-distributed with ${p}$ degrees of freedom (Casella and Berger,
1990, p381), thus asymptotically free of
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$. It follows that an
asymptotic level-$\alpha$ confidence set is
$\Gamma_{\alpha}(\bm{z}_{obs})=\big{\\{}\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits:\bm{z}_{obs}=G(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits,\bm{u}^{*}),-2\log\\{\lambda_{\theta}(\bm{z}_{obs})\\}\leq
q_{\chi^{2}_{p}}(\alpha),\exists\,\bm{u}^{*}\in{\cal
U}\big{\\}}\subseteq\widetilde{\Gamma}_{\alpha}(\bm{z}_{obs})\equiv\big{\\{}\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits:-2\log\\{\lambda_{\theta}(\bm{z}_{obs})\\}\leq
q_{\chi^{2}_{p}}(\alpha)\big{\\}},$ where $q_{\chi^{2}_{p}}(\alpha)$ is the
$\alpha$-th quantile of the $\chi^{2}_{p}$ distribution and
$\widetilde{\Gamma}_{\alpha}(\bm{z}_{obs})$ is the confidence set obtained by
the regular LRT method. By Corollary 2, $\Gamma_{\alpha}(\bm{z}_{obs})$
preserves the optimality of $\widetilde{\Gamma}_{\alpha}(\bm{z}_{obs})$ when
$\widetilde{\Gamma}_{\alpha}(\bm{z}_{obs})$ is optimal.

(b) We now revisit Example 3, where the likelihood function
$L(\theta|\bm{y}^{obs})=\prod\limits_{i=1}^{n}\big{[}I\\{-1<y_{i}-\theta\leq
1\\}/2\big{]}\propto I\big{\\{}\max\limits_{1\leq i\leq
n}|y_{i}-\theta|<1\big{\\}}.$ The LRT statistic is
$\lambda_{\theta}(\bm{z}_{obs})=I\big{\\{}\max\limits_{1\leq i\leq
n}|y_{i}-\theta|<1\big{\\}}$ and $-2\log\\{\lambda_{\theta}(\bm{Y})\\}$ does
not converge to a $\chi^{2}$ distribution. However,
$\lambda_{\theta}(\bm{y}_{obs})$ is decreasing in
$\widetilde{T}(\bm{y}_{obs},\theta)=\max\limits_{1\leq i\leq
n}|y_{i}-\theta|$. So, the LRT rejects $H_{0}$ with a large
$\widetilde{T}(\bm{y}_{obs},\theta)$, which leads to a level-$\alpha$
confidence set
$\widetilde{\Gamma}_{\alpha}(\bm{y}_{obs})=\\{\theta:\max\limits_{1\leq i\leq
n}|y_{i}-\theta|<\alpha^{1/n}\\}$. Although we could directly define our
nuclear mapping function through
$\widetilde{T}(\bm{y}_{obs},\theta)=\max\limits_{1\leq i\leq n}|y_{i}-\theta|$
to obtain a repro sample confidence interval, we consider here an alternative
choice of a vector nuclear mapping function
$T(\bm{u},\theta)=(u_{(1)},u_{(n)})$, where $u_{(k)}$ is the $k$th order
statistic of a sample from $U(-1,1)$. Let $c_{\alpha}\in(0,1)$ be a solution
of $({c_{\alpha}+1})^{n}-2^{n-1}c_{\alpha}^{n}=2^{n-1}{(1-\alpha)}$. We can
show $\mathbb{P}\\{(U_{(1)},U_{(n)})\in B_{\alpha}\\}=\alpha$ for
$B_{\alpha}=(-1,-c_{\alpha})\times(c_{\alpha},1)$. Thus, a level-$\alpha$
confidence set is
$\Gamma_{\alpha}(\bm{y}_{obs})=\\{\theta:y_{(n)}-1<\theta<y_{(1)}+1;$
$y_{(1)}-\theta<-c_{\alpha},y_{(n)}-\theta>c_{\alpha}\\}=\big{(}\max\big{\\{}y_{(1)}+c_{\alpha},y_{(n)}-1\big{\\}},\min\big{\\{}y_{(n)}-c_{\alpha},y_{(1)}+1\big{\\}}\big{)}$.
We have conducted a numerical study with true $\theta_{0}=0$ and $n=5,20,200$,
respectively. In all cases, the coverage rates of both confidence intervals
are right on target around $95\%$ in 3000 repetitions. The intervals by the
repro sample method are consistently shorter than those obtained using the LRT
method across all sample sizes, as evidenced in Figure 2. Although a LRT is
uniformly most powerful for a simple-versus-simple test by the Neyman-Pearson
Lemma, it is not the case for the two-sided test in this example. we are able
to explore and use the repro samples method to obtain a better confidence
interval.

Figure 2: The side-by-side box plots comparing the widths of the two 95%
confidence intervals discussed in Example 4, for $n=5,20$ and $200$,
respectively; Number of repetitions $=3000$.

###  3.2  Define a nuclear mapping and control nuisance parameters

Frequently researchers are interested in making inference only for a certain
parameter, say $\eta$, and the remaining parameters
$\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits$ in
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits=(\eta,\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits^{\top})^{\top}$
are considered nuisance parameters. There is a large number of literature on
this topic where profiling and conditional inference remain to be two major
techniques that handle nuisance parameters (Cox and Reid, 1987; Davison, 2003,
and references therein). These techniques typically lead to test statistics,
say $\widetilde{T}(\bm{z},\eta)$, that do not involve the nuisance parameter
$\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits$ in expressions. Often times
however the distributions of these test statistics
$\widetilde{T}(\bm{z},\eta)=\widetilde{T}(G_{z}((\eta,\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits^{\top})^{\top},\bm{u}),\eta)$
still depend on the unknown nuisance parameter
$\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits$ and thus not readily
available. To obtain the distributions of $\widetilde{T}(\bm{z},\eta)$,
researchers often invoke large samples and assume that there exists a
consistent plug-in estimator of the nuisance parameter; e.g., Reid and Cox
(2015); Chuang and Lai (2000). In this subsection, however, we use a similar
profiling idea to provide an alternative approach. The advantage of our
approach is that it can give us a confidence set with a guaranteed finite-
sample coverage rate.

Suppose $T_{a}(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$
is a mapping function from ${\cal U}\times\Theta\to\mathbb{R}^{d}$. For
instance, we can often make
$T_{a}(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)=\widetilde{T}(\bm{z},\eta)$,
a test statistic mentioned above using either a profile or a conditional
approach. Or, more generally
$T_{a}(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ can be
any nuclear mapping function
$T(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ used to
make a joint inference about
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits=(\eta,\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits^{\top})^{\top}$.
Suppose we have an artificial data point
$\bm{z}^{\prime}=G_{z}(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$
and let $\Gamma_{\alpha^{\prime}}(\bm{z}^{\prime})$ be a
level-$\alpha^{\prime}$ confidence set for
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$, where
$\Gamma_{\alpha^{\prime}}(\bm{z})$ has the expression of (5) or (10) but with
$\bm{z}_{obs}$ replaced by the artificial data
$\bm{z}^{\prime}=G_{z}(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$
and the nuclear mapping function replaced by
$T_{a}(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$. Now,
beside
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits=(\eta,\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits^{\top})^{\top}$,
let us consider another copy of parameter vector
$\widetilde{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}}\nolimits=(\eta,\widetilde{\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}}\nolimits^{\top})^{\top}\in\Theta$
that has the same $\eta$ but perhaps a different
$\widetilde{\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}}\nolimits\not=\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits$.
To explore the impact of a potentially wrong nuisance parameter
$\widetilde{\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}}\nolimits\not=\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits$,
we define

$\nu(\bm{u};\eta,\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits,\widetilde{\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}}\nolimits)=\inf_{\alpha^{\prime}}\\{\alpha^{\prime}:\widetilde{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}}\nolimits\in\Gamma_{\alpha^{\prime}}(\bm{z})\\}=\inf_{\alpha^{\prime}}\\{\alpha^{\prime}:\widetilde{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}}\nolimits\in\Gamma_{\alpha^{\prime}}(G_{z}(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits))\\}.$
(13)

From Corollary 1, we can interpret
$1-\nu(\bm{u};\eta,\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits,\widetilde{\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}}\nolimits)$
as a ‘$p$-value’ when we test the null hypothesis that the data $\bm{z}$ is
from
$\widetilde{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}}\nolimits=(\eta,\widetilde{\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}}\nolimits^{\top})^{\top}$
but in fact
$\bm{z}=G_{z}(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$
with
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits=(\eta,\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits^{\top})^{\top}$.
Our proposed profile nuclear mapping function for the target parameter $\eta$
is defined as

$\displaystyle
T^{P}(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)=\min_{\widetilde{\xi}}\,\,\nu(\bm{u};\eta,\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits,\widetilde{\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}}\nolimits).$
(14)

Although
$T^{P}(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ still
depends on
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits=(\eta,\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits^{\top})^{\top}$,
including the unknown nuisance parameter
$\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits$, it is dominated by a
$U(0,1)$ random variable as stated in the lemma below. See Appendix I for a
proof of the lemma.

######  Lemma 2.

Suppose $T^{P}(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$
is a profile nuclear mapping function from ${\cal U}\times\Theta\to[0,1]$
defined above. Then, we have
$\mathbb{P}\left\\{T^{P}(\bm{U},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\leq\alpha\right\\}\geq\alpha.$

Based on the above lemma and following the construction outlined in Section 1,
we define

$\displaystyle\Xi_{\alpha}(\bm{z}_{obs})=\bigg{\\{}\eta:\exists\,\bm{u}^{*}\in{\cal
U}\,\text{and}\,\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits,$
$\displaystyle\qquad\text{ s.t. }\begin{pmatrix}\eta\\\
\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits\end{pmatrix}\in\Theta,\bm{z}_{obs}=G_{z}\left(\begin{pmatrix}\eta\\\
\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits\end{pmatrix},\bm{u}^{*}\right),\,T^{P}\left(\bm{u}^{*},\begin{pmatrix}\eta\\\
\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits\end{pmatrix}\right)\leq\alpha\bigg{\\}}.$
(15)

The next theorem states that $\Xi_{\alpha}(\bm{z}_{obs})$ is a level-$\alpha$
confidence set of the target parameter $\eta$. A proof is provided in
Appendix.

######  Theorem 3.

Suppose that
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}=(\eta_{0},\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits_{0}^{\top})^{\top}$
is the true parameter with
$\bm{Z}=G_{z}(\bm{U},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0})$
and $\Xi_{\alpha}(\cdot)$ is defined in (3.2), we have
$\mathbb{P}\left\\{\eta_{0}\in\Xi_{\alpha}(\bm{Z})\right\\}\geq\alpha$.

In cases when the distribution of
$T_{a}(\bm{U},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ is
known, the function
$\nu(\cdot;\eta,\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits,\widetilde{\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}}\nolimits)$
in (13) and the nuclear mapping function
$T^{P}(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ have
explicit formulas. When the distribution of
$T_{a}(\bm{U},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ is not
available, we can use a Monte-Carlo method to get
$\nu(\cdot;\eta,\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits,\widetilde{\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}}\nolimits)$
for each given
$(\eta,\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits,\widetilde{\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}}\nolimits)$
as described in the lemma below. See Appendix I for a proof.

######  Lemma 3.

Let $D_{{\cal S}_{\theta}}(\bm{t})$ be the empirical depth function that is
computed based on the Monte-Carlo points in ${\cal S}_{\theta}$, where ${{\cal
S}_{\theta}}=\\{T_{a}(\bm{u}^{s},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits),\bm{u}^{s}\in{\cal
V}\\}$ is the Monte-Carlo set defined in Algorithm 1 for a given
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits=(\eta,\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits^{\top})^{\top}$.
Let $F_{\mathcal{V}|D}(s)=\frac{1}{|{\cal V}|}\sum_{\bm{u}^{s}\in{\cal
V}}I\big{\\{}D_{{\cal
S}_{\theta}}\big{(}T_{a}(\bm{u}^{s},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\big{)}\leq
s\big{\\}}$ be the empirical CDF of $D_{{\cal
S}_{\theta}}\big{(}T_{a}(\bm{u}^{s},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\big{)},$
then

$\nu(\bm{u};\eta,\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits,\widetilde{\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}}\nolimits)=\inf_{\\{\bm{u}^{*}:G_{z}(\bm{u}^{*},\tilde{\theta})=G_{z}(\bm{u},\tilde{\theta})\\}}\left\\{1-F_{\mathcal{V}|D}\left(D_{{\cal
S}_{\tilde{\theta}}}\big{(}T_{a}(\bm{u}^{*},\tilde{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}}\nolimits)\big{)}\right)\right\\},$

where
$\tilde{\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}}\nolimits=(\eta,\widetilde{\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}}\nolimits^{\top})^{\top}.$

Using Lemma 3, we can evaluate
$T^{P}(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ via
(14), and then obtain the confidence set in (3.2). Note that, for any given
$(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$, we can
express
$\nu(\bm{u};\eta,\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits,\widetilde{\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}}\nolimits)$
as a function of
$\widetilde{\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}}\nolimits$ only, so we can
use a build-in optimization function in R (or another computing package) to
obtain $T^{P}(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$
in (14). The following continuation of Example 2 illustrates constructing
confidence interval (3.2) for a parameter of interests in presence of nuisance
parameters. The case study example in Section 4 provides another illustration
of a more difficult inference problem in a normal mixture model where the
target parameter is the number of unknown components.

######  Example 2 (Continuation).

We continue Example 2(b) to obtain a confidence interval for a single
regression coefficient $\theta_{k}(\zeta)$ in the censored quantile regression
model (with the remaining coefficients as nuisance parameters),
$k=0,\ldots,p$. We consider
$T_{a}(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ be the
$(k+1)$-th element of the vector
$T(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)=\big{\\{}\sum\limits_{i=1}^{n}\bm{x}_{i}\bm{x}_{i}^{T}I({\bm{x}}_{i}^{T}\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits>0)\big{\\}}^{-1}\sum\limits_{i=1}^{n}{\bm{x}}_{i}I{({\bm{x}}_{i}^{T}\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits>0)I(u_{i}\leq\tau)}$
that is previously used in Example 2(b) to make joint inference on the entire
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$. Following Theorem 3
and Lemma 3, a confidence set for the single parameter $\theta_{k}(\zeta)$ is

$\displaystyle\Xi_{\alpha}(\bm{y}_{obs})=\left\\{\theta_{k}:\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits=(\theta_{k},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{(-k)}),g(\bm{y}_{obs},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits,\bm{u}^{*})=0,F_{\mathcal{V}|D}\big{(}T_{a}(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\big{)}\geq
1-\alpha,\exists\bm{u}^{*}\in{\cal U}\right\\},$

where $F_{\mathcal{V}|D}(s)$ is defined in Lemma 3. Table 3 reports the
performance of the above confidence interval for $\theta_{2}(\zeta)$ from
simulation studies on the same setting as in Example 2(B). We observe that for
$\zeta=0.5$, repro samples achieves the desirable coverage while the bootstrap
method by Bilias et al. (2000) consistently under covers. The widths of the
confidence intervals from the two approaches are similar. For $\zeta$ =0.9,
both approaches well cover the true parameter, but the bootstrap confidence
intervals appear to be much wider.

Quantile | Error Distribution | repro samples | Bootstrap
---|---|---|---
Coverage | Width | Coverage | Width
$\zeta=0.50$ | Normal | 0.944(0.007) | 0.747(0.006) | 0.920(0.009) | 0.727(0.007)
Mixture | 0.950(0.007) | 0.832(0.007) | 0.935(0.008) | 0.839(0.008)
Heteroscedastic | 0.940(0.008) | 0.826(0.007) | 0.921(0.009) | 0.844(0.009)
$\zeta=0.90$ | Normal | 0.979(0.005) | 0.813(0.007) | 0.983(0.004) | 4.190 (0.424)
Mixture | 0.973(0.005) | 1.407(0.013) | 0.973(0.005) | 3.260(0.312)
Heteroscedastic | 0.945(0.007) | 1.149(0.009) | 0.967(0.006) | 4.124(0.439)

Table 3:  The model and error distributions are identical to those in Table 2.
Here we compare the coverage rates and widths of confidence intervals for the
regression coefficient of $x_{i2}$ by the repro samples approaches and the
bootstrap approach (Bilias et al., 2000). The coverage rates and widths are
computed using 1000 repetitions.

###  3.3 Finding candidate set for target parameter

Sometimes when an explicit expression of the confidence set
$\Gamma_{\alpha}(\bm{z}_{obs})$ is not available, one can use Algorithm 1 to
search through the parameter space $\Theta$ and construct
$\Gamma_{\alpha}(\bm{z}_{obs})$. Finding a suitable candidate set to limit the
search can greatly improve the computing efficiency of the algorithm. Here, a
candidate set, say $\widehat{\Theta}=\widehat{\Theta}(\bm{z}_{obs})$ refers to
a subset of the parameter space $\Theta$ that is pre-screened using the
observed data $\bm{z}_{obs}$. In this section, we first provide a general
result of the impact of using a data-dependent candidate set on the coverage.
We then discuss the special case of a discrete target parameter in which we
can take advantage of a many-to-one mapping that is unique in a repro samples
setup to create an effective pre-screen candidate set to significantly reduce
our computing cost.

The candidate set $\widehat{\Theta}=\widehat{\Theta}(\bm{z}_{obs})$ may be
obtained using different pre-screening methods. The key is that we need to
overcome the issue of a “double usage” of the data but still ensure an at
least approximate coverage of the final output confidence set

$\displaystyle\Gamma_{\alpha}^{\prime}(\bm{z}_{obs})=\widehat{\Theta}(\bm{z}_{obs})\cap\Gamma_{\alpha}(\bm{z}_{obs})$
$\displaystyle\qquad=\big{\\{}\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits:\exists\,\bm{u}^{*}\in{\cal
U},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits\in\widehat{\Theta},\,\hbox{s.t.}\,\bm{z}_{obs}=G_{z}({\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits},\bm{u}^{*}),\,T(\bm{u}^{*},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\in
B_{\alpha}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\big{\\}},$
(16)

we impose a condition that

$\mathbb{P}\big{\\{}\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}\in\widehat{\Theta}(\bm{Z})\big{\\}}=1-o(\delta^{\prime}),$
(17)

where $\delta^{\prime}>0$ is a small number as described in Theorem 1. Here,
the “double usage” refers to the fact that we use the observed data in
constructing both the candidate set $\widehat{\Theta}(\bm{z}_{obs})$ and
confidence set $\Gamma_{\alpha}(\bm{z}_{obs})$. The corollary below contains a
result showing that the output set $\Gamma_{\alpha}^{\prime}(\bm{z}_{obs})$ is
still an (approximate) level-$\alpha$ confidence set. See Appendix I for a
proof of the corollary.

######  Corollary 3.

Let
$\bm{Z}=G_{z}(\bm{U},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0})$
and $\widehat{\Theta}=\widehat{\Theta}(\bm{Z})$ is a candidate set that
satisfies (17). Then, we have
$\mathbb{P}\big{\\{}\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}\in\Gamma_{\alpha}^{\prime}(\bm{Z})\big{\\}}\geq\alpha-o(\delta^{\prime}).$

When the parameter space $\Theta$ is continuous and the mapping function
$T(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ is
continuous in $\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits$, we may
use a grid search method to cover the entire $\Theta$, as discussed in Remark
2.6. A data-dependent candidate set $\widehat{\Theta}$ can help further narrow
down to the grid search; cf., Example Example 2 (Continuation) and also
Michael et al. (2019).

An intriguing yet unique situation in the repro samples development is when
the parameter space $\Theta$ is discrete. We can take advantage of a “many-to-
one” mapping that often exists in the setup to get a data-dependent candidate
set $\widehat{\Theta}$ that is much smaller than $\Theta$. We describe the
above key idea and a generic approach below. Section 4 includes a case study
with with greater details to provide a specific illustration on how to obtain
this kind of candidate set $\widehat{\Theta}.$

For convenience of discussion, we rewrite (2) as
$\bm{Z}=G_{z}(\kappa_{0},\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits_{0},\bm{U})$,
where $\kappa_{0}$ is the true value of a discrete parameter with a parameter
space ${\it\Upsilon}$ and, if applicable,
$\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits_{0}$ are the values of other
(nuisance) parameters. Write
$\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits=(\kappa,\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits^{\top})^{\top}$.
Our goal is to find a candidate set
$\widehat{\it\Upsilon}=\widehat{\it\Upsilon}(\bm{z}_{obs})$ such that
$\widehat{\it\Upsilon}(\bm{z}_{obs})$ is significantly smaller than
${\it\Upsilon}$ and $\widehat{\it\Upsilon}$ satisfies (17).

Let us first examine a “many-to-one” mapping often inherited in a repro sample
setup involving a discrete parameter $\kappa$. In particular, based on the
equation
$\bm{z}_{obs}=G_{z}(\kappa_{0},\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits_{0},\bm{u}^{rel}),$
we can rewrite $\kappa_{0}$ as a solution of an optimization problem
$\kappa_{0}=\arg\min\limits_{\kappa}\min\limits_{\xi}\|\bm{z}_{obs}-G_{z}(\kappa,\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits,\bm{u}^{rel})\|$
or, more generally,

$\kappa_{0}=\arg\min\limits_{\kappa}\min\limits_{\xi}L\big{(}\bm{z}_{obs},G_{z}(\kappa,\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits,\bm{u}^{rel})\big{)},$
(18)

where $L\big{(}\bm{z},\bm{z}^{\prime}\big{)}$ is a continuous loss function
that measures the difference between two copies of data $\bm{z}$ and
$\bm{z}^{\prime}$. However, we do not know (observe) $\bm{u}^{rel}$. By
replacing $\bm{u}^{rel}$ with a repro copy $\bm{u}^{*}$ in (18), we get

$\kappa^{*}=\arg\min\limits_{\kappa}\min\limits_{\xi}L\big{(}\bm{z}_{obs},G_{z}(\kappa,\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits,\bm{u}^{*})\big{)},$
(19)

which maps a value $\bm{u}^{*}\in{\cal U}$ to a $\kappa^{*}\in{\it\Upsilon}$.
Typically, ${\cal U}$ is uncountable and ${\it\Upsilon}$ is countable, so the
mapping in (19) is “many-to-one”; That is, many $\bm{u}^{*}\in{\cal U}$
correspond to one $\kappa^{*}\in{\it\Upsilon}$. We are particularly interested
in the subset
$S=\big{\\{}\bm{u}^{*}:\kappa_{0}=\arg\min\limits_{\kappa}\min\limits_{\xi}L\big{(}\bm{z}_{obs},G_{z}(\kappa,\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits,\bm{u}^{*})\big{)}\big{\\}}$
$\subset{\cal U}$, where for any $\bm{u}^{*}\in S$ the mapping in (19)
produces $\kappa^{*}=\kappa_{0}$. Since $\bm{u}^{rel}\in S$,
$S\not=\emptyset$. We assume $S\subset{\cal U}$ is a nontrivial set containing
many elements that can not be ignored.

Now suppose $\mathbb{P}\big{\\{}\bm{U}\in S\big{|}\bm{z}_{obs}\big{\\}}>0$ has
a nontrivial probability. Then, we can use a Monte-Carlo method to obtain a
candidate set $\widehat{\Upsilon}=\widehat{\Upsilon}(\bm{z}_{obs})$ such that
$\mathbb{P}\big{\\{}\kappa_{0}\not\in\widehat{\it\Upsilon}(\bm{Z})\big{\\}}=o(\delta^{\prime})$.
The idea is to simulate a sequence of $\bm{u}^{c}\sim\bm{U}$, say, ${\cal
V}_{c}=\big{\\{}\bm{u}_{1}^{c},\ldots,\bm{u}_{N}^{c}\big{\\}}$. When $N=|{\cal
V}_{c}|$ is large enough, ${\cal V}_{c}\cap S\not=\emptyset$ with a high
(Monte-Carlo) probability. That is $\kappa_{0}\in\widehat{\it\Upsilon}_{{\cal
V}_{c}}(\bm{z}_{obs})$ with a high (Monte-Carlo) probability, where
$\widehat{\it\Upsilon}_{{\cal V}_{c}}(\bm{z}_{obs})$ is our candidate set

$\widehat{\it\Upsilon}_{{\cal
V}_{c}}(\bm{z}_{obs})=\left\\{\kappa^{*}=\arg\min\limits_{\kappa}\min\limits_{\xi}L\big{(}\bm{z}_{obs},G_{z}(\kappa,\mathop{\xi\kern-4.49997pt\hbox{$\xi$}}\nolimits,\bm{u}^{c})\big{)}\,\,\big{|}\,\,\bm{u}^{c}\in{\cal
V}_{c}\right\\}.$ (20)

We show next that
$\mathbb{P}\big{\\{}\kappa_{0}\in\widehat{\it\Upsilon}_{{\cal
V}_{c}}(\bm{Z})\big{\\}})=1-o(\delta^{\prime})$, where
$\delta^{\prime}=\delta^{\prime}({\footnotesize|{\cal V}_{c}|})\to 0$ as
$|{\cal V}_{c}|\to\infty$.

Let $\bm{U}^{*}$ be an independent copy of $\bm{U}$. Formally, we assume the
following condition:

  * (N1)

For any $\bm{u}^{rel}\in\mathcal{U}$, there exists a neighborhood of
$\bm{u}^{rel}$, say $S_{\cal N}(\bm{u}^{rel})$, such that

$S_{\cal N}(\bm{u}^{rel})\subseteq
S=\big{\\{}\bm{u}^{*}:\kappa_{0}=\arg\min\limits_{\kappa}\min\limits_{\theta}L\big{(}G_{z}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0},\kappa_{0},\bm{u}^{rel}),G_{z}(\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits,\kappa,\bm{u}^{*})\big{)}\big{\\}}$

(21)

and $\mathbb{P}_{(\bm{U}^{*},\bm{U})}\left\\{\bm{U}^{*}\in S_{\cal
N}(\bm{U})\right\\}\geq P_{\cal N}>0$, for a positive number $P_{\cal N}$.

Although the neighborhood $S_{\cal N}$ and the number $P_{\cal N}$ are problem
specific, Condition (N1) is typically satisfied (or nearly satisfied) when
$L(\cdot,\cdot)$ is a continuous loss function and $|{\cal U}|$ is
uncountable. The following lemma suggests that, under Condition (N1) and as
$|{\cal V}_{c}|\to 0$, $\widehat{\it\Upsilon}_{{\cal V}_{c}}(\bm{z}_{obs})$ in
(20) satisfy condition (17). Here, the probability $\mathbb{P}_{\bm{U},{\cal
V}_{c}}(\cdot)$ refers to the joint distribution of $\bm{U}$ and $|{\cal
V}_{c}|$ Monte-Carlo $\bm{U}$’s in ${\cal V}_{c}$. See Appendix I for a proof
of the lemma.

######  Lemma 4.

Suppose Condition (N1) holds and $\widehat{\it\Upsilon}(\bm{z}_{obs})$ is
defined in (20). Then,

$\displaystyle\mathbb{P}_{\bm{U},{\cal
V}_{c}}\left\\{\kappa_{0}\not\in\widehat{\it\Upsilon}_{{\cal
V}_{c}}(\bm{Z})\right\\}\leq(1-P_{\cal N})^{|{\cal V}_{c}|}.$

Based on the lemma, we propose the following theorem. Here, the probability
$\mathbb{P}_{\bm{U}|{\cal V}_{c}}(\cdot)$ refers to the condition distribution
of $\bm{U}$ given the set ${\cal V}_{c}$. See Appendix I for a proof.

######  Theorem 4.

Suppose Condition (N1) holds and $\widehat{\it\Upsilon}(\bm{z}_{obs})$ is
defined in (20). Then, there exists a positive number $c$ such that we have,
in probability, (a) $\mathbb{P}_{\bm{U}|{\cal
V}_{c}}\left\\{\kappa_{0}\in\widehat{\it\Upsilon}_{{\cal
V}_{c}}(\bm{Z})\right\\}=1-o_{p}\left(e^{-c|{\cal V}_{c}|}\right)$ and (b)
$\mathbb{P}_{\bm{U}|{\cal
V}_{c}}\big{\\{}\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits_{0}\in\Gamma_{\alpha}^{\prime}(\bm{Z})\big{\\}}\geq\alpha-
o_{p}\left(e^{-c|{\cal V}_{c}|}\right)$.

In the theorem, the small $o_{p}\left(e^{-c|{\cal V}_{c}|}\right)$ term tends
to $0$ as the size of the Monte-Carlo set $|{\cal V}_{c}|\to\infty$, and the
result holds for any fixed (data) sample size $n<\infty$. So the theorem is a
small sample result. Also, the “in probability” statement regards the Monte-
Carlo randomness of ${\cal V}_{c}$. Based on Theorem 4, we have a clean
interpretation for the frequentist coverage statement: if we use the proposed
Monte-Carlo method to obtain the candidate set $\widehat{\it\Upsilon}_{{\cal
V}_{c}}(\bm{Z})$, we have a near $100\%$ chance (i.e., in probability
statement in ${\cal V}_{c}$) that the coverage provability of the obtained set
$\Gamma_{\alpha}^{\prime}(\bm{Z})$ is no less than $\alpha$ (i.e., frequency
statement in $\bm{U}$).

We will further discuss and provide an extension of Condition (N1) and Theorem
4 in Section 4.3.

##  4 Case study example: normal mixture with unknown number of components

Mixture distributions are a commonly used model in data science. However, a
challenging task remains to be quantifying the uncertainty and making
inference for the unknown number of components, say $\tau_{0}$, in a finite
mixture distribution. Although several existing procedures provide point
estimates for $\tau_{0}$, the question of how to assess the uncertainty of
estimation of $\tau_{0}$ is still open. This difficulty is due in large part
to the fact that the unknown $\tau_{0}$ is a discrete natural number and
typical tools such as the large sample CLT are no long applicable to the point
estimators of $\tau_{0}$. In this section, we demonstrate how to use the repro
samples method to solve this open question in a normal mixture model by
constructing a finite sample level-$\alpha$ confidence set for the unknown
$\tau_{0}$. This case study example also sheds light on how to address
inference problems of discrete parameters and handle nuisance parameters in
general cases.

###  4.1 Notations: normal mixture model

Suppose $Y_{i}\sim N\big{(}\mu^{(0)}_{k},(\sigma^{(0)}_{k})^{2}\big{)}$, one
of the $\tau_{0}$ models in ${\cal
M}=\\{N\big{(}\mu^{(0)}_{k},(\sigma^{(0)}_{k})^{2}\big{)},$
$k\in\\{1,\ldots,\tau_{0}\\}\\}$, for $i=1,\ldots,n.$ Write
${\bm{Y}}=(Y_{1},\ldots,Y_{n})^{\top}$. Here, the model parameters are
${\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits}=(\tau,\bm{M}_{\tau},{\bm{\mu}}_{\tau},\bm{\sigma}_{\tau}^{2})^{\top}$,
where ${\bm{M}}_{\tau}=({\bm{m}}_{1},\ldots,{\bm{m}}_{n})^{\top}$ is the
unknown (unobserved) $n\times\tau$ membership assignment matrix, where
${\bm{m}}_{i}=(m_{i1},\ldots,m_{i\tau})^{\top}$ and $m_{ik}$ is the indicator
of whether $Y_{i}$ is from the $k$th distribution $N(\mu_{k},\sigma^{2})$,
$1\leq i\leq n$, $1\leq k\leq\tau$. We denote by the corresponding true values
$\tau_{0}$, $\bm{M}_{0}=({\bm{m}}^{(0)}_{1},\ldots,{\bm{m}}^{(0)}_{n})^{\top}$
$=$ $\big{(}(m^{(0}_{11},\ldots,m^{(0)}_{1\tau})^{\top},\dots,(m^{(0)}_{n1},$
$\ldots,$ $m^{(0)}_{n\tau})^{\top}\big{)}^{\top}$,
${\bm{\mu}}_{0}=(\mu_{1}^{(0)},\ldots,\mu_{{}_{\tau_{0}}}^{(0)})^{\top}$,
$\bm{\sigma}_{0}=(\sigma^{(0)}_{1},\dots,\sigma^{(0)}_{\tau_{0}})^{\top}$ and
${\bm{\theta}}_{0}$ $=$ $(\tau_{0},\bm{M}_{0},{\bm{\mu}}_{0},$
$\bm{\sigma}_{0})^{\top}$, respectively. Because we are interested in
inference of the unknown number of components $\tau<n,$
${\bm{\mu}}_{\tau}=(\mu_{1},\ldots,\mu_{\tau})^{\top}$ and
$\bm{\sigma}_{\tau}=(\sigma_{1},\dots,\sigma_{\tau})^{\top}$ are therefore
unknown nuisance parameters of the $\tau$ components’ means and variances.

We can re-express
$Y_{i}=\sum_{k=1}^{\tau}m^{(0)}_{ik}(\mu_{k}^{(0)}+\sigma^{(0)}_{k}U_{i})$,
where $U_{i}\sim N(0,1)$, in a vector form,

$\displaystyle{\bm{Y}}={\bm{M}}_{0}{\bm{\mu}}_{0}+\hbox{diag}(\bm{U})\bm{M}_{0}{\bm{\sigma}_{0}}={\bm{M}}_{0}{\bm{\mu}}_{0}+\hbox{diag}(\bm{M}_{0}{\bm{\sigma}_{0}})\bm{U},$
(22)

where $\bm{U}=(U_{1},\ldots,U_{n})^{\top}$. The realized version of (22) is
${\bm{y}_{obs}}={\bm{M}}_{0}{\bm{\mu}}_{0}+\hbox{diag}({\bm{u}}^{rel}){\bm{M}}_{0}\bm{\sigma}_{0},$
where $\bm{M}_{0}$ is the realized membership assignment matrix that generates
$\bm{y}_{obs}$. A typical setup
$Y_{i}\sim\sum_{k=1}^{\tau}p_{k}N(\mu_{k},\sigma_{k}^{2})$ further assumes
that ${\bm{m}}_{i}=(m_{i1},\ldots,m_{i\tau})^{\top}$ is a random draw from
multinomial distribution $\text{Multinomial}(\tau;p_{1},\ldots,p_{\tau})$ with
unknown proportion parameters $(p_{1},\ldots,p_{\tau})$,
$\sum_{k=1}^{\tau}p_{k}=1$.

The literature confirms the identifiability issue in a mixture model that
multiple sets of parameters could possibly generate the same data (Liu and
Shao, 2003). To account for this identifiability issue, we re-define the true
parameters as those that correspond to the ones with the smallest $\tau$,
i.e.,

$\displaystyle(\tau_{0},\bm{M}_{0},{\bm{\mu}}_{0},\bm{\sigma}_{0}^{2})=\underset{\\{(\tau,\bm{M}_{\tau},{\bm{\mu}}_{\tau},\bm{\sigma}^{2})|{\bm{M}}_{\tau}{\bm{\mu}}_{\tau}={\bm{M}}_{0}{\bm{\mu}}_{0},\,{\bm{M}}_{\tau}{\bm{\sigma}}_{\tau}={\bm{M}}_{0}{\bm{\sigma}}_{0}\\}}{\arg\min}\tau,$
(23)

when the identifiability issue arises. Here, $\bm{M}_{0}$ and ${\bm{\mu}}_{0}$
on the right hand side are the true parameter values that generates
$\bm{y}^{obs}.$ With a slight abuse of notation, we still use the same
notations on the left hand side to denote the ones with the smallest $\tau$.
We further assume that
$(\tau_{0},\bm{M}_{0},{\bm{\mu}}_{0},\bm{\sigma}_{0}^{2})$ defined in (23) is
unique. We also assume that $0<\sigma_{\min}\leq\sigma_{\max}<\infty,$ where
$\sigma_{\min}=\|\bm{\sigma}_{0}\|_{0}$ and
$\sigma_{\max}=\|\bm{\sigma}_{0}\|_{\infty}.$ That is, among the $\tau_{0}$
components, none of the variance $\sigma_{k}^{2}$ can completely dominate
another $\sigma_{k^{\prime}}^{2}$, for any $k\not=k^{\prime}$. The inference
we develop is on the set of parameters defined in (23) with the smallest
$\tau_{0}$.

The inference problem for the target parameter $\tau_{0}$ is complicated.
Here, we have both types of nuisance parameters of varying sizes: discrete
nuisance parameters $\bm{M}_{\tau}$ that are potentially high-dimensional and
continuous nuisance parameters $(\bm{\mu}_{\tau},\bm{\sigma}_{\tau})$ whose
lengths change with $\tau$. To have an effective repro samples approach that
is computationally efficient, we handle the two types of nuisance parameters
in two different ways. For the discrete nuisance parameter $\bm{M}_{\tau}$ we
take advantage of the many-to-one mapping described in Section 3.3 to obtain a
manageable candidate set. For the continuous nuisance parameters
$(\bm{\mu}_{\tau},\bm{\sigma}_{\tau})$, we utilize a specific feature of the
normal mixture distribution, namely that, given $\tau$, we can find sufficient
statistics for $(\bm{\mu}_{\tau},\bm{\sigma}_{\tau})$ and avoid searching
through all possible values of $(\bm{\mu}_{\tau},\bm{\sigma}_{\tau})$. In
Section 4.2, we define a nuclear mapping function
$T^{P}(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ that
minimizes the impact of the nuisance parameters. In Section 4.3, we construct
an effective candidate set of $(\tau,\bm{M}_{\tau})$ over the corresponding
discrete parameter space to greatly improve our computing efficiency.

###  4.2 Nuclear mapping function for given $(\tau,\bm{M}_{\tau})$

In most derivations of this subsection, $\tau$ is fixed. For notation
simplicity, we suppress the subscripts in
$(\bm{\mu}_{\tau},\bm{\sigma}_{\tau})$ and simply write
$(\bm{\mu},\bm{\sigma})$. Note that the vectors $\bm{\mu}$ and $\bm{\sigma}$
have the length of $\tau$ that changes when $\tau$ varies. Also, let
$\bm{H}_{\tau}={\bm{M}}_{\tau}({\bm{M}}_{\tau}^{\top}{\bm{M}}_{\tau})^{-1}{\bm{M}}_{\tau}^{\top}$,
$\bm{H}_{\tau}^{(k)}=\bm{H}_{\tau}[\mathcal{I}_{k},\mathcal{I}_{k}]$ be the
submatrix of $\bm{H}_{\tau}$ corresponding to the $k$th component with
$\mathcal{I}_{k}=\\{i:m_{ik}=1\\}$ being the index set corresponding component
$k$, and $\bm{Y}^{(k)}=\bm{Y}[\mathcal{I}_{k}]$ be the subvector of $\bm{Y}$
that is corresponding to the $k$th component.

Corresponding to
$\bm{Y}=\bm{M}_{0}\bm{\mu}_{0}+\hbox{diag}(\bm{M}_{0}{\bm{\sigma}}_{0})\bm{U}$
in (22), we can write $\bm{Y}^{\prime}$ as

$\bm{Y}^{\prime}=\bm{M}_{\tau}\bm{\mu}+\hbox{diag}(\bm{M}_{\tau}{\bm{\sigma}})\bm{U}=\bm{M}_{\tau}\mathbb{A}_{\tau}(\bm{Y}^{\prime})+\hbox{diag}\left\\{\bm{M}_{\tau}\mathbb{B}_{\tau}(\bm{Y}^{\prime})\right\\}\mathbb{C}_{\tau}(\bm{Y}^{\prime}),$

where
$\mathbb{A}_{\tau}(\bm{Y}^{\prime})=(\bm{M}_{\tau}^{\top}\bm{M}_{\tau})^{-1}\bm{M}_{\tau}^{\top}\bm{Y}^{\prime}$,

$\mathbb{B}_{\tau}(\bm{Y}^{\prime})=\begin{pmatrix}\|(I-\bm{H}_{\tau}^{(1)})\bm{Y}^{{}^{\prime}(1)}\|\\\
\vdots\\\
\|(I-\bm{H}_{\tau}^{(\tau)})\bm{Y}^{{}^{\prime}(\tau)}\|\end{pmatrix}\,\,\hbox{and}\,\,\mathbb{C}_{\tau}(\bm{Y}^{\prime})=\begin{pmatrix}\frac{(I-\bm{H}_{\tau}^{(1)})\bm{Y}^{{}^{\prime}(1)}}{\|(I-\bm{H}_{\tau}^{(1)})\bm{Y}^{{}^{\prime}(1)}\|}\\\
\vdots\\\
\frac{(I-\bm{H}_{\tau}^{(\tau)})\bm{Y}^{{}^{\prime}(\tau)}}{\|(I-\bm{H}_{\tau}^{(\tau)})\bm{Y}^{{}^{\prime}(\tau)}\|}\end{pmatrix}=\begin{pmatrix}\frac{(I-\bm{H}_{\tau}^{(1)})\bm{U}^{(1)}}{\|(I-\bm{H}_{\tau}^{(1)})\bm{U}^{(1)}\|}\\\
\vdots\\\
\frac{(I-\bm{H}_{\tau}^{(\tau)})\bm{U}^{(\tau)}}{\|(I-\bm{H}_{\tau}^{(\tau)})\bm{U}^{(\tau)}\|}\end{pmatrix}.$

For the given $(\tau,\bm{M}_{\tau})$,
$\\{\mathbb{A}_{\tau}(\bm{Y}^{\prime}),\mathbb{B}_{\tau}(\bm{Y}^{\prime})\\}$
is the sufficient statistics for $(\bm{\mu},\bm{\sigma})$ and
$\mathbb{C}_{\tau}(\bm{Y}^{\prime})$ is the remaining normalized “noises”
whose distribution is free of nuisance parameters $(\bm{\mu},\bm{\sigma})$. By
its expression, $\mathbb{C}(\bm{Y}^{\prime})$ depends only on $\bm{U}$, so we
also write $\mathbb{C}(\bm{Y}^{\prime})=\mathbb{C}(\bm{U})$. To develop a
nuclear mapping for $\tau$ while removing the impact of the nuisance
parameters $(\bm{\mu},\bm{\sigma})$, we consider below a conditional
probability given
$\\{\mathbb{A}_{\tau}(\bm{Y}^{\prime}),\mathbb{B}_{\tau}(\bm{Y}^{\prime})\\}$,
the sufficient statistics for $(\bm{\mu},\bm{\sigma})$.

Let $\hat{\tau}(\bm{Y}^{\prime})$ be any reasonable estimates of $\tau$ based
on $\bm{Y}^{\prime}.$ For instance, we can use a BIC-based estimator

$\displaystyle\hat{\tau}(\bm{Y}^{\prime})=\arg\min_{\tau\in{\it\Upsilon}}\left[\left\\{-2\max_{(\bm{M}_{\tau},{\bm{\mu}},\bm{\sigma})}\log\ell(\bm{M}_{\tau},{\bm{\mu}},\bm{\sigma}|\bm{Y}^{\prime})\right\\}+2\tau\log(n)\right],$

where $\ell(\bm{M}_{\tau},{\bm{\mu}},\bm{\tau}|\bm{Y}^{\prime})$ is the log-
likelihood function of $(\bm{M}_{\tau},{\bm{\mu}},\bm{\tau})$ given
$\bm{Y}^{\prime}.$ We also denote the conditional probability
$\mathcal{P}(k|\bm{a},\bm{b})=\mathbb{P}\big{\\{}\hat{\tau}(\bm{Y}^{\prime})=k\big{|}\mathbb{A}_{\tau}(\bm{Y}^{\prime})=\bm{a},\mathbb{B}_{\tau}(\bm{Y}^{\prime})=\bm{b}\big{\\}}$
and let

$\displaystyle\mathcal{F}_{\theta}(w|\bm{a},\bm{b})=\sum_{\\{k:\mathcal{P}(k|\bm{a},\bm{b})>\mathcal{P}(w|\bm{a},\bm{b})\\}}\mathcal{P}(k|\bm{a},\bm{b}).$

Conditional on
$\big{\\{}A_{\tau}(\bm{Y}^{\prime})=\bm{a},B_{\tau}(\bm{Y}^{\prime})=\bm{b}\big{\\}}$,
$\mathcal{F}_{\theta}(w|\bm{a},\bm{b})$ does not involve often unknown
$(\bm{\mu},\bm{\sigma})$, as stated by the lemma below. See Appendix I for a
proof of the lemma.

######  Lemma 5.

The function defined above $\mathcal{F}_{\theta}(w|\bm{a},\bm{b})$ is free of
$(\bm{\mu},\bm{\sigma})$ and so we can express
$\mathcal{F}_{\theta}(w|\bm{a},\bm{b})=\mathcal{F}_{(\tau,M_{\tau})}(w|\bm{a},\bm{b}).$

Next we follow Section 3.2 to remove the impact of the remaining latent
discrete nuisance parameter $\bm{M}_{\tau}$. Specifically, let
$T_{a}(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)=(\mathbb{A}_{\tau}(\bm{y}^{\prime}),\mathbb{B}_{\tau}(\bm{y}^{\prime}),$
$\mathbb{C}_{\tau}(\bm{y}^{\prime})),$ where
$\bm{y}^{\prime}=\bm{M}_{\tau}\bm{\mu}+\hbox{diag}(\bm{M}_{\tau}{\bm{\sigma}})\bm{u}$,
and $T_{a}(\bm{u},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ it
is a mapping from ${\cal U}\times\Theta\to{\cal T}$. We define a set in $\cal
T$

$\displaystyle B_{\theta}(\alpha)$
$\displaystyle=\left\\{(\bm{a},\bm{b},\bm{c}):\mathcal{F}_{\theta}(\hat{\tau}(\bm{y}^{*})|\bm{a},\bm{b})\leq\alpha\,\,\text{for}\,\,\bm{y}^{*}=\bm{M}_{\tau}\bm{a}+\hbox{diag}\\{\bm{M}_{\tau}\bm{b}\\}\bm{c}\right\\}$
$\displaystyle=\left\\{(\bm{a},\bm{b},\bm{c}):\mathcal{F}_{(\tau,\bm{M}_{\tau})}(\hat{\tau}(\bm{y}^{*})|\bm{a},\bm{b})\leq\alpha\,\,\text{for}\,\,\bm{y}^{*}=\bm{M}_{\tau}\bm{a}+\hbox{diag}\\{\bm{M}_{\tau}\bm{b}\\}\bm{c}\right\\}$
$\displaystyle\overset{def}{=}B_{(\tau.\bm{M}_{\tau})}(\alpha).$

This set only depends on the given $(\tau,\bm{M}_{\tau})$ and is free of
nuisance parameters $(\bm{\mu},\bm{\sigma})$. Theorem 5 below states that
$T_{a}(\bm{U},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)$ is
inside $B_{(\tau.\bm{M}_{\tau})}(\alpha)$ with a probability greater than
$\alpha$, both conditionally and unconditionally. See Appendix I for a proof.

######  Theorem 5.

For any $(\bm{a},\bm{b})$, the conditional probability

$\displaystyle\mathbb{P}\\{T_{a}(\bm{U},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\in
B_{(\tau.\bm{M}_{\tau})}(\alpha)|\mathbb{A}_{\tau}(\bm{Y}^{\prime})=\bm{a},\mathbb{B}_{\tau}(\bm{Y}^{\prime})=\bm{b}\\}\geq\alpha.$
(24)

Therefore,
$\mathbb{P}\\{T_{a}(\bm{U},\mathop{\theta\kern-4.49997pt\hbox{$\theta$}}\nolimits)\in
B_{(\tau.\bm{M}_{\tau})}(\alpha)\\}\geq\alpha.$

For $\bm{y}^{\prime}=G(\bm{u},\tau,\bm{M}_{\tau},\bm{\mu},\bm{\sigma})$ and
any $(\tau,\widetilde{\bm{M}_{\tau}})$, following (13), we define

$\displaystyle\nu(\bm{u},\tau,\bm{M}_{\tau},\bm{\mu},\bm{\sigma},\widetilde{\bm{M}_{\tau}})=\tilde{\nu}(\bm{y}^{\prime},\tau,\bm{M}_{\tau},\widetilde{\bm{M}_{\tau}})=\inf_{\alpha^{\prime}}\\{\alpha^{\prime}:(\tau,\widetilde{\bm{M}}_{\tau})\in\Gamma_{\alpha^{\prime}}(\bm{y}^{\prime})\\},$

where
$\Gamma_{\alpha}(\bm{y}^{\prime})=\\{(\tau,\bm{M}_{\tau}):T_{a}(\bm{u},\bm{\theta})\in
B_{(\tau.\bm{M}_{\tau})}(\alpha)\\}=\\{(\tau,\bm{M}_{\tau}):(\mathbb{A}_{\tau}(\bm{y}^{\prime}),$
$\mathbb{B}_{\tau}(\bm{y}^{\prime}),\mathbb{C}_{\tau}(\bm{y}^{\prime}))\in
B_{(\tau.\bm{M}_{\tau})}(\alpha)\\}$. Then the (profile) nuclear mapping
function for $\tau$ is

$\displaystyle
T^{P}(\bm{u},\tau,\bm{M}_{\tau},\bm{\mu},\bm{\sigma})=\widetilde{T}^{P}(\bm{y}^{\prime},\tau)=\min_{\widetilde{\bm{M}}_{\tau}}\tilde{\nu}(\bm{y}^{\prime},\tau,\bm{M}_{\tau},\widetilde{\bm{M}_{\tau}}).$

According to Lemma 2 and (3.2), we define

$\displaystyle\Xi_{\alpha}(\bm{y}_{obs})$
$\displaystyle=\big{\\{}\tau:\exists\bm{u}^{*}\in{\cal U}\mbox{ and
$(\bm{M}_{\tau},\bm{\mu},\bm{\sigma})$ s.t.
}\bm{y}_{obs}=\bm{\mu}+\hbox{diag}\\{\bm{M}_{\tau}\bm{\sigma}\\}\bm{u}^{*},$
$\displaystyle\qquad\qquad
T^{P}(\bm{u}^{*},\tau,\bm{M}_{\tau},\bm{\mu},\bm{\sigma})\leq\alpha\big{\\}}$
(25)

From Theorem 3 and a direct calculation, the corollary below states that
$\Xi_{\alpha}(\bm{y}_{obs})$ can be further simplified and is a level $\alpha$
confidence set for $\tau_{0}$. See Appendix I for a proof.

######  Corollary 4.

The set defined in (25) can be further expressed as

$\displaystyle\Xi_{\alpha}(\bm{y}_{obs})=\\{\tau:\exists\widetilde{\bm{M}}_{\tau},\mathcal{F}_{\tau,\widetilde{\bm{M}}_{\tau}}(\hat{\tau}(\bm{y}^{obs})|\widetilde{\mathbb{A}}_{\tau}(\bm{y}_{obs}),\widetilde{\mathbb{B}}_{\tilde{\tau}}(\bm{y}_{obs}))\leq\alpha\\},$
(26)

where $\widetilde{\mathbb{A}}_{\tau}(\bm{y}_{obs})$ and
$\widetilde{\mathbb{B}}_{\tau}(\bm{y}_{obs})$ are defined by replacing
$\bm{M}_{\tau}$ with $\widetilde{\bm{M}}_{\tau}$ in the definition of
$\mathbb{A}_{\tau}(\bm{y}_{obs})$ and $\mathbb{B}_{\tau}(\bm{y}_{obs}).$ Also,
$\mathbb{P}\big{\\{}\tau_{0}\in\Xi_{\alpha}(\bm{Y})\big{\\}}\geq\alpha$, so
$\Xi_{\alpha}(\bm{y}_{obs})$ is a level $\alpha$ confidence set for
$\tau_{0}$.

Corollary 4 suggests that, to construct $\Xi_{\alpha}(\bm{y}_{obs}),$ we only
need to calculate the function
$\mathcal{F}_{\tau,\widetilde{\bm{M}}_{\tau}}(\hat{\tau}(\bm{y}^{obs})|\bm{a},\bm{b})$
in (26) for each pair of $(\tau,\widetilde{\bm{M}}_{\tau})$ and
$\\{\bm{a},\bm{b}\\}=\\{\widetilde{\mathbb{A}}_{\tau}(\bm{y}_{obs}),\widetilde{\mathbb{B}}_{\tilde{\tau}}(\bm{y}_{obs})\\}$.
This can be achieved with Monte-Carlo approach as follows. Given
$(\tau,\widetilde{\bm{M}}_{\tau})$ and $(\bm{a},\bm{b}),$ we simulate a large
number of $\bm{u}^{s}\sim\bm{U}$ and the collection of $\bm{u}^{s}$’s forms a
set $\mathcal{V}$. Let

$\displaystyle\bm{y}^{s}_{\widetilde{\bm{M}}_{\tau},\bm{a},\bm{b}}=\widetilde{\bm{M}}_{\tau}\bm{a}+\hbox{diag}\\{\widetilde{\bm{M}}_{\tau}\bm{b}\\}\mathbb{C}_{\tau}(\bm{u}^{s}).$

The following theorem validates that we can approximate the conditional
distribution of $\bm{Y}$ given $\mathbb{A}_{\tau}(\bm{Y})=\bm{a}$ and
$\mathbb{B}_{\tau}(\bm{Y})=\bm{b}$ with the Monte-Carlo sampling distribution
of $\bm{y}^{s}_{\widetilde{\bm{M}}_{\tau},\bm{a},\bm{b}}.$ See Appendix I for
a proof.

######  Lemma 6.

Given a $(\tau,\widetilde{\bm{M}}_{\tau}),$ let
$\widetilde{\bm{Y}}=\widetilde{\bm{M}}_{\tau}{\bm{\mu}}_{\tau}+\hbox{diag}(\bm{U})\widetilde{\bm{M}}_{\tau}{\bm{\sigma}_{\tau}}$
and
$\bm{Y}^{s}_{\widetilde{\bm{M}}_{\tau},\bm{a},\bm{b}}=\widetilde{\bm{M}}_{\tau}\bm{a}+\hbox{diag}\\{\widetilde{\bm{M}}_{\tau}\bm{b}\\}\bm{C}_{\tau}(\bm{U}^{s}),$
where $\bm{U}^{s}\sim\bm{U},$ then

$\displaystyle\left\\{\widetilde{\bm{Y}}\middle|\widetilde{\mathbb{A}}_{\tau}(\widetilde{\bm{Y}})=\bm{a},\widetilde{\mathbb{B}}_{\tau}(\widetilde{\bm{Y}})=\bm{b}\right\\}\quad\sim\quad\bm{Y}^{s}_{\widetilde{\bm{M}}_{\tau},\bm{a},\bm{b}}.$

When $|\cal V|$ is large enough, the function
$\mathcal{F}_{\tau,\widetilde{\bm{M}}_{\tau}}(w|\bm{a},\bm{b})$ in (26) can be
approximated by the Monte-Carlo version

$\displaystyle\mathcal{F}_{\tau,\widetilde{\bm{M}}_{\tau}}^{s}(w|\bm{a},\bm{b})\overset{def}{=}\sum_{\\{k:{\mathcal{P}}_{\mathcal{V}}^{s}(k|\bm{a},\bm{b})>{\mathcal{P}}_{\mathcal{V}}^{s}(w|\bm{a},\bm{b})\\}}{\mathcal{P}}_{\mathcal{V}}^{s}(k|\bm{a},\bm{b}),$

since based on Lemma 6,
$\mathcal{F}_{\tau,\widetilde{\bm{M}}_{\tau}}^{s}(w|\bm{a},\bm{b})\to\mathcal{F}_{\tau,\widetilde{\bm{M}}_{\tau}}(w|\bm{a},\bm{b}),$
as $|\cal V|\to\infty$. Here,
$\mathcal{P}_{\mathcal{V}}^{s}(k|\bm{a},\bm{b})=\frac{1}{|\mathcal{V}|}\sum_{\bm{u}^{s}\in\mathcal{V}}I(\hat{\tau}(\bm{y}^{s}_{\widetilde{\bm{M}}_{\tau},\bm{a},\bm{b}})=k)$
and
$\bm{y}^{s}_{\widetilde{\bm{M}}_{\tau},\bm{a},\bm{b}}=\widetilde{\bm{M}}_{\tau}\bm{a}+\hbox{diag}\\{\widetilde{\bm{M}}_{\tau}\bm{b}\\}\mathbb{C}_{\tau}(\bm{u}^{s})$,
for $\bm{u}^{s}\in{\cal V}$.

###  4.3 Finding candidates for $(\tau_{0},\bm{M}_{0})$

We follow Section 3.3 to construct a data-driven candidate set
$\widehat{\it\Upsilon}(\bm{y}_{obs})$ for $\kappa_{0}=(\tau_{0},\bm{M}_{0})$
based on the observed data $\bm{y}_{obs}.$ Corresponding to (19), we use a
modified BIC as our lost function and define our many-to-one mapping for
$\kappa^{*}=(\tau^{*},\bm{M}^{*}_{\tau})$ as

$(\tau^{*},\bm{M}^{*}_{\tau^{*}})=\arg\min\limits_{(\tau,\bm{M}_{\tau})}\min\limits_{(\bm{\mu},\bm{\sigma})}\left\\{n\log\bigg{(}\frac{\|\bm{y}_{obs}-{\bm{M}}_{\tau}{\bm{\mu}}_{\tau}-\hbox{diag}(\bm{M}_{\tau}{\bm{\sigma}_{\tau}})\bm{u}^{*}\|^{2}+1}{n}\bigg{)}+2\lambda\tau\log(n)\right\\}.$

(27)

By (20), a data-dependent candidate set for $\kappa_{0}=(\tau_{0},\bm{M}_{0})$
is

where ${\cal V}_{c}=\\{\bm{u}_{1},\ldots,\bm{u}_{N}\\}$ is a collection of
simulated $\bm{u}_{j}^{c}\sim N(\bm{0},\bm{I})$, $j=1,\ldots,N$.

Then based on (3.3) and (26), we define the confidence set

(28)

One complication under the mixture model setting is that the equation (21) in
(N1) may not hold when $\bm{u}^{rel}$ hits certain directions (i.e., the two
vectors $\hbox{diag}(\bm{M}_{0}\sigma_{0})\bm{u}^{rel}$ and
$\\{I-\bm{M}_{\tau^{\prime}}(\bm{M}_{\tau^{\prime}}^{\top}\bm{M}_{\tau^{\prime}})^{-1}\bm{M}_{\tau^{\prime}}^{\top}\\}\bm{M}_{0}\bm{\mu}_{0}$
are on the same direction for some $\bm{M}_{\tau^{\prime}}\not=\bm{M}_{0}$,
$\tau^{\prime}\leq\tau_{0}$) or $\|\bm{u}^{rel}\|$ is extremely large.
Fortunately, the probabilities of these events are arbitrarily small,
resulting in the arbitrarily small difference $\delta$ between the probability
bound in the following Theorem 6 and that in Theorem 4. In the theorem,
$C_{\min}=\underset{\small\\{\bm{M}_{\tau}:\,|\tau|\leq|\tau_{0}|,\tau\neq\tau_{0}\\}}{\min}\|(\bm{I}-\bm{M}_{\tau})\bm{M}_{0}\bm{\mu}_{0}\|>0.$
See Appendix I for a proof.

######  Theorem 6.

Suppose $n-\tau_{0}>4.$ Let the distance metric between $\bm{U}$ and
$\bm{U}^{*}$ be
$\rho(\bm{U},\bm{U}^{*})=\frac{\|(\bm{I}-\bm{H}_{\bm{U}^{*}})\bm{U}\|}{\|\bm{U}^{*}\|},$
then for any arbitrarily small $\delta>0,$ there exist $\gamma_{\delta}>0,c>0$
and also an interval of positive width
$\Lambda_{\delta}=[\frac{\gamma_{\delta}^{1.5}}{\log(n)/n},$
$\frac{\log\\{0.5C^{2}_{\min}\sigma_{\max}\gamma_{\delta}/\sigma_{\min}+1\\}}{2\tau_{0}\log(n)/n}]$
such that when $\lambda\in\Lambda_{\delta},$ (a) $\mathbb{P}_{\bm{U}|{\cal
V}_{c}}\\{(\tau_{0},\bm{M}_{0})\in\widehat{\it\Upsilon}_{{\cal
V}_{c}}(\bm{Y})\\}$ $\geq 1-\delta-o_{p}(e^{-c|\mathcal{V}^{c}|})$, (b)
$\mathbb{P}_{\bm{U}|{\cal
V}_{c}}\\{(\tau_{0},\bm{M}_{0})\in\Xi^{\prime}_{\alpha}(\bm{Y})\\}\geq\alpha-\delta-
o_{p}(e^{-c|\mathcal{V}^{c}|}),$ where $\Xi^{\prime}_{\alpha}(\cdot)$ is
defined in (28).

Remark 4.1 Compared to the standard BIC, the modified BIC many-to-one mapping
(27) has an extreme term that involves the repro copy $\bm{u}^{*}$. We
modified the term inside the log function so that the mapping (27) can give us
$(\tau_{0},M_{0})$ when $\bm{u}^{*}$ equals $\bm{u}^{rel}$. We also allow the
penalty on the model size to change by using $\lambda\log(n)(2\tau)$ instead
of $\log(n)(2\tau)$. This change adds some flexibility to help control the
finite sample coverage probability as shown in Theorem 6. In practice, we
usually still use $\lambda=1,$ since it is often in its required range for a
small $\delta$. For example, to achieve a coverage error $\delta=O(e^{-n}),$
we can take $\gamma_{\delta}=O(1/n)$. Considering $C_{\min}$ is usually
$O(n),$ it is likely that
$\frac{\gamma_{\delta}^{1.5}}{\log(n)/n}<1<\frac{\log\\{0.5C^{2}_{\min}\gamma_{\delta}\sigma_{\max}/\sigma_{\min}+1\\}}{2\tau_{0}\log(n)/n}$.

##  5 Numerical Studies: Real and simulated data

###  5.1 Real data analysis: red blood cell sodium-lithium countertransport
(SLC) data

In this section we consider a data set that consists of Red blood cell sodium-
lithium countertransport (SLC) activity measurements from 190 individuals
(Dudley et al., 1991). SLC measurement is known to be correlated with blood
pressure, and perceived as an essential cause of hypertension by some
researchers. Moreover, SLC is usually easier to study than blood pressure that
could be influenced by numerous environmental and genetic factors (Dudley et
al., 1991). Roeder (1994) and Chen et al. (2012) both analyzed the SLC data
set of 190 individuals using a normal mixture model, and perform a hypothesis
testing on the number of components $H_{0}:\tau_{0}=k$ vs $H_{1}:\tau_{0}>k$
for a small integer $k.$ Roeder (1994) concluded that $\tau_{0}=3$ is the
smallest value of $\tau_{0}$ not rejected by the data assuming variance of
each component is the same. Chen et al. (2012) argued instead that dropping
the assumption of equal variances, the hypothesis of $\tau_{0}=2$ would not be
rejected and thus a good fit for the data.

We use the repro samples approach developed in Section 4 to re-analyze the SLC
data. Following Roeder (1994); Chen et al. (2012), we assume that the SLC
measurements follow a normal mixture distribution with unknown number of
components $\tau_{0}$ and obtain a 95% confidence set for $\tau_{0}$ that is
$\\{2,3,4\\}.$ After that, we fit the normal mixture model with $\tau=2,3,4$
on the SLC data with an EM algorithm, respectively. We provide the estimated
mean, standard deviation and weight of each component in Table 4. The
confidence set obtained by our repro samples approach suggests that both
hypothesis $\tau_{0}=2$ and $\tau_{0}=3$ are plausible. Since, in genetics, a
two-component mixture distribution corresponds to a simple dominance model and
a three-component corresponds to an additive model (Roeder, 1994), evidently
we can not rule out either of the two competing genetic models based on the
SLC data alone.

An advantage of the proposed repor samples method is that its confidence set
provides both a lower bound and an upper bound for plausible values of
$\tau_{0}$ supported by the data set, whereas inverting one-sided tests
proposed by Roeder (1994) and Chen et al. (2012) can only lead to one sided
confidence set. One formidable challenge of finding an upper bound is the fact
that any distribution can be arbitrarily well approximated by a normal mixture
distribution, therefore $H_{0}:\tau_{0}=k$ for a large $k$ can never be
rejected. Here we manage to tackle this challenge by using repro samples to
find a proper candidate set as described in Section 4.3. The upper bound of
the confidence set for $\tau_{0}$ we find for the SLC data is $4.$ In Figure
3, we plot the histogram of the SLC measurements together with estimated
normal mixture densities with number of components $\tau=2,3,4$ using EM
algorithm. It appears that only with $\tau=4$, the normal mixture model could
capture the three spikes in the histogram, while still smoothly fitting the
rest of the histogram. Indeed, a normal mixture model with four components
should demand further investigation, since it effectively represents the data
without over fitting it.

$\tau_{0}$ | Component Mean | Component SD | Weight
---|---|---|---
2 | (0.2206,0.3654) | (0.0571,0.1012) | (0.7057,0.2943)
3 | (0.1887,0.4199,0.2809) | (0.0414,0.0886,0.0474) | (0.4453,0.168,0.3866)
4 | (0.1804,0.3351,0.2556,0.4403) | (0.0362,0.0359,0.0268,0.086) | (0.4018,0.1742,0.2941,0.1299)

Table 4: Estimated mean, standard deviation and weight of each component of
the normal mixture distributions for the SLC data. Figure 3: Histogram of 190
SLC measurements, including three estimated normal mixture density curves with
number of components $\tau_{0}=2,3,4.$

We also conducted a Bayesian analysis on the SLD data using R package mixAK
(Komárek and Komárková, 2014), where the Bayesian inference procedure on
Gaussian mixture model with unknown number of components is discussed in
Richardson and Green (1997). We implement four different priors on $\tau$ in
our analysis: uniform, truncated Poisson(0.4), truncated Poisson(1) and
truncated Poisson(5) over $\tau=1,\dots,10.$ The priors used on $\bm{\mu}$ and
$\bm{\sigma}$ are the default priors (i.e., Gaussian and inverse Gamma priors)
in mixAK. The posterior distributions for the four priors are plotted in
Figure 4, with corresponding 95% credible sets $\\{2,3\\}$, $\\{2\\}$,
$\\{2,3\\}$, $\\{2,3,4\\}$, respectively. It turns out the posterior
distributions and their corresponding 95% credible intervals (sets) of $\tau$
are quite different for different priors. Even in the cases with the same
credible set using uniform and truncated Poisson(1) priors, their posterior
distributions are quite different. It is clear that the Bayesian inference of
$\tau$ is very sensitive to the choice of prior distribution. A simulation
study in the next section affirms this observation. In addition, the
simulation study further suggests that the Bayesian procedure is not suited
for the task of recovering the true value $\tau_{0}$ (under repeated
experiments) and Bayesian outcomes are also affected greatly by the default
priors on $\bm{\mu}$ and $\bm{\sigma}$.

Figure 4: Posterior distributions of $\tau$ under four different priors: (a)
uniform; (b) Poisson($\mu=0.4$); (c) Poisson($\mu=1$); (d) Poisson($\mu=5$);
All priors are truncated outside the set $\\{\tau:\tau=1,\dots,10\\}$.

###  5.2 Simulation Study: Inference and Recovery of True $\tau_{0}$

To further demonstrate the performance of the proposed confidence set, we
carry out simulation studies in which the simulated data are generated using
the estimated normal mixture model with $\tau_{0}=3$ and $\tau_{0}=4$ in Table
4, respectively. We keep the sample size at $n=190$. For each simulated data,
we apply the proposed repro samples method. We repeat the simulation and
analysis 200 times, and subsequently summarize the 200 level-95% confidence
sets using a bar plot in Figure 5, for $\tau_{0}=3$ and $\tau_{0}=4$,
respectively.

Figure 5: Bar plots of the 200 level 95% repro samples confidence sets
obtained in the simulation study (200 repetitions). The upper panel is for the
truth $\tau_{0}=3$; lower for the truth $\tau_{0}=4.$

It appears that for both settings of $\tau_{0}=3$ and $\tau_{0}=4$, the
proposed confidence set for $\tau_{0}$ achieves coverage rates of 99% and
96.5% respectively, both greater than 95%. Note that the parameter space for
$\tau_{0}$ is discrete, therefore, unlike confidence intervals for continuous
parameters, it is not always possible to achieve an exact confidence level of
95%. For simulated data with $\tau_{0}=3,$ the average size of the proposed
confidence set is about 3, with a standard error of 0.05. In fact, we observe
from Figure 5 that for the majority of the 200 simulations, our approach
produces a confidence set of $\\{2,3,4\\}.$ As for $\tau_{0}=4,$ more than 80%
of the times, we obtain a confidence set of either $\\{2,3,4\\}$ or
$\\{2,3,4,5\\}$, making the average size of the confidence set around 3.65,
with a standard error of 0.05.

We compare the performance of the existing methods under the same simulation
settings. First, the classical BIC method recovers the true value of
$\tau_{0}=3$ only 13 times (6.5%) out of the 200 repetitions. The number of
times drops to $0$ (out of 200 repetitions) when the true $\tau_{0}=4.$ The
results are not surprising, since with penalty terms point estimators of
$\tau_{0}$ are often biased towards smaller $\tau$ values in practice, even
though BIC estimator is a consistent estimator. These two extremely low rates