Adversarial Meta-Learning of Gamma-Minimax Estimators That Leverage Prior Knowledge

Bayes estimators are well known to provide a means to incorporate prior knowledge that can be expressed in terms of a single prior distribution. However, when this knowledge is too vague to express with a single prior, an alternative approach is needed. Gamma-minimax estimators provide such an approach. These estimators minimize the worst-case Bayes risk over a set $\Gamma$ of prior distributions that are compatible with the available knowledge. Traditionally, Gamma-minimaxity is defined for parametric models. In this work, we define Gamma-minimax estimators for general models and propose adversarial meta-learning algorithms to compute them when the set of prior distributions is constrained by generalized moments. Accompanying convergence guarantees are also provided. We also introduce a neural network class that provides a rich, but finite-dimensional, class of estimators from which a Gamma-minimax estimator can be selected. We illustrate our method in two settings, namely entropy estimation and a prediction problem that arises in biodiversity studies.


Introduction
A variety of principles can be used to guide the search for a suitable statistical estimator.Asymptotic efficiency (Pfanzagl, 1990), minimaxity (Wald, 1945) and Bayes optimality (Berger, 1985) are popular examples of such principles.Defining the performance criteria underlying these principles requires specifying a model space, that is, a collection of possible data-generating mechanisms known to contain the true, underlying distribution.
It is often desirable to incorporate prior information about the true data-generating mechanism into a statistical procedure.This might be done differently in different statistical paradigms.For frequentist methods, such as those based on the asymptotic efficiency or minimax principle, the primary way to incorporate this information is via the definition of the model space.In the Bayesian paradigm, such information may be represented by further specifying a prior distribution (or prior for short) over the model space and aiming for an estimator that minimizes the induced Bayes risk.However, in many cases, there may be several priors that are compatible with the available information, especially in the context of rich model spaces.
The Gamma-minimax paradigm, proposed by Robbins (1951), provides a principled means to overcome the challenge of specifying a single prior distribution.Under this paradigm, a statistician first specifies a set Γ of all priors that are consistent with the available prior information and subsequently seeks an estimator that minimizes the worst-case Bayes risk over this set of priors.The Gamma-minimax paradigm may be viewed as a robust version of the Bayesian paradigm that is less sensitive to misspecification of a prior distribution (Vidakovic, 2000).When it is infeasible to specify a prior due to the complexity of the model space, the Gamma-minimax paradigm may also be viewed as a feasible substitute for the Bayesian paradigm.The Gamma-minimax paradigm is closely related to Bayes and minimax paradigms: when the set of priors consists of a single prior, a Gamma-minimax estimator is Bayes with respect to that prior; when the set Γ of priors is the entire set of possible prior distributions, a Gamma-minimax estimator is also minimax.
Gamma-minimax estimators have been studied for a variety of problems.Some explicit forms of Gamma-minimax estimators have been obtained.For example, Olman and Shmundak (1985) studied Gamma-minimax estimation of the mean of a normal distribution for the set of symmetric and unimodal priors on an interval and obtained an explicit form when this interval is sufficiently small.Eichenauer-Herrmann (1990) generalized this result to more general parametric models and Eichenauer-Herrmann et al. (1994) obtained a further generalization with the requirement of symmetry on the priors dropped.Chen et al. (1988) studied Gamma-minimax estimation for multinomial distributions and the set of priors with bounded mean.Chen et al. (1991) studied Gamma-minimax estimation for one-parameter exponential families and the set of priors that place certain bounds on the first two moments.These results do not deal with general model spaces, such as semiparametric or nonparametric models, and general forms of the set of priors that may not directly impose bounds on prior moments of the parameters of interest.One reason for this lack of generality might be that, in the existing literature, Gammaminimaxity is defined only for parametric models.However, an issue with parametric models is that they often fail to contain the true data-generating mechanism, in which case output from the aforementioned statistical procedures may no longer be interpretable.Another possible reason is that it is typically intractable to analytically derive Gamma-minimax estimators, even for parametric models.
To overcome this lack of analytical tractability, meta-learning algorithms to compute a minimax or Gamma-minimax estimator have been proposed.Still, most of these works focus on parametric models.
For example, Nelson (1966) and Kempthorne (1987) each proposed an algorithm to compute a minimax estimator.Bryan et al. (2007) and Schafer and Stark (2009) proposed an algorithm to compute an approximate confidence region of optimal expected size in the minimax sense.Noubiap and Seidel (2001) proposed an iterative algorithm to compute a Gamma-minimax decision for the set of priors constrained by generalized moment conditions.More recent works explored computing estimators under more general models.For example, Luedtke et al. (2020a) introduced an approach, termed Adversarial Monte Carlo meta-learning (AMC), for constructing minimax estimators.In the special case of prediction problems with mean-squared error, Luedtke et al. (2020b) studied the invariance properties of the decision problem and their implications for AMC.
In this paper, we make the following contributions: 1. We propose iterative adversarial meta-learning algorithms for constructing Gamma-minimax estimators for a general model space and class of estimators.We further provide convergence guarantees for these algorithms.
To our best knowledge, this is the first algorithm to compute Gamma-minimax estimators under general models, including infinite-dimensional models.We also show that, for certain problems, there is a unique Gamma-minimax estimator and our computed estimator converges to this estimator as the number of iterations increases to infinity.
Like the approach proposed in Noubiap and Seidel (2001), we consider sets of priors characterized by (in)equality constraints on prior generalized moments and our proposed iterative algorithm involves solving a discretized Gamma-minimax optimization problem in each intermediate step.However, we explicitly describe algorithms to solve these minimax problems, which facilitates the use of our approach by practitioners.When the space of estimators can be parameterized by a Euclidean space and gradients are available, we propose to use a gradient-based algorithm or a stochastic variant thereof.When gradients are unavailable, we propose to instead use fictitious play (Brown, 1951;Robinson, 1951) to compute a stochastic estimator, which is a mixture of deterministic estimators belonging to some specified collection.We also provide a convergence result that is applicable even when this collection has infinite cardinality.This is in contrast to the results in Robinson (1951), which require that each player has only finitely many possible deterministic strategies.
2. We propose a Markov chain Monte Carlo (MCMC) method to construct the approximating grids defining the discretized Gamma-minimax problems used in our iterative scheme.
Like the approach proposed in Noubiap and Seidel (2001), our proposed iterative algorithm relies on increasingly fine finite grids over the model space.However, since we allow the model space to be high or even infinite-dimensional, randomly adding points to the grid may lead to unacceptably slow convergence.
To overcome this challenge, we propose to use MCMC to efficiently construct such grids.
Our theoretical results allow for many different choices of classes of estimators.Our final contribution concerns the introduction of one such class: 3. We introduce a new neural network architecture that can be used to parameterize statistical estimators and argue that this class represents an appealing class to optimize over.
For this final point, we build on existing works in adversarial learning (e.g., Goodfellow et al., 2014;Luedtke et al., 2020a,b) and extreme learning machines (Huang et al., 2006b).Thanks to the universal approximation properties of neural networks (e.g., Hornik, 1991;Csáji, 2001) and extreme learning machines (Huang et al., 2006a), we also show that both of these parameterizations can achieve good performance for sufficiently large networks.Furthermore, inspired by pre-training (e.g., Erhan et al., 2010) and transfer learning (e.g., Torrey and Shavlik, 2009), we recommend leveraging knowledge of existing estimators as inputs to the network in settings where this is possible.Under such choices of the space of estimators, we can expect to obtain a useful estimator even if the associated nonconvex-concave minimax problems prove to be difficult.This paper is organized as follows.In Section 2, we introduce the framework of Gamma-minimax estimation and regularity conditions that we assume throughout the paper.In Section 3, we describe our proposed iterative adversarial meta-learning algorithms.In Section 4, we discuss considerations when choosing hyperparameters in the algorithms.In Section 5, we demonstrate our method in three simulation studies.We conclude with a discussion in Section 6. Proof sketches of key results are provided in the main text, and complete proofs can be found in the appendix.We also provide a table summarizing the frequently used symbols in Table 7 in the appendix.The code for our simulations is available on GitHub (Qiu, 2022).

Problem setup
Let M be a separable Hausdorff space of data-generating mechanisms that contains the truth P 0 and is equipped with a metric ρ.Under a data-generating mechanism P ∈ M, let X * ∈ X * denote the random data being generated, where X * is the space of values that the random data takes.Let C denote a known coarsening mechanism such that the observed data X = C(X * ) ∈ X , where X is the space of observed data.In some cases, the coarsening mechanism will be the identity map, whereas in other settings, such as those in which missing, censored or truncated data is present, the coarsening mechanism will be nontrivial (e.g., Birmingham et al., 2003;Gill et al., 1997;Heitjan and Rubin, 1991;Heitjan, 1993Heitjan, , 1994)).Let D denote the space of estimators (or decision functions) equipped with a metric ϱ.In practice, for computational feasibility, we will mainly consider an estimator space D that contains estimators parameterized by a Euclidean space such as linear estimators or neural networks, and approximates a more general space D 0 , for example, the space of all estimators satisfying certain smoothness conditions.We discuss considerations concerning the choice of D in Section 4.2 and note that our proposed methods may be applied to broader estimator classes.We treat D as fixed throughout this paper.Let R : D × M → R denote a risk function that measures the performance of an estimator under a data-generating mechanism such that smaller risks are preferable.We suppose throughout that M and D are equipped with the topologies induced by ρ and ϱ, respectively.
We now present three examples in which we formulate statistical decision problems in the above form.The first example is a general example of point estimation.We use this example to illustrate how the Gamma-minimax estimation framework naturally fits into many statistical problems.The other two examples are more concrete and we will study them in the simulations and data analyses.
Example 1 (Point estimation).Suppose that M is a statistical model, which may be parametric, semiparametric, or nonparametric (Bickel et al., 1993).The data X * consists of n independently and identically distributed (iid) random variables O i , i = 1, . . ., n, following the true distribution P 0 ∈ M.
We set C to be the identity function so that X = X * .We wish to estimate an aspect Ψ(P 0 ) ∈ R of P 0 .Then, we can consider D to be a set of X → R functions and the mean-squared error risk ii) Cumulative distribution function at a point o: Example 2 (Predicting the expected number of novel categories to be observed in a new sample).Suppose that M consists of multinomial distributions with an unknown number of categories.Let an iid random sample of size n be generated from the true multinomial distribution, so that X * is a multiset containing the number X k of observations in each category k.Suppose that only categories with nonzero occurrences are observed, so that X is a left-truncated version of X * .In other words, X is the multiset C(X * ) = {X k : X k > 0}.Then, we may wish to predict the number of new categories that would be observed if a new sample of size m were collected.This problem has been extensively studied in the literature, with applications in microbiome data, species taxonomic surveys, and assessment of vocabulary size, among other areas (e.g., Shen et al., 2003;Bunge et al., 2014;Orlitsky et al., 2016).This prediction problem can be formulated in our framework.For each P ∈ M, let p k (k = 1, . . ., K P ) be the probability of category k, and Ψ(P )(X * ) be ), the expected number of new observed categories given the current full data X * .We consider D to be a set of X → R functions and set the risk to be the . This prediction problem is known to be intrinsically difficult when the future sample size m is greater than the observed sample size n, and we might expect prior information to substantially improve prediction.
Example 3 (Entropy estimation).Consider the same data-generating mechanism and observed data as in Example 2. We may wish to estimate Shannon entropy (Shannon, 1948) Ψ(P ) = − K P k=1 p k log p k , a measure of diversity.We consider D to be a set of X → R functions and set the risk to be the mean-squared error, that is, R(d, P ) = E P [{d(X) − Ψ(P )} 2 ].Jiao et al. (2015) proposed a rate-minimax estimator.Thus, in contrast to Example 2, this is an example of a statistical problem with a satisfactory solution.For these problems, we might not expect prior information to substantially improve estimation.
We now define Gamma-minimaxity within our decision-theoretic framework.We assume that M is equipped with the Borel σ-field B and let Π denote the set of all probability distributions on the measurable space (M, B).We also assume that, for any d ∈ D and any π ∈ Π, P → R(d, P ) is πintegrable.The Bayes risk corresponding to an estimator d and a prior π is defined as r : (d, π) → R(d, P ) π(dP ).Let Γ ⊆ Π be the set of priors such that all π ∈ Γ are consistent with the available prior information.An estimator is called a Γ-minimax estimator if it is in the set argmin d∈D sup π∈Γ r(d, π). (1) Throughout the rest of this paper, we assume the existence of this solution set and other solution sets to minimax problems, and that sup π∈Γ r(d, π) is finite for any d ∈ D.
In this paper, we consider the case in which Γ is characterized by finitely many generalized moment conditions, that is, where each Φ k : M → R is a prespecified function that extracts an aspect of a data-generating mechanism and c k ∈ R is a prespecified constant.The validity of our proposed template to find a Γ-minimax estimator in Section 3.1 does not require Γ to take this form, but our proposed algorithms in Sections 3.2 and 3.3 are computationally feasible for such constraints because these linear constraints lead to linear programs, which can be solved efficiently (e.g., Jiang et al., 2020).In principle, more general constraints can be handled by using suitable minimax problem solvers.Such constraints were considered in Noubiap and Seidel (2001) and can represent a variety of forms of prior information.For example, with Φ k = ±Ψ κ for some κ ≥ 1, Γ imposes bounds on prior moments of Ψ(P ); with Φ k (P ) = ±1(Ψ(P ) ∈ I) for some known interval I, Γ imposes bounds on the prior probability of Ψ(P ) lying in I. Similar prior information on aspects of P 0 other than Ψ(P 0 ) can also be represented.In addition, since an equality can be equivalently expressed by two inequalities, Γ may also impose equality constraints on prior generalized moments.Such information is commonly used to choose prior distributions in Bayesian settings (Sarma and Kay, 2020).Since we do not require specifying a parametric model or specifying an entire prior distribution for any finite-dimensional summary of P 0 , specifying a set Γ of prior distributions in the above form is no more difficult -and often easier -than specifying a single prior distribution, as would be required in a Bayesian approach.
3 Proposed meta-learning algorithms to compute a Γ-minimax estimator Since both the model space M and the estimator space D may be infinite, it is computationally infeasible to directly solve the minimax problem (1) defining a Γ-minimax estimator.Similarly to Noubiap and Seidel (2001), our general strategy is to discretize M and thus consider prior distributions with discrete supports.Once the supports of prior distributions are discrete, the optimization over prior distributions only involves finitely many parameters, namely the probability masses at support points, and thus is computationally possible.We will show that, when the grid is sufficiently fine, a solution to the discretized minimax problem is close to a solution to the original minimax problem.
Our proposed algorithm consists of two main steps.The first step is to discretize the model space M and consider an approximating grid M ℓ instead of the original complicated model space M.This discretization is illustrated in Fig. 1.We will describe M ℓ in more detail in Section 3.1.In the second step, we consider a set Γ ℓ of priors with support contained M ℓ and compute a Γ ℓ -minimax estimator.
We will describe two classes of algorithms to solve this discretized minimax problem in Sections 3.2 and 3.3, respectively. .An example of a prior distribution is displayed as black bars with their heights being proportional to the probability masses.

Grid-based approximation of Γ-minimax estimators
We first define the discretization of the model space M that we will use.Let {M ℓ } ∞ ℓ=1 be an increasing sequence of finite subsets of M such that ℓ=1 is an increasingly fine grid over M. Since M is separable, such an {M ℓ } ∞ ℓ=1 necessarily exists.Define Algorithm 1 describes how the grids M ℓ are used to compute an approximately Γ-minimax estimator in our proposed algorithms.We will show that the approximation error decays to zero as ℓ grows to infinity.Here and in the rest of the algorithms in the paper, for any real-valued function f , when we assign argmin x f (x) or argmax x f (x) to a variable, we arbitrarily pick a minimizer or maximizer if there are multiple optimizers.In practice, the user may stop the iteration at some ℓ and use a Γ ℓ -minimax estimator d * ℓ as the output estimator.We discuss the stopping criterion in more detail at the end of this section.
Algorithm 1 Iteratively approximate a Γ-minimax estimator over an increasingly fine grid.
We note that the minimax problem in Line 3 of Algorithm 1 is nontrivial to solve, and therefore we propose two algorithms that can solve this minimax problem in Sections 3.2 and 3.3.
We assume that the following conditions hold throughout the rest of the paper.
Condition 1.There exists a limit point d * ∈ D of the sequence {d * ℓ } ∞ ℓ=1 .Condition 1 holds if the sequence {d * ℓ } ∞ ℓ=1 eventually falls in a compact set.For example, if D is a space of neural networks and we take ϱ to be the Euclidean norm in the coefficient space, then we expect Condition 1 to hold if all coefficients are restricted to fall in a bounded set, which is a common restriction in theoretical analyses of neural networks (see, e.g., Goel et al., 2016;Zhang et al., 2016;Eckle and Schmidt-Hieber, 2019) and often leads to desirable generalization bounds (see, e.g., Bartlett, 1997;Bartlett et al., 2017;Neyshabur et al., 2017).Our theoretical results hold for any limit point d * in Condition 1, even if there is more than one of them.
Condition 2. The mapping d → R(d, P ) is continuous at d * for all P ∈ M.
Condition 2 also often holds.For example, when parameterized using neural networks, all estimators are continuous functions of coefficients for common activation functions such as the sigmoid or the rectified linear unit (ReLU) (Glorot et al., 2011) function, and therefore d → R(d, P ) is continuous everywhere.
We next present a sufficient condition to ensure that d * is Γ-minimax, so that d * ℓ is approximately Γ-minimax for sufficiently large ℓ.
Condition 3. We assume that there exists an increasing sequence {Ω ℓ } ∞ ℓ=1 of subsets of M such that 1.
We note that, in contrast to M ℓ , Ω ℓ may be an infinite set.We may expect Condition 3 to hold in many cases, especially when P → R(d, P ) is continuous at each d ∈ D and the grid M ℓ contains a variety of distributions that are consistent with prior information represented by Γ.We illustrate this point with two counterexamples in Appendix A. We will check the plausibility of Condition 3 for Example 2 in our simulation and data analysis in Section 5.1 for exemplar prior information; an almost identical argument shows the plausibility of Condition 3 for Example 3.
We now present the theorem on Γ-minimaxity of d * .
Theorem 1 (Validity of grid-based approximation).Under Conditions 1-3, d * is Γ-minimax and To prove Theorem 1, we utilize a result in Pinelis (2016) to establish that r sup (d, Γ) can be approximated arbitrarily well by a discrete prior in Γ for any d ∈ D. This is a key ingredient in the proof of Lemma 1, which states that, for any d ∈ D, r sup (d, Γℓ ) converges to r sup (d, Γ).Then, we show that the sequence {r sup (d * ℓ , Γ ℓ )} ∞ ℓ=1 is nondecreasing and upper bounded by inf d∈D r sup (d, Γ), which is less than or equal to the Γ-maximal Bayes risk r sup (d * , Γ) of the limit point d * of {d * ℓ } ∞ ℓ=1 in Condition 1.Therefore, r sup (d * ℓ , Γ ℓ ) converges to a limit.We finally use a contradiction argument to prove that this limit is greater than or equal to r sup (d * , Γ), which implies Theorem 1.
We have the following corollary on the uniqueness of the Γ-minimax estimator and the convergence of {d * ℓ } ∞ ℓ=1 for certain problems.is the unique Γ-minimax estimator and We prove Corollary 1 by establishing that d → r sup (d, Γ) is strictly convex.
In practice, the user also needs to specify a stopping criterion for Algorithm 1.In Noubiap and Seidel (2001), the authors recommended computing or approximating r sup (d * ℓ , Γ) and stop if r sup (d * ℓ , Γ) is sufficiently close to r sup (d * ℓ , Γ ℓ ).However, the procedure to approximate r sup (d * ℓ , Γ) in that work relies on the compactness of M, but we do not want to assume this condition because it may restrict the applicability of the method.Therefore, we propose to use the following alternative criterion: stop if for a prespecified tolerance level ϵ > 0. This criterion was proposed but not recommended in Noubiap and Seidel (2001) because it does not guarantee that ℓ is far from being Γ-minimax.In contrast, we recommend this criterion for our proposed methods because we allow more flexibility in model specification, that is, M need not be compact.We discuss this issue in more detail in Section 4.1.
We finally remark that r sup (d, Γ ℓ ) may be difficult to evaluate exactly.Since the risk is often an expectation, we recommend approximating r sup (d, Γ ℓ ) for any given d via Monte Carlo as follows: first, estimate risks R(d, P ) for all P ∈ M ℓ with a large number of Monte Carlo runs; second, estimate the corresponding least favorable prior π d,ℓ ∈ argmax π∈Γ ℓ r(d, π) using the estimated risks; third, estimate the risks R(d, P ) (P ∈ M ℓ ) again with independent Monte Carlo runs, and, finally, calculate r(d, π d,ℓ ) with the estimated risks and the estimated least favorable prior.Using two independent estimates of the risk can remove the positive bias that would otherwise arise due to using the same data to estimate the risks and the least favorable prior.

Computation of an estimator on a grid via stochastic gradient descent with max-oracle
In this section, we present methods to compute a Γ ℓ -minimax estimator, which corresponds to Line 3 in Algorithm 1. Gradient descent with max-oracle (GDmax) and its stochastic variant (SGDmax), which were presented in Lin et al. (2020), can be used to solve general minimax problems in Euclidean spaces.
We focus on SGDmax in the main text and present GDmax in Appendix B. To apply these algorithms to find a Γ ℓ -m inimax estimator, we need to assume that D can be parameterized by a subset of a Euclidean space, that is, that for any d ∈ D, there exists a real vector-valued coefficient β ∈ R D such that d may be written as d(β).For example, D may be a neural network class.More discussions on the parameterization of D can be found in Section 4.2.In this section, in a slight abuse of notation, we define R(β, P ) := R(d(β), P ), r(β, π) := r(d(β), π) and r sup (β, Γ ℓ ) := r sup (d(β), Γ ℓ ) for a coefficient β ∈ R D , a data-generating mechanism P ∈ M and a prior π ∈ Γ.We assume that β → R(β, P ) is differentiable for all P ∈ M, and hence so is β → r(β, π) for all π ∈ Γ.
It is often the case that R(β, P ) is expressed as an expectation.In this case, R(β, P ) may instead be approximated using Monte Carlo techniques.With ξ being an exogenous source of randomness according to law Ξ, let R(β, P, ξ) be an unbiased approximation of R(β, P ) with Stochastic maximization: use a stochastic procedure to find where the expectation is over the randomness in stochastic maximization (e.g., variants of stochastic gradient ascent).

5:
Stochastic gradient descent: We next present two conditions needed for the validity of Algorithm 2.
Note that Condition 4 differs from Condition 2 in that the former relies on the parameterization of D in a Euclidean space equipped with the Euclidean norm, while the latter may rely on a different metric on D such as an L 2 -distance.
Under these conditions, using the results in Lin et al. (2020), we can show that SGDmax yields an approximation to a local minimum of β → r sup (β, Γ ℓ ) when the algorithms' hyperparameters are suitably chosen.Before we formally present the theorem, we introduce some definitions related to the local optimality of a potentially nondifferentiable and nonconvex function.
x is an ϵ-stationary point in expectation, we may conclude that it is an ϵ-stationary point with high probability by Markov's inequality.Lemma 3.8 in Lin et al. (2020) shows that an ϵ-stationary point of f is close to a point x ′ at which f has at least one small subgradient for small ϵ, so that f (x ′ ) is close to a local minimum.In other words, if an algorithm outputs an estimator d = d( β) such that β is an ϵ-stationary point of β → r sup (β, Γ ℓ ), then we know that r sup ( β, Γ ℓ ) is close to a local minimum of We next present the validity result for Algorithm 2.
, and is thus close to a local minimum of β → r sup (β, Γ ℓ ) with high probability.
The assumption that the batch size J = 1 is purely for convenience since increasing J corresponds to decreasing variance σ 2 .To run Algorithm 2 in practice, the user only needs to specify tuning parameters in Line 1 and all other constants in Theorem 2 need not be known.In general, a small learning rate η, a stringent accuracy ζ, and a large batch size J make Algorithm 2 likely to eventually reach an approximation of a local minimum of β → r sup (β, Γ ℓ ), but computation time might increase.Similar to most numeric optimization algorithms, fine-tuning is needed to achieve a balance between convergence guarantee and computation time, but a conservative choice of tuning parameters would typically result in convergence at the cost of computation time.
We note that Line 3 in Algorithm 2 may be inconvenient to implement because linear program solvers often do not use stochastic optimization.Therefore, we propose to use a convenient variant (Algorithm 6 in Appendix B), where the stochastic maximization step (Line 3 in Algorithm 2) is replaced by solving a linear program where the objective is approximated via Monte Carlo.This variant has similar validity under similar conditions.We also note that the two uniform Lipschitz continuity conditions (4 and 5) heavily rely on the fact that M ℓ is finite and the compactness of a set containing the coefficients.
Nevertheless, the latter compactness restriction is common in theoretical analyses of neural networks (see, e.g., Goel et al., 2016;Zhang et al., 2016;Eckle and Schmidt-Hieber, 2019).Moreover, these two conditions are sufficient conditions for the validity of the gradient-based methods, namely SGDmax, our variant of SGDmax and GDmax; a guarantee similar to these validity results might hold when two conditions are violated.
We finally remark that other algorithms similar to SGDmax can be applied, for example, (stochastic) gradient descent ascent with projection (Lin et al., 2020), (stochastic) mirror descent ascent, or accelerated (stochastic) mirror descent ascent (Huang et al., 2021).It is of future research interest to develop gradient-based methods to solve minimax problems with convergence guarantees under weaker conditions.

Computation of an estimator on a grid via fictitious play
The algorithms in Section 3.2 may be convenient in many cases, but the requirements such as parameterization of the space D of estimators in a Euclidean space, differentiability of the risk function R with respect to the coefficients β, and uniform Lipschitz continuity may be restrictive for certain problems.
In this section, we propose an alternative algorithm, fictitious play, that avoids these requirements.We also present its convergence results.Brown (1951) introduced fictitious play as a means to find the value of a zero-sum game, that is, the optimal mixed strategy for both players and their expected gains.Robinson (1951) then proved that fictitious play can be used to iteratively solve a two-player zero-sum game for a saddle point that is a pair of mixed strategies where both players have finitely many pure strategies.Our problem of finding a Γ-minimax estimator may also be viewed as a two-player zero-sum game where one player chooses a prior from Γ and the other player chooses an estimator from D. If we assume that, for the Γ-minimax problem at hand, the pair of both players' optimal strategies is a saddle point, which holds in many minimax problems (e.g., v. Neumann, 1928;Fan, 1953;Sion, 1958), then fictitious play may also be used to find a Γ-minimax estimator.Since Γ may be too rich to allow for feasible implementation of fictitious play, we propose to use this algorithm to find a Γ ℓ -minimax estimator.
In the fictitious play algorithm in Robinson (1951), the two players take turns to play the best pure strategy against the mixture of the opponent's historic pure strategies, and the final output is a pair of mixtures of the two players' historic pure strategies.Since this algorithm aims to find minimax mixed strategies, we consider stochastic estimators.That is, consider the Borel σ-field F over D and let Π denote the set of all probability distributions on the measurable space (D, F).We define D to be the space of stochastic estimators with each element taking the following form: first draw an estimator from D according to a distribution ϖ ∈ Π with an exogenous random mechanism and then use the estimator to obtain an estimate based on the data.Note that we may write any d ∈ D as d(ϖ) for some ϖ ∈ Π.
We consider estimators in D throughout this section, with the definition of Γ-minimaxity extended in the natural way, so that d ; we similarly extend all other definitions from Section 2. We assume that there exists In other words, (d Algorithm 3 presents the fictitious play algorithm for finding a Γ ℓ -minimax estimator in D. Note that Γ ℓ is convex, and hence π always lies in Γ ℓ throughout the iterations.In practice, we may initialize ϖ as a point mass at an initial estimator in D. In addition, similarly to Robinson (1951), we may replace Line 5 with d † (t) ← argmin d∈D r(d, π (t) ), that is, minimizing the Bayes risk with the most recently updated prior rather than with the previous prior.
Algorithm 3 Fictitious play to compute a Γ ℓ -minimax stochastic estimator , where δ(d) denotes a point mass at d ∈ D.
We next present a convergence result for this algorithm.
Theorem 3 (Validity of fictitious play (Algorithm 3)).Assume that there exists a compact subset D of for all t and Consequently, the Γ ℓ -maximal risk of d(ϖ (t) ) converges to the Γ ℓ -minimax risk, that is, Robinson (1951) proved a similar case for two-player zero-sum games where each player has finitely many pure strategies.In contrast, in our problem, each player may have infinitely many pure strategies.A natural attempt to prove Theorem 3 would be to consider finite covers of D and Γ ℓ , that is, D = I i=1 D i and Γ ℓ = J j=1 Π j , such that the range of r(d, π) in each D i and Π j is small (say less than ϵ), bin pure strategies into these subsets, and then apply the argument in Robinson (1951) to these bins.The collection of D i and Π j may be viewed as finitely many approximated pure strategies to Γ ℓ and D up to accuracy ϵ, respectively.Unfortunately, we found that this approach fails.The problem arises because Robinson (1951) inducted on I and J, and, after each induction step, the corresponding upper bound becomes twice as large.Unlike the case with finitely many pure strategies that was considered in Brown (1951) and Robinson (1951), as the desired approximation accuracy ϵ approaches zero, the numbers of approximated pure strategies, I and J, may diverge to infinity, and so does the number of induction steps.Therefore, the resulting final upper bound is of order 2 I+J ϵ and generally does not converge to zero as ϵ tends to zero.To overcome this challenge, we instead control the increase in the relevant upper bound after each induction step more carefully so that the final upper bound converges to zero as ϵ decreases to zero, despite the fact that I and J may diverge to infinity.
We remark that, because Line 5 of Algorithm 3 typically involves another layer of iteration in addition to that over t, this algorithm will often be more computationally intensive than SGDmax.Nevertheless, Algorithm 3 provides an approach to construct Γ ℓ -minimax estimators in cases where these other algorithms cannot be applied, for example, in settings where the risk is not differentiable in the parameters indexing the estimator or uniform Lipschitz conditions fail.In our numerical experiments, we have implemented this algorithm in the context of mean estimation (Appendix C).

Considerations in implementation
4.1 Considerations when constructing the grid over the model space Though this guarantee holds for all such sequences {M ℓ } ∞ ℓ=1 , in practice, judiciously choosing this sequence of grids of distributions can lead to faster convergence.In particular, it is desirable that the least favorable prior Γ ℓ puts mass on some of the distributions in M ℓ \M ℓ−1 since, if this is not the case, then d * ℓ will be the same as d * ℓ−1 .While we may try to arrange for this to occur by adding many new points when enlarging M ℓ−1 to M ℓ , it may not be likely that any of these points will actually modify the least favorable prior unless they are carefully chosen.
To better address this issue, we propose to add grid points using a Markov chain Monte Carlo (MCMC) method.Our intuition is that, given an estimator d, the maximal Bayes risk is likely to significantly increase if we add distributions that (i) have a high risk for d, and (ii) are consistent with prior information so that there exists some prior such that these distributions lie in a high-probability region.We propose to use the MCMC algorithm to bias the selection of distributions in favor of those with the above characteristics.Let τ : M → [0, ∞) denote a function such that τ (P ) > τ (P ′ ) if P is more consistent with prior information than P ′ .For example, given a prior mean µ of some real-valued summary Ψ(P ) of P and an interval I that contains Ψ(P ) with prior probability at least 95%, we may choose τ : P → ϕ(Ψ(P )), where ϕ is the density of a normal distribution that has mean µ and places 95% of its probability mass in I.We call τ a pseudo-prior.Then, with the current estimator being d, we wish to select distributions P for which R(d, P )τ (P ) is large.We may use the Metropolis- Hastings-Green algorithm (Metropolis et al., 1953;Hastings, 1970;Green, 1995) to draw samples from a density proportional to P → R(d, P )τ (P ).We then let M ℓ be equal to the union of M ℓ−1 and the set containing all unique distributions in this sample.
Details of the proposed scheme are provided in Algorithm 4. To use this proposed algorithm, we rely on it being possible to define a sequence of parametric models { Ωℓ } ∞ ℓ=1 such that M := ∪ ∞ ℓ=1 Ωℓ is dense in M ℓ -this is possible in many interesting examples (see, e.g., Chen, 2007).When combined with separability of M, this condition enables the definition of an increasing sequence of grids of distributions {M ℓ } ∞ ℓ=1 such that, for each ℓ, M ℓ ⊆ M. The following theorem on distributional convergence follows from that for the Metropolis-Hastings-Green algorithm (see Section 3.2 and 3.3 of Green, 1995).Propose a distribution P ′ ∈ M from P (t−1)

5:
With probability p accept , accept P ′ and P (t) ← P ′ 6: if P ′ is not accepted then 7: Theorem 4 (Validity of MCMC algorithm (Algorithm 4)).Suppose that P → R(d * ℓ−1 , P )τ (P ) is bounded and integrable with respect to some measure µ on M and let L denote the probability law on M whose density function with respect to µ is proportional to this function.Suppose that the MCMC is constructed such that the Markov chain is irreducible and aperiodic.Then, P (t) converges weakly to L as t → ∞.
Therefore, if L corresponds to a continuous distribution with nonzero density over the parameter space of M, then Theorem 4 implies that ∞ ℓ=1 M ℓ is dense in M, as required by Algorithm 1. Implementing Algorithm 4 relies on the user making several decisions.These decisions include the choice of the pseudo-prior τ and the technique used to approximate the risk R(d, P ) to a reasonable accuracy.Fortunately, regardless of the decisions made, Theorem 1 suggests that r sup (d * ℓ , Γ ℓ ) ↗ min d∈D r sup (d, Γ) for a wide range of sequences {M ℓ } ∞ ℓ=1 .Indeed, all that theorem requires on this sequence is that the grid M ℓ becomes arbitrarily fine as ℓ increases.Though the final decisions made are not important when ℓ is large, we still comment briefly on the decisions that we have made in our experiments, First, we have found it effective to approximate R(d, P ) via a large number of Monte Carlo draws.Second, in a variety of settings, we have also identified, via numerical experiments, candidate pseudo-priors that balance high risk and consistency with prior information (see Sections 5.1 and 5.2 for details).

Considerations when choosing the space of estimators
It is desirable to consider a rich space D 0 of estimators to obtain an estimator with low maximal Bayes risk, and thus good general performance.However, to make numerically constructing these estimators computationally feasible, we usually have to consider a restricted space D of estimators.This approximation is justified because, if estimators in D can approximate the Gamma-minimax estimator in D 0 well, then we expect the resulting excess maximal Bayes risk is small.Feedforward neural networks (or neural networks for short) are natural options for the space of estimators because of their universal approximation property (e.g., Hornik, 1991;Csáji, 2001;Hanin and Sellke, 2017;Kidger and Lyons, 2020).However, training commonly used neural networks can be computationally intensive.Moreover, a space of neural networks is typically nonconvex, and hence it . . .may be difficult to find a global minimizer of the maximal Bayes risk even if the risk is convex in the estimator.Therefore, the learned estimator might not perform well.
To help overcome this challenge, we advocate for utilizing available statistical knowledge when designing the space of estimators.We call estimators that take this form statistical knowledge networks.
In particular, if a simple estimator is already available, we propose to use neural networks with such an estimator as a node connected to the output node.An example of such an architecture is presented in Fig. 2. In this sample architecture, each node is an activation function such as the sigmoid or the rectified linear unit (ReLU) (Glorot et al., 2011)  The dashed lines are the Bayes risks after an update in the prior to increase the Bayes risk.The two horizontal lines are the Bayes risk of the sample mean (dashed) and d * (solid), respectively, for π * .For ease of visualization, in subfigures (a) and (b), the Bayes risks are plotted every 50 iterations; in subfigures (c) and (d), the Bayes risks are plotted every 200 iterations; subfigure (d) contains the part in subfigure (c) after 500 iterations.
We note that we might overcome the challenge of nonconvexity and local optimality by using an extreme learning machine (ELM) (Huang et al., 2006b) to parameterize the estimator.ELMs are neural networks for which the weights in hidden nodes are randomly generated and are held fixed, and only the weights in the output layer are trained.Thus, the space of ELMs with a fixed architecture and fixed hidden layer weights is convex.Like traditional neural networks, ELMs have the universal approximation property (Huang et al., 2006a).In addition, Corollary 1 may be applied to an ELM so that the Γ ℓminimax estimator may converge to the Γ-minimax estimator.As for traditional neural networks, we may incorporate knowledge of existing statistical estimators into an ELM.
We finally remark that, besides computational intensity when constructing (i.e., learning) a Γminimax estimator, another important factor to be considered when choosing D is the computational intensity to evaluate the learned estimator at the observed dataset.This is another reason for our choosing neural networks or ELMs as the space of estimators.Indeed, existing software packages (e.g., Paszke et al., 2019) make it easy to leverage graphics processing units to efficiently evaluate the output of neural networks for any given input.Therefore, if the existing estimator being used is not too difficult to compute, then estimators parameterized using similar architectures to that displayed in Figure 2 will be able to be computed efficiently in practice.This efficiency may be especially important in settings where the estimator will be applied to many datasets, so that the cost of learning the estimator is amortized and the main computational expense is evaluating the learned estimator.

Simulations and data analyses
We illustrate our methods in Examples 1-3.A toy example of Example 1 is presented in Appendix C. We focus on the more complex Examples 2 and 3 in this section.

Prediction of the expected number of new categories
We apply our proposed method to Example 2. In the simulation, we set the true population to be an infinite population with the same categories and same proportions as the sample studied in Miller and Wiegert (1989), which consists of 1088 observations in 188 categories.This setting is the same as the simulation setting in Shen et al. (2003).We set the sample size to be n = 100 and the size of the new sample to be m = 200.In this setting, the expected number of new categories in the new sample unconditionally on the observed sample, namely Φ(P 0 ) := E P0 [Ψ(P 0 )(X * )], can be analytically computed and equals 48.02.We note that this quantity can also be computed via simulation: (i) sample n and m individuals with replacement from the dataset in Miller and Wiegert (1989), (ii) count the number of new categories in the second sample, and (iii) repeat steps (i) and (ii) many times and compute the average.
We note that a traditional Bayesian approach would require specifying a prior on M, including the total number of categories and the proportion of each category, which may be difficult in practice.
We check the plausibility of Condition 3 in this context.We take the strongly informative prior information as an example.Take Ω ℓ to be the collection of multinomial distributions with at most ℓ categories.It is obvious that Let d ∈ D be fixed and π ∈ Γℓ be a fixed prior with finite support, that is, π = J j=1 q j δ(Q j ) where δ(•) denotes the point mass distribution, Q j ∈ Ω ℓ , q j > 0 and J j=1 q j = 1.Let ϵ > 0 be an arbitrary small number such that is dense in M and Φ is continuous, there exists a sufficiently large i such that, for every distribution Q j , there exists P j ∈ M i ∩ Ω ℓ satisfying the following:  2016); SCL: the estimator proposed in Shen et al. (2003).The arrows from data (X 1 , . . ., X n ) to the OSW and SCL estimators are omitted from this graph.
Take π i to be J j=1 q j δ(P j ).Then it is easy to verify that 40,55] implies that Φ(P j ) ∈ [40, 55] and therefore We design the architecture of the neural network estimator as in Fig. 4. We choose two existing estimators (referred to as the OSW and SCL estimators, respectively) proposed by Orlitsky et al. (2016) and Shen et al. (2003) as human knowledge inputs to the architecture.We use the ReLU activation function.There are 50 hidden nodes in the first hidden layer.We initialize the neural network that we train to output the average of these two existing estimators.
We use Algorithm 4 to construct M ℓ .There are 2000 grid points in M 1 , and we add 1000 grid points each time we enlarge the grid.When generating M 1 , we chose the starting point to be a distribution P (0) with 146 categories and Φ(P (0) ) = 49.9.The choice of this starting point P (0) was quite arbitrary.
We first generated a sample from P 0 and treated it as data from a pilot study.We then came up with a distribution P (0) such that five random samples generated from P (0) all appear qualitatively similar to the pilot data.In practice, this starting point can be chosen based on prior knowledge.Our chosen grid sizes for Algorithm 4 were quite arbitrary.For M 1 , the generated distributions P (t) appear similar for all t, and thus the initial grid size 2000 and the increment size 1000 appeared sufficient.Smaller grid sizes would simply lead to more iterations in Algorithm 1, which effectively increases the grid size.We selected the log pseudo-prior as a weighted sum of two log density functions: (i) a normal distribution with the mean being the midpoint of the interval constraint on the prior mean of Φ(P ) and central 95% probability interval being the interval with at least 95% prior probability, (ii) a negative-binomial distribution of the total number of categories with success probability 0.995 and 2 failures until the Bernoulli trial is stopped so that the mode and the variance are approximately 200 and 8 × 10 4 , respectively.These logdensities are provided weights 30 and 10, respectively.We selected the weights based on the empirical observation that distributions with only a few categories tend to have high risks, but these distributions are relatively inconsistent with prior information and may well be given almost negligible probability weight in a computed least favorable prior, thus contributing little to computing a Γ-minimax estimator.
We chose the aforementioned weights so that Algorithm 4 can explore a fairly large range of distributions and does not generate too many distributions with too few categories.
We use Algorithm 6 with learning rate η = 0.005 and batch size J = 30 to compute Γ ℓ -minimax estimators.The number of iterations is 4,000 for Γ 1 and 200 for Γ ℓ (ℓ > 1).The stopping criterion in Algorithm 1 is that the estimated maximal Bayes risk with 2000 Monte Carlo runs does not relatively increase by more than 2% or absolutely increase by more than 0.0001.We chose the aforementioned tuning parameters based on the prior belief that at least one of OSW and SCL estimators should perform reasonably well, but the performance of SGDmax (Algorithm 6) and Algorithm 4 might be sensitive to tuning parameters.Thus, the network we used is neither deep nor wide.We chose a moderately small learning rate and a large number of iterations for SGDmax.Our chosen learning rate and chosen number of iterations led to a trajectory of estimated Bayes risks that approximately reached a plateau with small fluctuations, suggesting that the obtained estimator is approximately Γ 1 -minimax (see Fig. 5).
In practice, such trajectory plots may help tune the learning rate and the number of iterations.
We also run additional simulations to investigate the sensitivity of our methods to tuning parameter selections.We present these simulations in Appendix D. The results suggest that our methods may be insensitive to tuning parameter selections.
We examine the performance of the OSW estimator, the SCL estimator and our trained Γ-minimax estimator by comparing their risks under our set data-generating mechanism computed with 20000 Monte Carlo runs.We also compare their Bayes risks under the computed prior from Algorithm 6 using the last and finest grid in the computation with 20000 Monte Carlo runs.We present the results in Table 1.In this simulation experiment, our Γ-minimax estimator substantially reduces the risk compared to two existing estimators.The Γ-minimax estimator also has the lowest Bayes risk in all cases.Therefore, incorporating fairly informative prior knowledge into the estimator may lead to a significant improvement in predicting the number of new categories.We expect similar substantial improvement for difficult or even ill-posed statistical problems by incorporating prior knowledge.
Fig. 5 presents the unbiased estimator of Bayes risks over iterations when computing a Γ 1 -minimax estimator.The Bayes risks appear to have a decreasing trend and to approach a liming value.Over iterations, the Bayes risks decrease by a considerable amount.The limiting value of the Bayes risks appears to be slightly higher than the risk of the computed Γ-minimax estimator under P 0 .This might indicate that P 0 is not an extreme distribution that yields a high risk.
We also apply the above methods to analyze the dataset studied in Miller and Wiegert (1989), which is used as the true population in the simulation.Based on this sample consisting of n = 1088 observations in 188 categories, we use various methods to predict the number of new categories that would be observed if another m = 2000 observations were to be collected.We train Gamma-minimax estimators using Figure 5: Estimated Bayes risks of the estimator over iterations when computing a Γ 1 -minimax estimator.The lines are unbiased estimates of the current Bayes risks (y-axis) with 30 Monte Carlo runs over iterations (x-axis).The two dashed horizontal lines are the risks of the OSW (upper) and the SCL (lower) estimators, respectively, under P 0 in the simulation.The solid horizontal line is the risk of the computed Γ-minimax estimator under P 0 .For clearness of visualization, the estimated Bayes risks are plotted every 50 iterations.
exactly the same tuning parameters as those in the above simulation, except that the starting point in Algorithm 4 has more categories.The predictions of all methods are presented in Table 2.The Γ-minimax estimator outputs a more similar prediction to the SCL estimator.This similarity appears different from our observation in the simulation, but can be explained by the fact that having more observations (n = 1088 vs n = 100; m = 2000 vs m = 200) decreases the variance of the number of new observed categories and thus lowers discrepancies between predictions from these methods.Since the SCL estimator outperforms the OSW estimator in the above simulation where this dataset is the true population, we expect the SCL estimator to achieve reasonably good performance in this application.
Moreover, given that the Γ-minimax estimators outperform the SCL estimator in the above simulation, we expect that 57 or 58 represents an improved prediction of the number of new categories as compared to the SCL prediction of 51 in the case where there is limited prior information available.
The computation time to compute an approximated Γ-minimax estimator was about five to seven  Miller and Wiegert (1989).The strength of prior information in Γ-minimax estimators is shown in brackets.

Estimator
Predicted # new categories OSW 72 SCL 51 Γ-minimax (strong) 57 Γ-minimax (weak) 57 Γ-minimax (almost none) 58 hours on an AWS EC2 instance (Amazon, 2019) with at least 4 vCPUs and at least 8 GiB of memory, depending on the number of times the grid was enlarged.As shown in Fig. 5, far few iterations are needed for SGDmax to output a good approximation of a Γ 1 -minimax estimator, which is itself quite close to Γ-minimax.Therefore, with suitably less conservative tuning parameters or more adaptive minimax problem solvers, the computation time might drastically decrease.Moreover, the computation time needed to evaluate the computed estimator at any sample is almost zero.

Estimation of the entropy
We also apply our method to estimate the entropy of a multinomial distribution (Example 3).The data-generating mechanism is the same as that described in Example 2, and the estimand of interest is Shannon entropy (Shannon, 1948), that is, Ψ(P 0 ) = − K k=1 p k log p k .In the simulation, we choose the same true population and the same sample size n = 100 as in Section 5.1.The true entropy Ψ(P 0 ) is 4.57.As a reference, the entropy of the uniform distribution with the same number of categorieswhich corresponds to the maximum entropy of multinomial distributions with the same total number of categories-is 5.24.Jiao et al. (2015) developed a minimax rate optimal estimator of the Shannon entropy, and we run this simulation to investigate the potential gain of computing a Gamma-minimax estimator in well-posed problems with satisfactory solutions.As in Section 5.1, we consider three sets of prior information: 1. Strongly informative: Prior mean of Ψ(P ) in [4.3, 4.7], ≥ 95% probability that Ψ(P ) lies in [4,5]; 2. Weakly informative: Prior mean of Ψ(P ) in [4,5], ≥ 95% probability that Ψ(P ) lies in [3.5, 5.5]; 3. Almost noninformative: Prior mean of Ψ(P ) in [3.7, 5.3], ≥ 95% probability that Ψ(P ) lies in [3,6].
The architecture of our neural network estimator is almost identical to that in Section 5.1 except that the existing estimator being used is the one proposed in Jiao et al. (2015) (referred to as the JVHW estimator), and we initialize the network to return the JVHW estimator.We use Algorithm 4 to construct M ℓ and Algorithm 6 to compute a Γ ℓ -minimax estimator.The tuning parameters in the algorithms are identical to those used in Section 5.1 except that, in Algorithm 6, (i) the learning rate is η = 0.001, and (ii) the number of iterations is 6,000 for Γ 1 .We change these tuning parameters because the JVHW estimator is already minimax in terms of its convergence rate (Jiao et al., 2015), and we may need to update the estimator in a more cautious manner in Algorithm 6 to obtain any possible improvement.The trajectories of the estimated Bayes risks (Fig. 6) all appear to approximately reach a plateau, suggesting that the obtained estimator approximately Γ 1 -minimax and that our choice of a smaller learning rate and a larger number of iterations is valid.Because of the additional complexity of the JVHW estimator, we ran our simulations on an AWS EC2 instance (Amazon, 2019) with 4 vCPUs and 32 GiB of memory.The computation time was ten to seventeen hours, depending on the number of times the grid was enlarged.The longer computation time than that described in Section 5.1 is primarily due to more iterations in SGDmax and the additional complexity of the JVHW estimator.
We compare the risk of the JVHW estimator and our trained Γ-minimax estimator under our set data-generating mechanism computed with 20000 Monte Carlo runs.We also compare their Bayes risk under the computed prior from Algorithm 6 using the last and finest grid in the computation with 20000 Monte Carlo runs.The results are summarized in Table 3.In this simulation experiment, our Γ-minimax estimator reduces the risk by a fair percentage compared with the JVHW estimator and achieves lower worst-case Bayes risk.According to these simulation results, we conclude that incorporating informative prior knowledge into the estimator may result in some improvement in estimating entropy.Thus, for well-posed statistical problems with satisfactory solutions, we expect mild or no substantial improvement and little deterioration from using a Gamma-minimax estimator.
Fig. 6 presents the unbiased estimator of Bayes risks over iterations when computing a Γ 1 -minimax estimator.With strongly informative prior information present, the Bayes risks appear to fluctuate without an increasing or decreasing trend at the beginning and decrease slowly after several thousand iterations.With weakly informative or almost no prior information, the Bayes risks also decrease slowly.
A reason may be that the JVHW estimator is already minimax rate optimal (Jiao et al., 2015).The computed Γ-minimax estimators also appear to be somewhat similar to the JVHW estimator: in the output layer of the three settings with different prior information, the coefficients for the JVHW estimator are 0.97, 0.90 and 0.89, respectively; the coefficients for the previous hidden layer are 0.17, 0.17 and 0.20, respectively; the intercepts are 0.06, 0.30 and 0.30, respectively.
We further use the above methods to estimate entropy based on the dataset used as the true pop- Figure 6: Estimated Bayes risks of the estimator over iterations when computing a Γ 1 -minimax estimator.The lines are unbiased estimates of the current Bayes risks (y-axis) with 30 Monte Carlo runs over iterations (x-axis).The horizontal lines are the risks of the JVHW (dashed) and the computed Γ-minimax (solid) estimators, respectively, under P 0 in the simulation.For clearness of visualization, the estimated Bayes risks are plotted every 100 iterations.(Jiao et al., 2015), we expect the JVHW estimator to have little room for improvement, which explains why the three Γ-minimax estimators perform similarly to the JVHW estimator.In other words, Gamma-minimax estimators appear to maintain, if not improve, the performance of the original JVHW estimator.

Discussion
We propose adversarial meta-learning algorithms to compute a Gamma-minimax estimator with theoretical guarantees under fairly general settings.These algorithms still leave room for improvement.As we discussed in Section 3.1, the stopping criterion we employ does not necessarily indicate that the maximal Bayes risk is close to the true minimax Bayes risk.In future work, it would be interesting to derive a better criterion that necessarily does indicate this near optimality.Our algorithms also require the user to choose increasingly fine approximating grids to the model space.Although we propose a heuristic algorithm for this procedure that performed well in our experiments, at this point, we have not provided optimality guarantees for this scheme.It may also be possible to improve our proposed algorithms to solve intermediate minimax problems in Section 3.1 by utilizing recent and ongoing advances from the machine learning literature that can be used to improve the training of generative adversarial networks.
We do not explicitly consider uncertainty quantification such as confidence intervals or credible intervals under a Gamma-minimax framework.Uncertainty quantification is important in practice since it provides more information than a point estimator and can be used for decision-making.In theory, our method may be directly applied if such a problem can be formulated into a Gamma-minimax problem.However, such a formulation remains unclear.The most challenging part is to identify a suitable risk function that correctly balances the level of uncertainty and the size of the output interval/region.
Though the risk function used in Schafer and Stark (2009) appears to provide one possible starting point, it is not clear how to extend this approach to nonparametric settings.
It is possible to allow the space of estimators to increase as the grid M ℓ increase.For example, we may specify an increasing sequence of estimator spaces {D ℓ } ∞ ℓ=1 whose limit is dense in a general space D 0 ; then, in Line 3 of Algorithm 1, we compute a Γ ℓ -minimax estimator in D ℓ , namely replace D with D ℓ .This sequence of estimators might converge to a Γ-minimax estimator in D 0 .One possible choice of D ℓ (ℓ > 1) in this approach is a space of statistical knowledge networks with the given estimator being the computed Γ ℓ−1 -minimax estimator in D ℓ−1 .It is of future interest to investigate the properties of such an approach.
In conclusion, we propose adversarial meta-learning algorithms to compute a Gamma-minimax estimator under general models that can incorporate prior information in the form of generalized moment conditions.They can be useful when a parametric model is undesirable, semi-parametric efficiency theory does not apply, or we wish to utilize prior information to improve estimation.

A Two counterexamples of Condition 3
We provide two counterexamples of Condition 3 to illustrate that this condition fails in extremely ill cases.
In the first counterexample, P → R(d, P ) is discontinuous: we set R(d, P * ) to be zero for a fixed P * ∈ M and R(d, P ) to be one for all other P ∈ M. If we choose the grid M ℓ to be dense in M but to never contain P * , then Condition 3 does not hold since r sup (d, Γℓ ) = 1 for sufficiently large ℓ such that P * ∈ Ω ℓ but r sup (d, Γi|ℓ ) = 0 for all i and ℓ.This issue can be resolved by choosing a continuous risk function.
In the second counterexample, M ℓ does not contain distributions that are consistent with prior information.Suppose that Γ = {π ∈ Π : Φ(P )π(dP ) = 0} where Φ(P ) := E P [X 2 ].In other words, it is known that the true data-generating mechanism P 0 must be a distribution that is a point mass at zero, and thus Γ also only contains a point mass at P 0 .If Φ(P ) ̸ = 0 for every P ∈ ∪ ∞ i=1 M i , then, even if ∞ ℓ=1 M ℓ is dense in M, Γi|ℓ = ∅ and thus Condition 3 does not hold.This issue can be resolved by rewriting the problem such that these hard constraints on M are incorporated into the specification of M rather than Γ.

B Additional gradient-based algorithms
If we can evaluate R(β, P ) exactly for all β ∈ H and P ∈ M ℓ , then the GDmax algorithm (Algorithm 5) may be used.Note that Line 3 can be formulated into a linear program, which can always be solved in polynomial time with an interior point method (e.g., Jiang et al., 2020) and often be solved in polynomial time with a simplex method (Spielman and Teng, 2004).
Theorem 5 (Validity of GDmax (Algorithm 5)).Under conditions in Theorem 2, in Algorithm 5, with Therefore, we propose a variant (Algorithm 6) by replacing this line with Lines 3-4 so that ordinary linear program solvers can be directly applied.The following theorem justifies this variant.
The validity of this variant of SGDmax is given in Theorem 6 below.

6:
Stochastic gradient descent: there exists a sufficiently large J ′ such that where the expectation is taken over π (t) and β (t−1) is fixed.Therefore, with the chosen parameters in Theorem 2, we may choose a sufficiently large J ′ so that ) and is thus close to a local minimum of β → r sup (β, Γ ℓ ) with high probability.
We prove Theorem 6 by showing that max π∈Γ ℓ r( The proof is essentially an application of empirical process theory to the study of an M-estimator.

C Additional simulation: mean estimation
In this appendix, we illustrate our proposed methods via simulation in a special case of Example 1, namely for estimating the mean of a distribution.We assume that M consists of all probability distributions defined on the Borel σ-algebra on [0, 1] and we observe X = (X 1 , X 2 , . . ., X n ), where X 1 , . . ., X n iid ∼ P 0 ∈ M.Here we take n = 10.The estimand is Ψ(P 0 ) = x P 0 (dx).We use the mean squared error risk introduced in Example 1. Suppose that we represent the prior information by Γ = {π ∈ Π : Ψ(P ) π(dP ) = 0.3}, which corresponds to the set of prior distributions in Π that satisfy an equality constraint on the prior mean of Ψ(P ).
We apply our method to three spaces of estimators separately.The first space, D linear , is the set of affine transformations of the sample mean, that is, As shown in Proposition 1 in Appendix E.5, there is an estimator d * in D linear that is Γ-minimax in the space of all estimators that are square-integrable with respect to all P ∈ M, so we consider this simple space to better compare our computed estimator with that theoretical Γ-minimax estimator.When computing a Γ-minimax estimator in D linear , we initialize the estimator to be the sample mean, that is, we let β 0 = 0 and β 1 = 1.
The second space, D skn (statistical knowledge network), is a set of neural networks designed based on statistical knowledge that includes the sample mean as an input.We consider this space to illustrate our proposal in Section 4.2.More precisely, we use the architecture in Fig. 7, which is similar to the deep set architecture (Zaheer et al., 2017;Maron et al., 2019) and is a permutation invariant neural . . . . . . . . . Figure 7: Architecture of the permutation invariant neural network estimator of the mean in D skn .X i : observation i in the sample; : the node that sums up all ancestor nodes.In the first two hidden layers, all input nodes are transformed by the same function.The arrows from the input nodes to the sample mean estimator are omitted from this graph.Each node in the hidden layers represents a vector.
network.We use such an architecture to account for the fact that the sample is iid.In this architecture, the sample mean node is used as an augmenting node to an ordinary deep set network and is combined with the output of that ordinary network in the fourth hidden layer to obtain the final output.Note that D skn ⊃ D linear .When computing a Γ-minimax estimator for this class, we also initialize the network to be exactly the sample mean, which is a reasonable choice given that the sample mean is known to be a sensible estimator.In this simulation experiment, we choose the dimensionality of nodes in each hidden layer in Fig. 7 as follows: each node in the first, second, third and fourth hidden layer represents a vector in R 10 , R 5 , R 10 and R, respectively.We do not use larger architectures because usually the sample mean is already a good estimator, and we expect to obtain a useful estimator as a small perturbation of this estimator.We also use the ReLU as the activation function.We did not use ELMs in this and the following simulations because we found that neural networks perform well.
The third space, D nn , is a set of neural networks that do not utilize knowledge of the sample mean.
We consider this space to illustrate our method without utilizing existing estimators.These estimators are also deep set networks with similar architecture as D skn in Fig. 7.The main difference is that the explicit sample mean node and the fourth hidden layer are removed.When computing a Γ-minimax estimator in D nn , we also randomly initialize the network, unlike D linear and D skn , in order not to input statistical knowledge.Because the ReLU activation function is used, D nn ⊃ D linear , and we do not expect that optimizing over D nn should not lead to a Γ-minimax estimator with worse performance than those in D linear and D skn .
To construct the grid M ℓ for this problem, we use a simpler method than Algorithm 4. As indicated by Lemma 6 in Appendix E.5, for estimators in D linear , Bernoulli distributions tend to have high risks since all probability weights lie on the boundary of [0, 1]; in addition, a prior π * for which d * is Bayes is a Beta prior over Bernoulli distributions.Therefore, we randomly generate 2000 Bernoulli distributions as grid points in M 1 .We also include two degenerate distributions in this grid, namely the distribution that places all of its mass at 0 and that which places all of its mass at 1.When constructing M ℓ from M ℓ−1 , we still add in more complicated distributions to make the grid dense in the limit: we first randomly generate 500 discrete distributions with support being those in M ℓ−1 ; then we randomly generate 10 new support points in [0, 1] and 1000 distributions with support points being the union of the new support points and the existing support points in M ℓ−1 .
When computing the Γ-minimax estimator, for each grid M ℓ , we compute the Γ ℓ -minimax estimator for all three estimator spaces with Algorithm 6.We set the learning rate η = 0.005, the batch size J = 50 and the number of iterations to be 200 for Γ ℓ (ℓ > 1).The number of iterations for Γ 1 is larger because, in our experiments, we saw that a Γ 1 -minimax estimator is already close to a Γ-minimax estimator, and using a large number of iterations in this step can improve the initial estimator substantially.For D linear and D skn , the number of iterations for Γ 1 is 2000; the corresponding number for D nn is 6000 to account for the lack of human knowledge input.We also use Algorithm 3 with 10000 iterations to compute a Γ ℓ -minimax estimator for D linear for illustration.In this setup, as described in Section 3.3, we take the average of the computed Γ-minimax stochastic estimator as the final output estimator in D linear .We do not apply Algorithm 3 to D skn or D nn because it is computationally intractable for these estimator spaces.
We set the stopping criterion in Algorithm 1 as follows.When Algorithm 6 is used to compute Γ ℓminimax estimators, we estimate both r sup (d * ℓ−1 , Γ ℓ ) and r sup (d * ℓ−1 , Γ ℓ−1 ) with 2000 Monte Carlo runs as described in Section 3.1; when Algorithm 3 is used, r sup (d * ℓ−1 , Γ ℓ ) and r sup (d * ℓ−1 , Γ ℓ−1 ) are computed exactly because R(d, P ) has a closed-form expression for all d ∈ D linear and P ∈ M ℓ .We set the tolerance ϵ to be equal to 0.0001 so that we stop Algorithm 1 if r sup (d * ℓ−1 , Γ ℓ ) − r sup (d * ℓ−1 , Γ ℓ−1 ) ≤ ϵ.After computation, we report the Bayes risk of the computed and theoretical Γ-minimax estimators under π * , the prior such that r(d * , π * ) = inf d∈D r sup (d, Γ).For the estimators in D linear , we further report their coefficients.We also report two coefficients of the computed estimator in D skn as follows.
Since D linear ⊆ D skn and we initialize the estimator to be the sample mean for D skn , we would expect that the bias β 0 and the weight of the sample mean β 1 in the output layer for the computed Γ-minimax estimator in D skn may correspond to those in D linear .Therefore, we also report these two coefficients β 0 and β 1 for D skn .This may not be the case for D nn because the sample mean is not explicit in its parameterization and all coefficients are randomly initialized, so we do not report any coefficients for D nn .
Table 5 presents the computation results.By Theorem 7 in Appendix E.5, these computed estimators are all approximately Γ-minimax since their Bayes risks for π * are all close to that of a theoretical Γminimax estimator.The coefficients β 0 and β 1 of the computed estimators in D linear and D skn are also close to a theoretically derived estimator.For the computed estimator in D skn , the weight of the other ancestor node in the output layer (i.e., the node in the 4th hidden layer in Fig. 7) is 0.000.Therefore, our computed Γ-minimax estimator in D skn is also close to a theoretically derived Γ-minimax estimator.
In our experiments, Algorithm 1 converged after computing a Γ 1 -minimax estimator except when ). Fig. 3 presents the Bayes risks (or its unbiased estimates) over iterations when computing a Γ 1minimax estimator.In all cases using Algorithm 6, the Bayes risks appear to decrease and converge.
When using Algorithm 3, the upper and lower bounds both converge to the same limit.The limiting values of the Bayes risks in all cases are close to r(d * , π * ) because Γ 1 can approximate π * well.

D Sensitivity analysis for tuning parameter selection
For the simulation in Section 5.1 with strongly informative prior information, we conduct three simulations to investigate the sensitivity of our proposed method to the selection of tuning parameters.In each simulation below, we vary one set of tuning parameters and rerun the algorithm to obtain an estimator.
In the first simulation, we vary the starting point of Algorithm 4 to construct the first grid M 1 .The new starting point is a distribution with 173 categories and Φ(P (0) ) = 61, and so this starting point is qualitatively different from the one chosen in the original simulation.In the second simulation, we vary the grid sizes: There are 500 grid points in M 1 and we add 500 grid points each time we enlarge the grid.In the third simulation, we chose a wider and deeper statistical knowledge network (see Fig. 8): Compared to the original simulation, we add one more hidden layer and increased the number of hidden nodes in the first two hidden layers to 100.As shown in Table 6, the results in these sensitivity simulations appear similar to that in Section 5.1 within the variation due to randomness in MCMC (Algorithm 4) and SGDmax (Algorithm 6).Moreover, there exists π ℓ ∈ Γ ℓ such that Using the fact that d * ℓ is Γ ℓ -minimax and the definition of r sup , it holds that Since this inequality holds for any ϵ > 0, we have that r sup (d * ℓ , Γ ℓ ) ≤ inf d∈D r sup (d, Γ).An almost identical argument shows that the sequence {r sup (d * ℓ , Γ ℓ )} ∞ ℓ=1 is nondecreasing.Therefore, this sequence converges to some limit R ≤ inf d∈D r sup (d, Γ) ≤ r sup (d * , Γ).
Proof of Corollary 1.We first establish the strict convexity of d → r(d, π) for any π ∈ Γ.We then establish the strict convexity of d → r sup (d, Γ).We then establish that there is a unique minimizer of d → r sup (d, Γ) and show that the desired result follows from Theorem 1.
Let d 1 , d 2 ∈ D and c ∈ (0, 1) be arbitrary, then by the convexity of D and the strict convexity of Therefore, d → r(d, π) is strictly convex for any π ∈ Γ.
Let d 1 , d 2 ∈ D be distinct and c ∈ (0, 1) be arbitrary.Since r sup (d, Γ) is attainable for any d ∈ D, there exists π ∈ Γ such that ) is strictly convex and D is convex, this function achieves exactly one minimum on D. By Theorem 1, any limit point d * of {d * ℓ } ∞ ℓ=1 is a minimizer of d → r sup (d, Γ), and so the sequence has a limit point, which is also the unique Γ-minimax estimator.

E.2 Proof of Theorems 2 & 5
We prove Theorems 2 and 5 by checking that Assumptions 3.1 and 3.6 in Lin et al. (2020) are satisfied and using Theorem E.3 and E.4 in Lin et al. (2020), respectively.Since Assumption 3.1 is satisfied by our construction of R, we focus on Assumption 3.6 for the rest of this section.
For each ℓ = 1, 2, . .., for any β 1 , β 2 ∈ H and π ∈ Γ ℓ , Condition 4 implies that Finally, it is straightforward to check that (i) π → r(β, π) is concave for any β ∈ H, and (ii) Γ ℓ is parameterized by a convex subset of a simplex in a Euclidean space, which is a convex and bounded set.
These results show that Assumption 3.6 in Lin et al. ( 2020) is satisfied for Algorithm 5 and 2.

E.3 Proof of Theorem 6
Proof of Theorem 6.Let π (t),0 denote a maximizer of π → r(β (t−1) , π).It holds that Note that the right hand side does not depend on t.Therefore, where E * stands for outer expectation.We may apply Corollary 9.27 in Kosorok (2008) to Here, X * stands for the minimal measurable majorant with respect to Ξ of a (possibly non-measurable) mapping X (van der Vaart and Wellner, 2000).

E.4 Proof of Theorem 3
Our proof of Theorem 3 builds on that of Robinson (1951).Major modifications are needed to allow for more general definitions that can accommodate for potentially infinite spaces of pure strategies and a more careful control on a bound on r(d(ϖ ) towards the end of the proof.In this appendix, we slightly abuse the notation and use D to denote the compact set D that contains all d † (t) (t = 1, 2, . ..).We first introduce the notion of cumulative Bayes risk functions.Under Algorithm 3, we let U 0 : D → R and V 0 : Γ ℓ → R be any two continuous functions such that and recursively define for d ∈ D and π ∈ Γ ℓ .Here, we let π † (t) ∈ argmax π∈Γ ℓ V t−1 (π) and d † (t) ∈ argmin d∈D U t−1 (d).Note that the choices of π † (t) and d (t) in Algorithm 3 corresponds to setting U 0 ≡ 0 and V 0 ≡ 0, in which case for some π (t) ∈ Γ and d(ϖ (t) ) ∈ D. Later in this section, we will also make use of U t and V t with other initializations U 0 and V 0 .
To make notations concise, we define min d∈D ′ U t := min d∈D ′ U t (d) for any D ′ ⊆ D, and define max D ′ U t , min Π ′ V t and max Π ′ V t (Π ′ ⊆ Γ ℓ ) similarly.We also drop the subscript denoting the set when the set is the whole space we consider, that is, D or Γ ℓ .Note that for any t 1 , t 2 = 1, 2, . .., under the setting of Algorithm 3 and (2), it holds that Therefore, to prove the first result in Theorem 3, it suffices to show that lim sup t→∞ (max V t −min U t )/t ≤ 0.
We next introduce additional definitions related to iterations.We say that π ∈ Γ ℓ is eligible in the Note that by Condition 1 and Lemma 2, given an arbitrary desired approximation accuracy ϵ > 0, D can be covered by finitely many compact subsets with the maximum variation of each subset bounded by ϵt for all t; by Condition 2, since Γ ℓ is parameterized by a compact subset of a simplex in a Euclidean space, Γ ℓ can also be covered by finitely many compact subsets with the maximum variation of each subset bounded by ϵt for all t.These covers can be viewed as discrete finite approximations to D and Γ ℓ , respectively.
All of the above definitions are associated with the space of estimators D and the set of priors Γ ℓ .We call {(U t , V t )} t a pair of cumulative Bayes risk functions constructed from the pair (D, Γ ℓ ) of the space of estimators and the set of priors, and will consider pairs of cumulative Bayes risk functions constructed from other pairs (D ′ , Π ′ ) of the space of estimators and the set of priors in the subsequent proof.We can define the above quantities similarly for such cases.
The following lemma gives an upper bound on the maximum variation of U s+t and V s+t over the corresponding entire space from which they are constructed after t iterations from s when essentially all parts of these spaces are eligible in [s, s + t].
Lemma 4. Suppose that {(U t , V t )} t is a pair of cumulative Bayes risk functions constructed from shows that max Π ′ V s+t − min Π ′ V s+t ≤ (2B + A)t.
The next lemma builds on the previous lemma and provides an upper bound on max V s+t − min U s+t under the same conditions.
Lemma 5.Under the same setup and conditions as in Lemma 4, max Π ′ V s+t −min D ′ U s+t ≤ (4B +2A)t.
Proof of Lemma 5. Summing the two inequalities in Lemma 4 and rearranging the terms, we have that max Π ′ V s+t − min D ′ U s+t ≤ (4B + 2A)t + min Π ′ V s+t − max D ′ U s+t .It therefore suffices to show that min Π ′ V s+t ≤ max D ′ U s+t .
Let τ := s + t.There exists π ′ ∈ Π ′ and a stochastic strategy d ′ ∈ D ′ such that U τ (d) = U 0 (d) + τ • r(d, π ′ ) and V τ (π) = V 0 (π) + τ • r(d ′ , π) for all d ∈ D ′ and all π ∈ Π ′ .Therefore, for this choice of π ′ and d ′ , using (3), min Proof of Theorem 3. It suffices to show that lim sup t→∞ (max V t − min U t )/t ≤ 0 by letting U 0 ≡ 0 and V 0 ≡ 0, which corresponds to Algorithm 3. Let ϵ > 0. Note that r is Lipschitz continuous by Lemma 2 and the fact that r(d, π) is an average of bounded risks with weights π.Furthermore, D and Γ ℓ are both compact.In addition, U 0 and V 0 are both continuous.Therefore, there exist covers D = I i=1 D i and Γ ℓ = J j=1 Π j such that (i) D i and Π j are all compact, and (ii) sup i,t MV t (D i )/t ≤ ϵ, sup j,t MV t (Π j )/t ≤ ϵ. (Note that I and J may depend on ϵ.)For index sets I ⊆ {1, 2, . . ., I} and J ⊆ {1, 2, . . ., J}, define D I := i∈I D i and Π J := j∈J Π j .We show that max V t − min U t ≤ Cϵt for an absolute constant C and all sufficiently large t via induction on the sizes of I and J .
Let {(U t , V t )} t be a pair of cumulative Bayes risk functions constructed from (D I , Π J ) where |I| = |J | = 1.By (5) and the fact that MV t (D I ) ≤ ϵt and MV t (Π J ) ≤ ϵt, we have that Therefore, max Π J V t − min D I U t ≤ 2ϵt.
Let ϵ ′ > 0 be arbitrary.Suppose that there exists t 0 such that, for any I ′ ⊆ I and J ′ ⊆ J such that I ′ ̸ = I or J ′ ̸ = J , for any pair of cumulative Bayes risk functions {(U t , V t )} t constructed from (D I ′ , Π J ′ ), it holds that max Π J ′ V t − min D I ′ U t ≤ ϵ ′ t for all t ≥ t 0 .We next obtain a slightly greater bound on max Π J V t − min D I U t for all sufficiently large t.
We first prove that if, for a given pair of cumulative Bayes risk functions {(U t , V t )} t constructed from (D I , Π J ), there exists i ′ ∈ I or j ′ ∈ J such that D i ′ or Π j ′ is not eligible in an interval [s, s + t 0 ], then max Suppose that D i ′ is not eligible in [s, s + t 0 ], then define U ′ t := U s+t and V ′ t := V s+t − max Π J V s + min D I U s for all t ≥ 0. It is straightforward to check that {(U ′ t , V ′ t )} t0 t=0 satisfies the recursive definition of a pair of cumulative Bayes risk functions constructed from (D I\{i ′ } , Π J ).By the induction hypothesis, max Π J V ′ t0 −min D I\{i ′ } U ′ t0 ≤ ϵ ′ t 0 .Therefore, max Π J V s+t0 −min D I U s+t0 = max Π J V ′ t0 −min D I\{i ′ } U ′ t0 + max Π J V s − min D I U s ≤ max Π J V s − min D I U s + ϵ ′ t 0 .Similar argument can be applied if Π j ′ is not eligible in [s, s + t 0 ].Now we obtain a bound on max Π J V t −min D I U t .Let t > t 0 , Q := ⌊t/t 0 ⌋ ≥ 1 and R := t/t 0 −Q ∈ [0, 1).
There are two cases.Case 1: There exists s 0 ≤ Q such that D i and Π j are eligible in [(s 0 − 1 + R)t 0 , (s 0 + R)t 0 ] for all i ∈ I and j ∈ J .Take s 0 to be the largest such integer.Then, repeatedly apply (6) to intervals [(s 0 + R)t 0 , (s 0 + 1 + R)t 0 ], [(s 0 + 1 + R)t 0 , (s 0 + 2 + R)t 0 ], . . ., [(Q − 1 + R)t 0 , (Q + R)t 0 ] = [t − t 0 , t] and we derive that max By Lemma 5, max Π J V (s0+R)t0 − min D I U (s0+R)t0 ≤ (4B + ϵ)t 0 .Therefore, max Case 2: There is no integer s 0 satisfying the condition in Case 1.Then, repeatedly apply (6) to intervals [Rt 0 , (1 + R)t 0 ], [(1 + R)t 0 , (2 + R)t 0 ], . . ., [(Q − 1 + R)t 0 , (Q + R)t 0 ] = [t − t 0 , t], we derive that max By the bound on the risk, max Π J V Rt0 ≤ BRt 0 and min D I U Rt0 ≥ −BRt 0 .Hence, max Thus, in both cases, it holds that max Π J V t − min D I U t ≤ (4B + ϵ)t 0 + ϵ ′ t for t > t 0 .Let C > 0 be any constant (which may depend on ϵ, the approximation error of the covers, that is, the bound on MV t /t).The following holds for any sufficiently large t, max In other words, we show that after increasing the size of either index set by 1, for all sufficiently large t, we obtain a bound on max Π J V t − min D I U t that grows by a multiplicative factor of (1 + C) relative to the original bound.
It takes finitely many, say N , steps to induct from the initial case where the sizes of both index sets are one to the case of interest with index sets {1, . . ., I} and {1, . . ., J}. (Note that N may also depend on ϵ through its dependence on I and J.) Take C = 1/N in (7) and we derive that, for all sufficiently where e is the base of natural logarithm.Since ϵ is arbitrary, we show that lim sup t→∞ (max V t − min U t )/t ≤ 0.

E.5 Derivation of Γ-minimax estimator of the mean in Section C
In this section, we show that, for the problem of estimating the mean in Section C, one Γ-minimax estimator lies in D linear .This is formally presented below.
Proposition 1.Let M consist of all probability distributions defined on the Borel σ-algebra on [0, 1].
Let X 1 , . . ., X n iid ∼ P 0 ∈ M and X = (X 1 , X 2 , . . ., X n ) be the observed data.We first present a theorem on a criterion of Γ-minimaxity.

Corollary 1 (
Convergence of Γ ℓ -minimax estimator).Suppose that D is a convex subset of a vector space, d → R(d, P ) is strictly convex for each P ∈ M, and r sup (d, Γ) is attainable for each d ∈ D in the sense that, for all d ∈ D, there exists a π ∈ Γ such that r(d, π) = r sup (d, Γ).Under Conditions 1-3, d * Under this condition and the further conditions that D is convex and d → R(d, P ) is convex for all P ∈ M, it is possible to use a Γ-minimax estimator over the richer class D of stochastic estimators to derive a Γ-minimax estimator over the original class D. Indeed, for any d(ϖ) ∈ D and P ∈ M, by Jensen's inequality, R(d(ϖ), P ) = R(d, P ) ϖ(dd) ≥ R(d(ϖ), P ) where d(ϖ) := d ϖ(dd) ∈ D is the average of the stochastic estimator d(ϖ); that is, the risk of d(ϖ) is never greater than that of d(ϖ).Therefore, we may use the fictitious play algorithm to compute d(ϖ * ℓ ) for each ℓ and further apply Algorithm 1 to compute d(ϖ * ).After that, we may take d(ϖ * ) as the final output deterministic estimator.

Figure 2 :
Figure 2: Example of neural network estimator architecture utilizing an existing estimator.The arrows from the input nodes to the existing estimator are omitted from this graph.

Figure 3 :
Figure3: Estimated Bayes risks of the estimator over iterations when computing a Γ 1 -minimax estimator.The lines are the current Bayes risks (y-axis) over iterations (x-axis) (unbiased estimates with 50 Monte Carlo runs for Algorithm 6; exact values for Algorithm 3).The solid lines are the Bayes risks after an update in the estimator to decrease the Bayes risk.The dashed lines are the Bayes risks after an update in the prior to increase the Bayes risk.The two horizontal lines are the Bayes risk of the sample mean (dashed) and d * (solid), respectively, for π * .For ease of visualization, in subfigures (a) and (b), the Bayes risks are plotted every 50 iterations; in subfigures (c) and (d), the Bayes risks are plotted every 200 iterations; subfigure (d) contains the part in subfigure (c) after 500 iterations.

Figure 4 :
Figure4: Architecture of the neural network estimator of the expected number of new categories.X k : number of categories with k observations; OSW: the estimator proposed inOrlitsky et al. (2016); SCL: the estimator proposed inShen et al. (2003).The arrows from data (X 1 , . . ., X n ) to the OSW and SCL estimators are omitted from this graph.

Figure 8 :
Figure 8: Architecture of the deeper and wider neural network estimator of the expected number of new categories.

Table 1 :
Risks and Bayes risks of estimators.R(d, P 0 ): risk of the estimator under the true data-generating mechanism P 0 .r(d, π * ): Bayes risk under prior π * , the computed prior from Algorithm 6 in the last and finest grid in the computation.

Table 2 :
Predicted number of new categories (rounded to the nearest integer) in a new sample with size 2000 based on the sample with size 1088 studied in

Table 3 :
Risks and Bayes risks of estimators.R(d, P 0 ): risk of the estimator under the true data-generating mechanism P 0 .r(d, π * ): Bayes risk under prior π * , the computed prior from Algorithm 6 in the last and finest grid in the computation.

Table 4 :
Miller and Wiegert (1989)n the sample with size 1088 studied inMiller and Wiegert (1989).The strength of prior information in Γ-minimax estimators is shown in brackets.
ulation in the simulation.The tuning parameters of the Γ-minimax estimators are exactly the same as those in the above simulation except that the starting point in Algorithm 4 has more categories.The estimates are presented in Table4.All methods produce almost identical estimates.Because the sample size is more than ten times the sample size in the simulation and the JVHW estimator is minimax rate optimal

Table 5 :
Coefficients and Bayes risks of estimators of the mean.Unrestricted space: the space of all estimators that are square-integrable with respect to all P ∈ M.

Table 6 :
Table similar to Table 1 for sensitivity analysis with strongly informative prior information.linear .Even in this exceptional case, the computed Γ 1 -minimax estimator is still approximately Γ-minimax.We think the algorithm does not stop then in these cases because of Monte Carlo errors when computing r sup (d * ℓ−1 , Γ ℓ ) and r sup (d * ℓ−1 , Γ ℓ−1 D ′ , Π ′ ).Suppose that D ′ = MV t (D i )/t ≤ A, sup j,t MV t (Π j )/t ≤ A for A < ∞.If all D i and Π j are eligible in [s, s + t], then max D ′ U s+t − min D ′ U s+t ≤ (2B + A)t and max Π ′ V s+t − min Π ′ V s+t ≤ (2B + A)t.Proof of Lemma 4. Without loss of generality, assume that d ∈ (argmax d∈D ′ U s+t ) D 1 .Since D 1 is eligible in [s, t], there exists t ∈ [s, s + t] such that (argmin d∈D ′ U t) D 1 ̸ = ∅.By the recursive definition of the sequence {U t } t in (4), the bound on the risk, and the assumption that sup i,t MV t (D i )/t ≤ A, we have that maxD ′ U s+t = U s+t ( d) ≤ U t( d) + B(s + t − t) ≤ min D ′ U t + At + B(s + t − t) ≤ min D ′ U t + (A + B)t.Letting d′ ∈ argmin d∈D ′ U s+t , by the bound on the risk, we can derive thatmin D ′ U s+t = U s+t ( d′ ) ≥ U t( d′ ) − B(s + t − t) ≥ min D ′ U t − Bt.Combine these two inequalities and we have that max D ′ U s+t − min D ′ U s+t ≤ (2B + A)t.An identical argument applied to the sequence {V t } t

Table 7 :
Summary of frequently used symbols Symbol P 0 True data-generating mechanism M Space of data-generating mechanisms containing P 0 X * and X = C(X * ) Full generated data and coarsened data D Space of candidate estimators or decision functions (e.g., neural networks) An increasing sequence of finite subsets of M Γ ℓ Set of priors in Γ with support in M ℓ r sup Worst-case Bayes risk function r sup: (d, Γ ′ ) → sup π∈Γ ′ r(d, π) d * ℓ Γ ℓ -minimax estimator in D d * A limit point of sequence {d * ℓ } ∞ ℓ=1 , which is Γ-minimax in D by Theorem 1 β(∈ R D )Coefficient of a finite-dimensional estimator (e.g., neural network) * = d(ϖ * ) Γ-minimax estimator in D