The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
98,475 characters
Empirical Bayes for compound adaptive experiments
\pagenumbering{gobble}
\title{Empirical Bayes for compound adaptive experiments}
\author{Karun Adusumilli\textsuperscript{*}}
\author{Jiaying Gu\textsuperscript{\textdagger}}
\author{Junfan Tao$^\ddagger$}
\begin{abstract}
We investigate Empirical Bayes (EB) methods in the context of compound adaptive experiments, where the arm distribution in each experiment follows a normal distribution with an unknown mean that we seek to estimate. There are two main EB strategies: $g$-modeling, which estimates the prior by maximizing the marginal likelihood, and $f$-modeling, which derives posterior means directly from the empirical distribution of the observations. We show that $g$-modeling continues to be a valid EB procedure even when it incorrectly assumes that data are collected exogenously; its validity does not depend on the particular sampling algorithm or on whether sample sizes are endogenous. In practice, one can apply standard
$g$-modeling techniques by acting as though the data were exogenously sampled. We extend regret guarantees from exogenous sampling to adaptively generated data. By contrast, naively applying the Tweedie formula based on the marginal density of the observed data, as in standard $f$-modeling, can produce biased rules under adaptive sampling. We corroborate the robustness of $g$-modeling through simulations with widely used adaptive algorithms and demonstrate its applicability using a real-world dataset consisting of multiple sequential experiments.\end{abstract}
\thanks{\textit{This version}: \today{}. \\
\textsuperscript{*}Department of Economics, University
of Pennsylvania; [email removed].\\
\textsuperscript{\textdagger}Department of Economics, University of Toronto;
[email removed].\\
$^{\ddagger}$Institute of Economic Research, Kyoto University;
[email removed].
}
\maketitle
\newpage
\pagenumbering{arabic}
\section{Introduction \protect\label{sec:Introduction}}
Recent years have seen remarkable advancements in the theory and application of adaptive experiments. These cutting-edge designs are now widely adopted across a range of disciplines, including online advertising \citep{russo2017tutorial}, dynamic pricing \citep{ferreira2018online}, drug discovery \citep{wassmer2016group}, public health \citep{athey2021shared}, and economic interventions \citep{kasy2019adaptive}. Compared to traditional randomized trials, adaptive experiments offer a more efficient and flexible framework for balancing welfare, ethical, and economic considerations. Their potential has been widely recognized---e.g., since the launch of the Critical Path Initiative in 2006, the FDA has actively promoted the adoption of adaptive designs in clinical trials to reduce costs and mitigate risks for participants.
While a rich literature has introduced innovative algorithms for implementing adaptive designs, fundamental challenges remain in post-experiment estimation. For example, technology firms routinely run thousands of adaptive experiments, commonly known as Online Controlled Experiments (OCEs), to assess and improve digital products and services. In a standard OCE, visitors to a website are randomly assigned either to a control group or to one of $K$ treatment groups, and the allocation probabilities are continually adjusted in response to interim outcomes. Estimating treatment effects in such adaptive frameworks, however, poses distinctive challenges. Conventional methods, like simple differences in sample means, break down because the adaptive data collection process violates the assumptions underlying traditional statistical procedures. In practice, Bayesian methods are therefore often used, where a prior is placed on the treatment effect, and the posterior mean is reported. This strategy, however, makes the results sensitive to the choice of prior and can yield biased estimates when the prior is misspecified.
These challenges raise important methodological questions: How can treatment effects be estimated efficiently and reliably in adaptive experimental settings? Can information across multiple experiments be aggregated to enhance estimation accuracy? Moreover, can we learn the prior distribution of the treatment effects directly from the data?
To address these questions, this paper investigates the use of Empirical Bayes (EB) methods in the context of compound adaptive experiments, where multiple adaptive experiments share a common but unknown distribution of effect sizes. As noted by \cite{efron2014two}, there are two main strategies for EB estimation: $g$-modeling, which estimates the prior by maximizing the marginal likelihood, and $f$-modeling, which computes the posterior means directly from the sample distribution of the observations. Our key contribution is to establish that $g$-modeling remains valid even under adaptive sampling, whereas $f$-modeling fails.
Remarkably, we show that \( g \)-modeling requires no knowledge of the data-generating algorithm. Standard \( g \)-modeling procedures can be applied as if the data were exogenously drawn and the sample means were normally distributed - even when neither assumption holds in reality. Moreover, each treatment arm can be analyzed independently, even when some or all arms are part of a common adaptive design. Despite these apparent misspecifications, we establish that \( g \)-modeling retains similar risk guarantees in adaptive settings as it does under independently and identically distributed (i.i.d.) data.
These two key properties---algorithm independence and the ability to treat each arm as an independent experiment---substantially extend the applicability of EB methods. For instance, \( g \)-modeling can estimate mean treatment effects within a single adaptive experiment, even when the design includes a large number of treatment arms. Such settings are common in practice; e.g., \citet{chapelle2011empirical} describe an application in display advertising where bandit algorithms allocate traffic across 5,910 possible ads, each functioning as a distinct treatment.
Of course, algorithm independence is not unique to $g$-modeling. As established by the likelihood principle (see Chapter 7 in \citealp{berger2013statistical}), any Bayesian estimation strategy also remains valid regardless of the sampling algorithm. Since $g$-modeling is simply Bayesian estimation with an estimated prior, its posterior estimates naturally inherit this robustness. The novel insight of this paper, however, is that the prior itself can be estimated in an algorithm-independent manner.
To explain why $g$-modeling consistently estimates the prior despite the misspecification of the likelihood, we introduce a novel interpretation of $g$-modeling as a moment-matching procedure. Specifically, we show that $g$-modeling aligns the sample moments of the posterior distribution with those of the estimated prior. At the population level, the prior is identified through moments implied by the law of iterated expectations, which requires the prior moments to equal the expected posterior moments. The moment-matching perspective suggests that $g$-modeling effectively solves for the prior by leveraging the sample analogs of these identifying conditions. Consequently, for parametric prior classes, the validity of $g$-modeling can be established using standard generalized method of moments arguments.
When the class of candidate priors is left unspecified, as in non-parametric maximum likelihood (NPMLE) procedure, we show that $g$-modeling estimates a probability distribution satisfying a self-consistency property: it equals the average posterior distribution corresponding to itself. This is again a sample version of the population requirement that the prior be equal to the expected posterior.
We also formally demonstrate the regret consistency of the NPMLE procedure. This relies on a novel extension of the seminal results of \cite{jiang2009general} and \cite{jiang2020general}, and may be of independent of interest.
As is standard in the EB literature, we derive our results on regret consistency assuming Gaussian outcomes. However, we argue that even with general parametric models the Gaussian likelihood naturally emerges as an approximation to the true likelihood within a local asymptotic regime. Effectively, the $g$-modeling framework continues to be valid if we use the sample means of the score functions in place of the sample means of the outcomes.
Beyond treatment effect estimation, knowledge of the prior can improve future decision making and algorithmic design. For instance, \cite{azevedo2020b} utilized EB methods to estimate the prior across multiple static Online Controlled Experiments (OCEs). Their findings revealed that the distribution of treatment effects is typically fat-tailed, leading them to propose a ``lean'' experimentation strategy--conducting a larger number of experiments, each with smaller sample sizes. Our methods allow us to extend this analysis to adaptive experiments.
As an illustration of our methods, we apply the $g$-modeling strategy to the ASOS digital experiments dataset \cite{liu2021datasets}, which comprises 61 adaptive experiments run by a business unit within ASOS, a fashion retail platform. The specific algorithms used to generate this data are proprietary and not publicly disclosed, making this an ideal setting to evaluate the robustness and applicability of our approach. Our results appear to suggest a thick tailed distribution of treatment effects, though some caution is warranted given the small sample of experiments in our analysis.
\@startsection{subsection}{2}{\z@}
{0pt}
{0.4ex}
{\normalfont\normalsize\bfseries}{Related literature}
There is a substantial and growing literature on the design of adaptive experiments; see \cite{lattimore2020bandit} and \cite{wassmer2016group} for comprehensive surveys. In contrast, research on the estimation of treatment effects following adaptive experiments remains relatively limited. While the standard sample mean estimator is consistent at parametric rates, it is typically asymptotically biased and fails to be normally distributed. One approach to addressing this issue is inverse propensity weighting, as proposed by \cite{hadad2021confidence}, which restores asymptotic normality. However, this method requires knowledge of the data-generating algorithm and applies only to specific classes of non-deterministic algorithms, excluding widely used methods such as Upper Confidence Bound (UCB) and Bayes-optimal algorithms.
Alternative approaches, such as those in \cite{deshpande2018accurate} and \cite{nie2018adaptively}, attempt to de-bias the sample mean estimator in adaptive experiments, producing asymptotically normal estimates but often at the cost of increased variance. Meanwhile, a large applied literature eschews frequentist methods in favor of Bayesian estimation. While Bayesian estimators are Bayes-optimal by construction, their performance is highly sensitive to the choice of prior, raising concerns about robustness in practical applications.
We contribute to this literature by developing methods to estimate the prior from the data and extending the Empirical Bayes framework to multiple adaptive experiments, thereby improving both the efficiency and robustness of estimation.
While the foundations of EB trace back to \cite{robbins1951asymptotically}, the methodology has recently experienced a resurgence, particularly in economics, driven by applications in labor economics and related fields. For comprehensive overviews, see \cite{efron2012large}, \cite{walters2024empirical}, \cite{koenker2024empirical} and \cite{koenker_gu_2026}. Beyond treatment effect estimation, the EB framework has broad applicability to various compound decision problems, including multiple testing \citep{efron2012large} and ranking \citep{gu2023invidious}.
\cite{efron2014two} describes $g$-modeling and $f$-modeling approaches to EB estimation. This paper provides a novel interpretation of $g$-modeling as a moment-matching procedure. This perspective is based on \cite{adusumilli2020unobserved}, itself an extension of the seminal work of \cite{neal1998view} on the variational interpretation of the EM algorithm.
The non-parametric maximum likelihood estimator (NPMLE), first introduced by \cite{kiefer1956consistency}, is a flexible approach to \( g \)-modeling that estimates the prior without imposing parametric assumptions on its form. In recent years, the NPMLE has gained traction, driven by computationally efficient algorithms proposed by \cite{koenker2014convex} and strengthened by new theoretical analyses from \cite{jiang2009general}, \cite{jiang2020general}, and \cite{polyanskiy2020self}. \cite{gilraine2020new} apply this method to non-parametrically estimate distributions of teacher value-added. In this paper, we develop new theoretical results demonstrating that the regret consistency of the NPMLE carries over to adaptive experimental settings.
Our results indicate that standard \( g \)-modeling approaches, which assume a Gaussian likelihood for the sample means, remain valid even for adaptively sampled data, despite the fact that the sample mean may no longer be Gaussian or even a sufficient statistic. Because sample sizes are random, this mis-specified likelihood naturally introduces heteroskedasticity across experiments. Recently, \cite{chen2022empirical} proposed the CLOSE framework for analyzing such heteroskedastic models. However, our findings indicate that CLOSE loses its validity under our setup of compound adaptive experiments. In our setup, sample sizes are independent of the unknown parameters when conditioned on observed data . The breakdown of CLOSE occurs for two reasons. First, CLOSE assumes a specific relationship between parameters and sample sizes, which may not hold under adaptivity. Second, and more critically, CLOSE treats sample sizes as parameters, whereas in adaptive experiments, they are functions of the data. As we demonstrate in Section \ref{subsec:CLOSE}, this misalignment causes CLOSE to apply Bayesian updating twice to the same data, leading to estimates that are insufficiently shrunk.
\section{Compound Adaptive Experiments\protect\label{sec:Setup}}
\@startsection{subsection}{2}{\z@}
{0pt}
{0.4ex}
{\normalfont\normalsize\bfseries}{Motivating example: The ASOS Digital Experiments dataset} \protect\label{subsec:Motivating example}
Online Controlled Experiments (OCEs) are web-based randomized controlled trials designed to evaluate and improve digital products and services. These experiments often incorporate adaptive stopping rules, enabling data collection to stop when pre-specified criteria are met. Technology companies routinely conduct thousands of OCEs daily to refine user experiences and optimize product performance. In a typical OCE, website users are randomly assigned to either a control group or one of \( K \) treatment groups. OCEs with \( K = 1 \) are commonly referred to as A/B tests.
Between 2019 and 2020, the global fashion retailer ASOS.com conducted
$78$ such OCEs, documented in the ASOS Digital Experiments Dataset \citep{liu2021datasets}. This dataset includes daily sample means for the treatment and control groups, alongside the number of samples per treatment arm. The retailer used an adaptive stopping criterion to determine the end of each experiment, but the algorithms used are proprietary and unknown (to us).
For simplicity, we assume each ASOS experiment involves a single treatment arm (i.e., it is an A/B test).\footnote{In reality, 17 experiments included multiple treatments; we drop these experiments, resulting in $n = 61$ experiments in total. See Section \ref{sec:Empirical illustration} for more details.} In each experiment, let $j$ index the periods of experimentation, and in each period a single treatment and control observation is drawn with equal proportion, but we only observe the difference $Y_{j,i} = Y_{j,i}^{(1)} - Y_{j,i}^{(0)}$ between the treated and control observation. The rationale for this assumption is discussed in Section \ref{subsec: Generalizing}. We assume that the outcome differences follow a normal distribution, \( Y_{j,i} \sim \mathcal{N}(\tilde{\theta}_i, \omega_i^2) \), where \( \omega_i^2 \) is known. The Gaussianity assumption is relaxed in Section \ref{subsec:local asymptotics}.
Let $A_{j,i} \in \{0,1\}$ indicate whether experiment $i$ stops in period $j$, with $A_{j,i} = 1$ denoting a stopping decision. The stopping rule is governed by a policy $\pi_{j,i}$, which maps past information to a probability of stopping. Formally, stopping occurs when $ A_{j,i} = \mathbb{I}\{U_{j,i} < \pi_{j,i}\}$, where $U_{j,i} \sim \textrm{Uniform}[0,1]$ is an exogenous random variable (i.e., independent of all past data) encoding policy randomization.
Define the information set at period $j$ as $\mathcal{I}_{j,i} \equiv \sigma\{ Y_{1,i}, U_{1,i}, \dots, Y_{j-1,i}, U_{j-1,i} \}.$ The stopping policy is thus a function
\[
\pi_{j,i}: \mathcal{I}_{j,i} \to [0,1].
\]
The policies are implicitly restricted to depend only on past outcome differences and not on the outcome levels themselves. As discussed in Section \ref{subsec: Generalizing}, this restriction holds for most stopping rules commonly used in A/B testing.
The stopping time for experiment \( i \) is denoted by \( N_i \), and the realized sample from each experiment is given by \( \mathcal{D}_i \equiv \mathcal{I}_{N_i + 1,i} \). The sample mean of the treatment effects is defined as
\[
\tilde{Z}_i = \frac{1}{N_i} \sum_{j=1}^{N_i} Y_{j,i}.
\]
In the ASOS dataset, the observed sample means are small (on the order of \( 10^{-3} \)), while the realized \( N_i \) values are large (on the order of \( 10^5 \)). Such relative magnitudes are common in online experimentation. \cite{deng2013improving} survey A/B testing practices in digital environments and report that treatment effects often amount to less than 1\% of the expected control outcome. Despite their small size, these effects can generate substantial revenue gains when implemented at scale, motivating companies to conduct large experiments to detect them with statistical precision.
To facilitate accurate estimation of these small treatment effects, we rescale them as $\theta_i = \sqrt{N} \tilde{\theta}_i$, where
\begin{equation} \label{eq:definition of N}
N := \frac{1}{n}\sum_{i=1}^n \mathbb{E}[N_i].
\end{equation}
Correspondingly, we also rescale $N_i$ and $\tilde{Z_i}$ as
\begin{align} \protect \label{Definition of t_i, Z_i}
\tau_i & := \delta_\tau(\mathcal{D} _i) = N_i/N, \textrm{ and} \nonumber \\
Z_i &:= \delta_z(\mathcal{D}_i) = \sqrt{N} \tilde{Z}_i.
\end{align}
This ensures that the orders of magnitude of $\tau_i, Z_i$ are stable even as $N$ increases.\footnote{The scaling factor $N$ is used only as a conceptual tool to justify focusing on small treatment effects. Our methods do not depend on knowing its value.}
To motivate our Empirical Bayes (EB) procedures, we start by discussing the estimation of $\theta_i$ in a Bayesian setting where each $\theta_i$ is an independent random draw from a known prior $G_0$.
\subsubsection*{The likelihood principle and Bayes estimation}
We can write the likelihood of the observed history $\mathcal{D} _i$ as
\begin{align} \protect \label{eq:derivation of likelihood}
&p(\mathcal{D} _i | \theta_i) = \prod_{j=1}^{N_i} p\Large(A_{j,i}, U_{j,i}, Y_{j,i}| \mathcal{I}_{j,i}, \theta_i \Large) \nonumber \\
& = \left[\prod_{j=1}^{N_i-1} p(A_{j,i} = 0, U_{j,i}|\mathcal{I}_{j,i}, \theta_i)\right]
p(A_{N_i,i} = 1, U_{N_i,i}|\mathcal{I}_{N_i,i}, \theta_i)\cdot \prod_{j=1}^{N_i} p(Y_{j,i}|A_{j-1,i} =0, U_{j,i}, \mathcal{I}_{j,i}, \theta_i) \nonumber\\
& = \left[ \prod_{j=1}^{N_i-1} 1\{U_{j,i}\geq \pi_{j,i} \} \right]1\{U_{N_i,i}< \pi_{N_i,i} \}\cdot
\prod_{j=1}^{N_i} p(Y_{j,i}|\theta_i) ,
\end{align}
where the last equality follows from the facts: (1) $A_{j,i}$ is determined by the policy $\pi_{j,i}$ as a function solely of the past history and exogenous randomization $U_{j,i}$ (the actions and the policy can not depend on $\theta_i$ as it is unknown); (2) $U_{j,i}$ is uniformly distributed, implying its density is 1; and (3) the distribution of $Y_{j,i}$ given $\theta_i$ is independent of the past history of observations because, conditional on continuing the experiment, the outcomes are just a random draw from $\mathcal{N} (\theta_i/\sqrt{N}, \omega_i^2)$.
Making use of the normality of outcomes, we can represent $p(\mathcal{D} _i \vert \theta_i)$ as
\begin{equation}\label{likelihoodPrinciple}
p(\mathcal{D} _i |\theta_i)= \pi_i(\mathcal{D} _i) \cdot \frac{1}{\sigma_i} \varphi \big(\frac{Z_i - \theta_i}{\sigma_i} \big)
\end{equation}
where $\varphi(\cdot)$ represents the standard normal density,
$$
\sigma_i^2 := \delta_\sigma(\mathcal{D} _i) = \frac{\omega_i^2}{\tau_i},
$$ and
$$
\pi_i(\mathcal{D} _i) := p(\mathcal{D} _i \vert 0) \cdot \sqrt{2\pi\sigma_i^2} \exp \left\{\frac{Z_i^2}{2\sigma_i^2} \right\}.
$$
Note that $\pi_i(\mathcal{D} _i)$ is independent of $\theta_i$.
Suppose instead that the samples were drawn exogenously. In a such a scenario, we would have
$$
Z_i \mid \theta_i \sim \mathcal{N}(\theta_i, \sigma_i^2),
$$
leading to the likelihood
\[
\frac{1}{\sigma_i} \varphi \bigg(\frac{Z_i - \theta_i}{\sigma_i} \bigg)
\]
for the observed data. We refer to the above as the `working likelihood' for our Empirical Bayes procedures. Equation (\ref{likelihoodPrinciple}) implies that the working likelihood differs from the true likelihood $p(\mathcal{D} _i\vert \theta_i)$ by a multiplicative constant $\pi_i(\mathcal{D} _i)$ that is independent of $\theta_i$.
Consider the Bayesian approach that assigns a prior \( G_0 \) to \( \theta_i \). Since the true and working likelihoods differ only by a multiplicative factor free from $\theta_i$, using either leads to the same posterior distribution over \( \theta_i \). This is because of the likelihood principle: Bayesian updating depends solely on the observed data and does not require detailed knowledge of the data-generating algorithm.
This independence from the underlying sampling process extends naturally to estimation. When the goal is to estimate the scaled treatment effect \( \theta_i \) under squared error loss, the optimal Bayes estimator is the posterior expectation, \( \mathbb{E}_{G_0}[\theta_i | \mathcal{D}_i] \), where \( \mathbb{E}_{G_0}[\cdot| \mathcal{D}_i] \) denotes the expectation conditioned on the observed data \( \mathcal{D}_i \) and the prior \( G_0 \). Using (\ref{likelihoodPrinciple}), the posterior expectation of \( \theta_i \) conditional on \( \mathcal{D}_i \) can expressed as:
\[
\mathbb{E}_{G_0}[\theta_i \mid \mathcal{D}_i] = \frac{\int \theta_i \frac{1}{\sigma_i} \varphi \left(\frac{Z_i - \theta_i}{\sigma_i} \right) \, dG_0}{\int \frac{1}{\sigma_i} \varphi \left(\frac{Z_i - \theta_i}{\sigma_i} \right) \, dG_0}.
\]
We observe that the posterior mean is identical to the one that would be obtained if the data were assumed to have been generated exogenously, irrespective of the actual adaptive nature of the sampling process.
The Bayesian approach to estimation is particularly appealing because it is independent of the data collection algorithm. However, it is highly sensitive to the choice of a subjective prior. To address this limitation, our EB methods, introduced in Section \ref{sec:Empirical Bayes}, offer a data-driven and robust alternative by estimating the prior directly from observed data. We also demonstrate that the $g$-modeling-based EB procedure retains a key advantage of Bayesian estimation: its independence from the data generation algorithm.
Before delving into the EB methodology, we first generalize the problem setup. This is followed by a detailed discussion of the crucial common prior assumption for $(\theta_1, \dots, \theta_n)$.
\@startsection{subsection}{2}{\z@}
{0pt}
{0.4ex}
{\normalfont\normalsize\bfseries}{Generalizing the setup: Compound adaptive experiments} \protect\label{subsec: Generalizing}
Consider a setting in which an analyst observes \( n \) experiments, each associated with an unknown parameter \( \theta_1, \dots, \theta_n \). Let \( \mathcal{D}_i \) denote the data generated from experiment \( i \) and \( \mathcal{D} \) the combined data across all experiments. Our methodology applies to any setting where the likelihood takes the form
\begin{equation}\label{likelihoodPrinciple:generalization}
p(\mathcal{D} \mid \theta_1,\dots,\theta_n) = c(\mathcal{D}) \cdot \prod_{i=1}^n \frac{1}{\sigma_i} \varphi \left(\frac{Z_i - \theta_i}{\sigma_i} \right),
\end{equation}
where \( Z_i \) and \( \sigma_i \) are scalar statistics that depend only on \( \mathcal{D}_i \).
Equation \eqref{likelihoodPrinciple:generalization} implies that the likelihood is proportional to the usual working likelihood for Empirical Bayes methods, where the data from each experiment \( i \) follows a normal model,
\begin{equation} \label{eq:working_likelihood}
Z_i \mid \theta_i \sim \mathcal{N}(\theta_i, \sigma_i^2),
\end{equation}
with observations across experiments being mutually independent. The formulation \eqref{likelihoodPrinciple:generalization} turns out to be broad enough to encompass nearly all known classes of compound adaptive experiments. Below, we discuss specific examples that illustrate its applicability.
\begin{example}[Multi-arm adaptive experiments]\label{example 1}
Consider a setting in which an analyst observes one or more adaptive experiments, each involving multiple treatment arms. Our methodology treats each treatment arm as an independent experiment, regardless of whether some or all arms belong to the same adaptive experiment. Let \( i \) index the treatment arms, and model the outcomes from arm \( i \) as $Y_{j,i} \sim \mathcal{N} \left(\tilde{\theta}_i, \omega_i^2\right)$,
where \( \tilde{\theta}_i \) represents the unknown mean outcome for arm \( i \), and \( \omega_i^2 \) is assumed to be known.
Let \( N_i \) denote the number of times arm \( i \) has been sampled, and define the observed data for arm \( i \) as \( \mathcal{D}_i := \{ Y_{1,i}, \dots, Y_{N_i,i} \} \), the set of all outcomes generated from that arm. Define \( \tilde{Z}_i \) as the sample mean of outcomes from arm \( i \).
Our framework accommodates a broad range of adaptive sampling algorithms, allowing them to vary across arms. Moreover, the sampling algorithm for a given arm \( i \) may depend on data from other arms or experiments. For instance, if Arms 1 and 2 belong to the same adaptive experiment, data from Arm 1 may influence how often Arm 2 is sampled. Alternatively, if the arms correspond to separate experiments, data from an earlier experiment may inform the algorithm used in a subsequent experiment. Our setup allows for both scenarios. The only substantive restriction is that the choice of sampling algorithms cannot directly depend on the true treatment effects \( \tilde{\theta}_1, \dots, \tilde{\theta}_n \).
Define $N$ as in (\ref{eq:definition of N}). In practice, \( N \) is typically large because the mean effects \( \tilde{\theta}_i \) are often close to each other, requiring large samples to reliably distinguish between treatment arms. As before, we employ a local-to-zero reparameterization, $\theta_i := \sqrt{N} \tilde{\theta}_i$, to account for such small difference, and also rescale \( N_i \) and \( \tilde{Z}_i \) to obtain the normalized variables \( \tau_i \) and \( Z_i \) as in \eqref{Definition of t_i, Z_i}. Our primary objective is to estimate the scaled treatment effects \( \theta_1, \dots, \theta_n \), postulated to be random draws from a prior \( G_0 \).
By applying arguments similar to those in \eqref{eq:derivation of likelihood}, we establish that the overall conditional data density follows the structure of \eqref{likelihoodPrinciple:generalization}.
\begin{prop} \protect \label{Prop:likelihood principle}
In the context of compound adaptive experimentation with multiple treatment arms indexed by \( i \), suppose the sampling algorithm for any given arm \( i \) is conditionally independent of \( (\theta_1, \dots, \theta_n) \), given the data from all other arms and an exogenous randomization. Then, the overall conditional data density \( p(\mathcal{D} \mid \theta_1, \dots, \theta_n) \) takes the form
\begin{equation}\label{likelihoodPrinciple:MAB}
p(\mathcal{D} \mid \theta_1, \dots, \theta_n) = c(\mathcal{D} ) \cdot \prod_{i=1}^n \frac{1}{\sigma_i} \varphi \left(\frac{Z_i - \theta_i}{\sigma_i} \right),
\end{equation}
for some function \( c(\mathcal{D} ) \) that does not depend on $\theta_1, \dots, \theta_n$.
\end{prop}
See Appendix \ref{sec:Appendix:A} for the formal proof.
\end{example}
\begin{example}[Panel data with missingness and attrition]
Consider a panel data setting where outcomes are modeled as \( Y_{j,i} = \theta_i + \epsilon_{j,i} \), with \(\epsilon_{j,i} \sim \mathcal{N}(0, \omega_i^2)\). Here, the goal is to estimate the parameters \(\theta_i\). Observations may be subject to missingness or attrition, assumed to occur at random given the history of past observations. The probabilities of missingness and attrition are unspecified and can vary across both time periods \(j\) and individuals \(i\). This scenario aligns with the structure of an adaptive multi-arm experiment, implying that Proposition \ref{Prop:likelihood principle} applies to this setting as well.
\end{example}
\begin{example}[Correcting for p-hacking]
P-hacking refers to a range of practices by which researchers selectively report statistically significant results. Here, we focus on specific forms of p-hacking that arise when the experimental protocol is mis-characterized or incompletely reported---even as all relevant data are disclosed.
Examples include optional stopping, where data collection continues until statistical significance is achieved, and selective subgroup analysis, where only significant results from specific subgroups are reported. These forms of p-hacking are equivalent to adaptive experiments with an unknown sampling algorithm, placing them squarely within the scope of our framework.\footnote{Other forms of p-hacking include selectively replacing or failing to disclose relevant data. These possibilities fall outside the remit of our methods.}
Bayesian methods are robust to selective reporting of this kind, but are sensitive to the choice of prior. Our Empirical Bayes approach enables estimation of the prior from a meta-analysis of potentially p-hacked studies. The resulting Bayes estimators inherently adjust for the distortions introduced by p-hacking.
\end{example}
\begin{example}[Multiple A/B tests]
The multiple A/B testing problem, introduced in Section \ref{subsec:Motivating example}, is a special case of multi-arm adaptive experiments in which each experiment consists of exactly two arms: a treatment and a control. While these experiments can therefore be analyzed similarly to Example \ref{example 1}, the assumption that mean effects across arms are drawn from the same prior may be less credible. This is because the control arm typically represents the status quo and is not necessarily exchangeable with the treatment arm. Moreover, in some cases, the same control arm may be used across multiple experiments, further complicating standard exchangeability assumptions.
A complete analysis of multiple A/B testing would require specifying a joint prior on the mean effects of both the treatment and control arms. However, multivariate EB procedures are more complex, and for tractability, we adopt a simpler approach, as in Section \ref{subsec:Motivating example}, based on the following assumptions:
\begin{enumerate}
\item The arms are sampled in equal (or known) proportions.
\item Only the outcome differences between the treatment and control arms are observed in each period.
\item The stopping time depends only on past treatment effect differences.
\end{enumerate}
The first assumption is standard for A/B tests and is known to hold in the ASOS example. The second assumption implies that we set aside information about the mean outcome in each individual arm and instead concentrate solely on the differences between arms. This strategy is standard in EB analysis, even though it entails some loss of efficiency because only part of the available data is used for estimation. It can be supported by an equivariance requirement stating that the decision rule should remain unchanged if a constant is added to all outcomes. The third assumption is unverifiable in the ASOS context, but in practice, almost all A/B testing stopping rules are based solely on the difference in sample means. This includes Wald's seminal sequential probability ratio test (SPRT; \citealt{wald1945statistical}), various group-sequential methods from the clinical trials literature \citep{wassmer2016group}, and the minimax-optimal stopping rule of \cite{adusumilli2022sample}. Moreover, it can be shown that stopping times based only on the difference in sample means form an asymptotically complete class.
Under these assumptions we can simplify the analysis by focusing on the difference in means and placing a prior directly on the treatment effects, as we did in Section \ref{subsec:Motivating example}. It then follows straightforwardly that an analog of Proposition \ref{Prop:likelihood principle} holds, where \( Z_i \) now represents the scaled difference in sample means at the conclusion of experiment \( i \).
\end{example}
\@startsection{subsection}{2}{\z@}
{0pt}
{0.4ex}
{\normalfont\normalsize\bfseries}{On the common prior assumption}
Throughout this paper, we adopt the standard Empirical Bayes (EB) assumption that \(\theta_1, \dots, \theta_n\) are i.i.d draws from a common prior. This assumption implies that any specific realization \((\theta_1, \dots, \theta_n)\) is just as likely as any of its permutations \((\theta_{\kappa(1)}, \dots, \theta_{\kappa(n)})\), ensuring exchangeability of the parameter vectors. As \cite{efron2012large} notes, this is closely tied to the question of `comparability': ``\textit{Empirical Bayes methods involve each case learning from the experience of others... To make this believable, the `others' have to be similar in nature to the case at hand, or at least not obviously dissimilar}''.
The common prior seems particularly natural in the context of analyzing multiple arms within the same adaptive experiment: after all, there is no a priori reason to expect that the ordering of the arms carries any intrinsic significance. However, when considering multiple adaptive experiments, the validity of this assumption depends on the context. For instance, it would be highly questionable to pool adaptive experiments on ad-targeting with those on drug discovery. In the context of the ASOS dataset, the various experiments were run by a single business unit within the company. Moreover, it is reasonable to assume that the algorithms used were independent of the parameters, given the observed data (as required for Proposition \ref{Prop:likelihood principle}). Therefore, in this setting, we consider the common prior assumption to be much more justifiable.
The term \textit{compound decision problem} refers to a setting in which an experimenter seeks to estimate \(\theta_1, \dots, \theta_n\) by minimizing an average frequentist mean-squared error criterion (i.e., averaged across the $n$ experiments, see Section \ref{sec:Theoretical results}). The compound decision problem is not associated with a prior and is generally different from the EB problem. However, the fundamental theorem of compound decisions \citep{robbins1951asymptotically} establishes that these two problems are equivalent if the experiments are exchangeable---meaning the conditional likelihood for experiment \(i\) has the same form across all experiments---and if the prior \(G_0\) is taken to be the empirical distribution \(n^{-1} \sum_i \delta_{\theta_i}\).
In our setting, the sampling algorithms are allowed to differ across experiments, so exchangeability of experiments may not hold.\footnote{Note, however, that this is separate from the exchangeability of parameters, as required for Empirical Bayes.} However, if all arms belong to the same adaptive experiment or if identical algorithms were used across experiments with constant outcome variances (\(\omega_i^2 = \omega^2\) for all \(i\)), then exchangeability of experiments is preserved too, and our methods also solve the compound decision problem.
\section{Empirical Bayes Methodology for Estimating the Prior\protect\label{sec:Empirical Bayes}}
Continuing with the general setup of Section \ref{subsec: Generalizing}, let $G_0$ denote the common prior for $\theta_1,\dots ,\theta_n$. For any given candidate prior \( G \), the marginal density of the data \( \mathcal{D} \) is given by
$$
p_{G}(\mathcal{D} ) = \int p(\mathcal{D} |\theta_1, \dots, \theta_n)dG^{(n)}(\theta _1,\dots,\theta_n),
$$
where $G^{(n)}(\cdot)$ denotes the product prior over $\theta_1,\dots ,\theta_n$ corresponding to the marginal $G$. Using (\ref{likelihoodPrinciple:generalization}), $p_G(\mathcal{D} )$ can be written as:
\begin{equation} \label{eq:marginal_density}
p_{G}(\mathcal{D}) = c(\mathcal{D} ) \prod_{i=1}^n \int \frac{1}{\sigma_i} \varphi \bigg(\frac{Z_i - \theta_i}{\sigma_i} \bigg) \, dG(\theta_i ) = c(\mathcal{D} ) \prod_{i=1}^n f_{G, \sigma_i}(Z_i),
\end{equation}
where
\begin{equation} \protect \label{eq:definition_of_f_G}
f_{G, \sigma_i}(Z_i) := \frac{1}{\sigma_i} \int \varphi \bigg(\frac{Z_i - \theta_i}{\sigma_i} \bigg) \, dG(\theta_i).
\end{equation}
The term \( f_{G, \sigma_i}(Z_i) \) represents the marginal density of \( Z_i \) derived from the working likelihood.
The core idea of the $g$-modeling approach is to estimate the true prior \( G_0 \) by maximizing the marginal likelihood \( p_{G}(\mathcal{D} ) \) over a candidate family of priors \( \mathcal{G} \). Different choices of \( \mathcal{G} \) define different $g$-modeling procedures. For example, \( \mathcal{G} \) could represent the class of all Gaussian priors with mean 0, resulting in a parametric $g$-modeling approach. Alternatively, \( \mathcal{G} \) could remain unrestricted, leading to the non-parametric maximum likelihood estimation (NPMLE) approach. In all cases, the $g$-modeling framework estimates \( G_0 \) as:
\[
\hat{G} = \operatorname*{arg\,max}_{G\in \mathcal{G}}\frac{1}{n} \ln p_{G}(\mathcal{D} ) = \operatorname*{arg\,max}_{G\in \mathcal{G}} \frac{1}{n} \sum_{i} \ln f_{G,\sigma_i}(Z_i).
\]
Importantly, \( f_{G,\sigma_i}(Z_i) \) is independent of the structure of the adaptive experiment. This independence implies that, from a computational perspective, $g$-modeling does not require any knowledge of how the dataset was collected.
Our key insight is that this independence also supports statistical guarantees for $g$-modeling. Under the assumptions stated in Section \ref{sec:Theoretical results}, maximizing the Gaussian working likelihood yields a regret-consistent plug-in Bayes rule even when using the `wrong' marginal likelihood \( \prod_i f_{G,\sigma_i}(Z_i) \), which presumes the data were generated exogenously.
Why does $g$-modeling remain valid under this mis-specification? As shown below, this robustness stems from a novel interpretation of $g$-modeling as a moment-matching procedure. Specifically, we show that $g$-modeling aligns the moments of the posterior with those of the estimated prior. Because the posterior is algorithm independent, this ensures validity of g-modeling in adaptive settings even under a mis-specified likelihood. The moment matching interpretation builds on the work of \cite{neal1998view}, as further developed in \cite{adusumilli2020unobserved}.
\@startsection{subsection}{2}{\z@}
{0pt}
{0.4ex}
{\normalfont\normalsize\bfseries}{$g$-modeling as moment matching\protect\label{sec:g modeling moment matching}}
The Donsker-Varadhan variational formula states that
\begin{equation}\label{DV_formula}
\ln f_{G,\sigma_i}(Z_i) = \max_{q_i(\cdot)} \bigg\{ \mathbb{E}_{q_i(\cdot)}\left[\ln \bigg\{ \frac{1}{\sigma_i} \varphi \big(\frac{Z_i - \theta_i}{\sigma_i} \big) \bigg\} \right] - \textrm{KL} \big(q_i(\cdot) \parallel G \big) \bigg\},
\end{equation}
where $q_i(\cdot)$ denotes some arbitrary probability distribution and $\textrm{KL}(P \parallel Q)$ denotes the Kullback-Leibler (KL) divergence between two probability measures $P$ and $Q$. The optimal value of $q_i(\cdot)$ in (\ref{DV_formula}) is just the posterior distribution $q^*_{i,G}(\cdot)$ corresponding to the prior $G$ and the working likelihood $
\varphi \big(\frac{Z_i - \theta_i}{\sigma_i} \big)/\sigma_i.
$ Based on the variational formula, we can rewrite the $g$-modeling optimization problem as
\begin{align}\label{EM_interpretation}
& \max_{G \in \mathcal{G}} \frac{1}{n}\sum_i \ln f_{G,\sigma_i}(Z_i) =\max_{G \in \mathcal{G}} \max_{\{q_i(\cdot)\}_i } \frac{1}{n} \sum_i\bigg\{ \mathbb{E}_{q_i (\cdot)}\left[\ln \bigg\{ \frac{1}{\sigma_i} \varphi \big(\frac{Z_i - \theta_i}{\sigma_i} \big) \bigg\} \right] - \textrm{KL} \big(q_i(\cdot) \parallel G \big) \bigg\}.
\end{align}
The two max-operations have an EM interpretation. Conditional on the choice of the `posteriors' $\{q_i(\cdot)\}_i$, maximization over the prior $G \in \mathcal{G}$ is equivalent to the M-step in an EM algorithm. The maximization over $\{q_i(\cdot)\}_i$ given $G$, which involving computing the posterior corresponding to the prior $G$, is equivalent to the E-step.
At the $g$-modeling estimate \(\hat{G}\), the E-step and M-step updates constitute a fixed point. At this fixed point, the E-step ensures that the distributions \(\{q_i(\cdot)\}_i\) are the posterior distributions \(\{q^*_{i,\hat{G}}(\cdot)\}_i\) corresponding to \(\hat{G}\). Consequently, the M-step update at this fixed point implies that
\begin{align}
\hat{G} &= \operatorname*{arg\,min}_{G \in \mathcal{G}} \frac{1}{n} \sum_{i=1}^{n}\textrm{KL}\left(q^*_{i,\hat{G}}(\cdot)\mid\mid G \right) \nonumber \\
& =\operatorname*{arg\,min}_{G \in \mathcal{G}}-\int\left(\frac{1}{n}\sum_{i=1}^{n}q^*_{i,\hat{G}}(\theta)\right)\ln g(\theta)d\nu(\theta) \nonumber \\
& =\operatorname*{arg\,min}_{G \in \mathcal{G}}\textrm{KL}\left(\bar{q}_{\hat{G}}(\cdot)\mid\mid G\right)
\label{M_step_problem}
\end{align}
where $g(\cdot)$ denotes the density of $G\in \mathcal{G}$ with respect to some dominating measure $\nu$, and
\[
\bar{q}_G(\cdot):=\frac{1}{n}\sum_{i=1}^{n}q^*_{i,G}(\cdot)
\]
denotes the average posterior given some prior $G$. We now explore the properties of the solution in two distinct settings: (1) when $\mathcal{G}$ is an exponential family, and (2) when $\mathcal{G}$ is unrestricted.
\subsubsection{Exponential family of priors}
For an exponential family of priors \( \mathcal{G} \), the solution to (\ref{M_step_problem}) is identified by matching the moments of the sufficient statistics \( u(\theta) \) between the prior and the average posterior distribution. This property, known as \emph{moment matching}, is a fundamental characteristic of exponential families \citep{bishop2006pattern}. Specifically, this condition implies that
\[
\mathbb{E}_{\hat{G}}[u(\theta)] = \mathbb{E}_{\bar{q}_{\hat{G}}}[u(\theta)],
\]
where \( \mathbb{E}_{\hat{G}} \) denotes the expectation with respect to the prior \( \hat{G} \), and \( \mathbb{E}_{\bar{q}_{\hat{G}}} \) denotes the expectation with respect to the average posterior \( \bar{q}_{\hat{G}} \).
For the Gaussian prior family \( \mathcal{G} = \{\mathcal{N}(0, \gamma^{-1}) : \gamma > 0 \} \), with a zero mean, the sufficient statistic is \( u(\theta) = \theta^2 \). The posterior distribution of \( \theta_i \) given \( \gamma \) and data is
\[
\theta_i \mid \mathcal{D} \sim \mathcal{N}\left(\frac{Z_i}{1 + \gamma \sigma_i^2}, \frac{\sigma_i^2}{1 + \gamma \sigma_i^2}\right).
\]
Moment matching for \( \theta^2 \) between the prior and the average posterior implies that the Empirical Bayes estimate \( \hat{\gamma} \) of the prior precision \( \gamma \) satisfies the equation
\begin{align}
\protect \label{eq:sample_moment_Gaussian_prior}
\frac{1}{n} \sum_{i=1}^n m(Z_i, \sigma_i; \hat{\gamma}) &=0, \quad \textrm{where} \\
m(Z_i, \sigma_i; \gamma) &:= \left( \frac{Z_i}{1 + \gamma \sigma_i^2} \right)^2 + \frac{\sigma_i^2}{1 + \gamma \sigma_i^2} - \frac{1}{\gamma} \nonumber
\end{align}
The above expression is the sample counterpart of the law of iterated expectations involving $\theta_i^2$, which states that the true prior \( G_0 := \mathcal{N}(0, \gamma_0^{-1}) \) satisfies
$\mathbb{E}_{G_0}\left[\mathbb{E}_{G_0}[\theta_i^2 \mid \mathcal{D}]\right] = \mathbb{E}_{G_0}[\theta_i^2]$.
This condition leads to the following necessary population moment condition for the true prior precision \( \gamma_0 \):
\begin{align}
\protect \label{eq:population_moment_Gaussian_prior}
\mathbb{E}_{G_0^{(n)}}\left[\frac{1}{n} \sum_{i=1}^n m(Z_i, \sigma_i; \gamma_0)\right] = 0,
\end{align}
where $\mathbb{E}_{G_0^{(n)}}[\cdot]$ denotes the marginal over $\mathcal{D} $ given the product prior $G_0^{(n)}$. Clearly, (\ref{eq:sample_moment_Gaussian_prior}) serves as the method of moments counterpart to the population moment condition (\ref{eq:population_moment_Gaussian_prior}). The population moment is only a necessary score condition. Showing that it vanishes at $\gamma_0$ does not prove uniqueness of its root. Identification instead follows from Lemma \ref{lem:prior-identification} in Appendix \ref{sec:Appendix:A}, which shows that the joint distribution of $(Z_i, \sigma_i)$ identifies the entire mixing distribution and hence identifies $\gamma_0$ in the Gaussian subfamily.
\subsubsection{Unconstrained family of priors}
When the class \(\mathcal{G}\) is unrestricted, as in the NPMLE procedure, (\ref{M_step_problem}) implies that \(\hat{G}\) satisfies a self-consistency property: it equals the average posterior distribution corresponding to itself. Formally, \(\hat{G} = \bar{q}_{\hat{G}}\).
This expression represents the sample counterpart of the `martingale' property of Bayesian updating, which asserts that the expected posterior must equal the prior. This in turn follows from the law of iterated expectations, which states that for any measurable function \(h(\cdot)\), $\mathbb{E}_{G_0}\left[\mathbb{E}_{G_0}[h(\theta_i) \mid \mathcal{D}]\right] = \mathbb{E}_{G_0}[h(\theta_i)].$
Of course this is again a necessary fixed-point characterization, not a standalone identification result. For example, every point mass is unchanged by Bayesian updating when it is itself used as the candidate prior, irrespective of the true data law. Identification is established in Lemma \ref{lem:prior-identification} in Appendix \ref{sec:Appendix:A}.
\subsubsection{Taking stock}
The above discussion explains why \(g\)-modeling is still valid even when we substitute the true marginal likelihood \( p_G(\mathcal{D}) \) with the working marginal likelihood \( \prod_i f_{G,\sigma_i}(Z_i) \). In any adaptive experiment, Bayesian updating ensures that the sequence of posterior distributions forms a martingale (in the space of probability measures). Consequently, the prior equals the expected posterior, even under adaptive stopping; this is a basic informational constraint--Bayes consistency--that every experiment must obey. Viewed this way, \( g \)-modeling amounts to applying a sample-based version of Bayesian consistency, which delivers identifying restrictions for the true prior via the law of iterated expectations. Because the posterior is invariant to the choice of algorithm by the likelihood principle, the moment-matching interpretation implies that \( g \)-modeling remains valid irrespective of the specific algorithm employed.
For parametric families \( \mathcal{G} \), the regret consistency of \( \hat{G} \) is an immediate consequence of standard generalized method-of-moments (GMM) reasoning. Demonstrating regret consistency for the nonparametric maximum likelihood estimator (NPMLE) is more involved, since it corresponds to satisfying an uncountable collection of moment conditions (see Section \ref{sec:Theoretical results} for the rigorous analysis). However, the basic intuition is essentially the same.
\@startsection{subsection}{2}{\z@}
{0pt}
{0.4ex}
{\normalfont\normalsize\bfseries}{Tweedie's formula and the failure of naive $f$-modeling}
The celebrated Tweedie's formula establishes a connection between the posterior expectation of $\theta_i$ and the derivative of the marginal working likelihood. Specifically, recalling the definition of $f_{G,\sigma_i}(Z_i)$ from (\ref{eq:definition_of_f_G}), some straightforward algebra gives
\begin{align} \label{Tweedie's formula}
\sigma_i^2 \nabla_z\ln f_{G,\sigma _i}(z) \big\rvert_{z=Z_i} &= \frac{\int (\theta_i - Z_i) \varphi \left(\frac{Z_i - \theta_i}{\sigma_i} \right)dG(\theta_i)}{\int \varphi \left(\frac{Z_i - \theta_i}{\sigma_i} \right)dG(\theta_i)} =\mathbb{E}_{G}[\theta_i|\mathcal{D} ] - Z_i.
\end{align}
Tweedie's formula lies at the core of the Empirical Bayes \( f \)-modeling approach, which seeks to estimate the `true' posterior mean \( \mathbb{E}_{G_0}[\theta_i \mid \mathcal{D}] \) using a non-parametric approximation to \( \nabla_z \ln f_{G_0,\sigma_i}(z) \). In the classical exogenous sampling framework, \(f_{G_0,\sigma_i}(z) \) equals the marginal density, $p(Z_i)$, of \( Z_i \), which serves as a sufficient statistic for the data. Consequently, \( \nabla_z \ln f_{G_0,\sigma_i}(z) \) can be non-parametrically estimated as the derivative of the log-marginal density $\ln p(Z_i)$ of \( Z_i \).
This approach, however, breaks down in the context of adaptive experimentation. Under adaptive sampling, $ f_{G_0,\sigma_i}(z)$ no longer equals the true marginal density of $Z_i$ because the true conditional distribution of $Z_i$ given $\theta_i$ is not really $\mathcal{N}(\theta_i, \sigma_i^2)$. In fact, as shown in \cite{adusumilli2021risk}, $Z_i$ even ceases to be a sufficient statistic for the data in adaptive settings. At a minimum, both $Z_i$ and $N_i$ are needed for Bayesian updating. As a result, the derivative of $\ln p(Z_i)$ need not estimate \( \nabla_z \ln f_{G_0,\sigma_i}(Z_i) \), rendering the \( f \)-modeling approach invalid under adaptive sampling.
\@startsection{subsection}{2}{\z@}
{0pt}
{0.4ex}
{\normalfont\normalsize\bfseries}{Heteroskedastic models and CLOSE} \protect \label{subsec:CLOSE}
Recently, \cite{chen2022empirical} proposed the CLOSE framework for analyzing heteroskedastic models of the form (\ref{eq:working_likelihood}). This approach assumes a location-scale dependence between \(\theta_i\) and \(\sigma_i\):
\[
\frac{\theta_i - \mathbb{E}[\theta_i \vert \sigma_i]}{\textrm{sd}[\theta_i \vert \sigma_i]} \Bigg\vert \sigma_i \sim G, \quad \text{independently of } \sigma_i.
\]
This effectively constrains \(p(\theta_i | \sigma_i)\), the conditional density of \(\theta_i\) given \(\sigma_i\).
In our setting, however, \(\sigma_i\) depends on the data through \(\tau_i\), implying that \(p(\theta_i | \sigma_i)\) represents a partial posterior, where only a subset of the data (\(\tau_i\)) has been used for Bayesian updating. In adaptive experiments, \(\tau_i\) is determined in a complex manner by \(\theta_i\), so there is no reason to expect this assumption to hold.
More critically, the CLOSE framework first estimates \(p(\theta_i | \sigma_i)\) and then treats it as a prior for an additional Bayesian update over the data. However, since computing \(p(\theta_i | \sigma_i)\) already involves updating over part of the data, this results in Bayesian updating being applied twice. Repeated updating on the same data will ultimately recover the marginal distribution of \(Z_i\). We would therefore expect CLOSE to lead to under-shrunk in the context of compound adaptive experiments, as compared to standard \(g\)-modeling.\footnote{This concern does not apply to the typical use-case of CLOSE, where $\tau_i$ is determined \textit{prior} to the experiment and may be correlated with $\theta_i$. In that scenario $\sigma_i$ is independent of the data given $\theta_i$.}
\@startsection{subsection}{2}{\z@}
{0pt}
{0.4ex}
{\normalfont\normalsize\bfseries}{Efficiency of $g$-modeling}
Because $g$-modeling is a maximum-likelihood method over the space of priors, we expect it to be efficient when nothing is known about the adaptive algorithms in use. In fact, information about these algorithms would enhance the estimation of $G_0$ only if the algorithms were explicit, known functions of the true $G_0$. In realistic settings, however, characterizing how the algorithms depend on $G_0$ would be cumbersome, if not infeasible. Thus, any loss of efficiency arising from ignorance of the algorithms is likely to be negligible.
\section{Theoretical Results\protect\label{sec:Theoretical results}}
In this section, we examine the theoretical properties of $g$-modeling under two distinct scenarios: (1) when the class of priors is restricted to be Gaussian, and (2) when the class of priors is unrestricted, as in the NPMLE approach. For both scenarios, our primary focus is to evaluate the performance of $g$-modeling methods in terms of the average Bayes risk criterion, comparing their efficacy against an oracle benchmark that has full knowledge of the true prior distribution, $G_0$.
We start with a formal definition of compound Bayes risk and regret.
\@startsection{subsection}{2}{\z@}
{0pt}
{0.4ex}
{\normalfont\normalsize\bfseries}{Compound Bayes risk and regret}
Recall that our setup consists of $n$ experiments with individual treatment effects $\bm{\theta} := (\theta_1, \dots ,\theta_n)$ that we aim to estimate. We view the parameters $\theta_i$ as i.i.d draws from an unknown prior $G_0$.
Let $\delta_i := \delta_{\theta}(\mathcal{D} _i)$ denote an estimator of $\theta_i$, and $\bm{\delta} := (\delta_1, \dots, \delta_n)$ the collection of estimators for each experiment. We define compound frequentist risk as
$$
R(\bm{\delta}, \bm{\theta}) = \mathbb{E} \left[ \frac{1}{n} \sum_{i=1}^n \vert \delta_i - \theta_i \vert^2 \ \right\vert \bm{\theta} \bigg].
$$
The compound Bayes risk criterion integrates $R(\bm{\delta}, \bm{\theta})$ over the joint prior $\bm{\theta} \sim G_0^{(n)}$:
$$
R(\bm{\delta}, G_0) = \mathbb{E}_{G_0^{(n)}} \left[ \frac{1}{n} \sum_{i=1}^n \vert \delta_i - \theta_i \vert^2 \ \right].
$$
Following \cite{jiang2020general}, we evaluate estimators $\bm{\delta}$ using the compound Bayes risk criterion. Consider an oracle who knows the true prior $G_0$. The oracle estimator of $\theta_i$ is clearly $\delta_i^* := \mathbb{E}_{G_0}[\theta_i \vert \mathcal{D} _i]$ and write $\bm \delta^* :=(\delta_1^*, \dots, \delta_n^*)$. The difference in compound Bayes risk between the oracle estimator and the proposed estimator $\bm{\delta}$ is known as regret:
$$
\mathcal{R}(\bm{\delta},G_{0})=R(\bm{\delta},G_{0})-R(\bm{\delta}^{*},G_{0}).
$$
From the form of $\delta_{i}^{*}$, some simple algebra indicates that
\begin{equation}
\protect \label{eq:regret formula}
\mathcal{R}(\bm{\delta},G_{0}) = \mathbb{E}_{G_0^{(n)}}\left[ \frac{1}{n} \sum_{i=1}^n \left\vert \delta_{i}-\delta_{i}^{*}\right\vert ^{2}\right].
\end{equation}
We call an estimator $\bm{\delta}$ regret consistent if the regret goes to 0 asymptotically as $n\to \infty$. Note that the MLE estimator of $\theta_i$ is given by $\delta_i^\textrm{mle} = Z_i$.
\@startsection{subsection}{2}{\z@}
{0pt}
{0.4ex}
{\normalfont\normalsize\bfseries}{Regret consistency of $g$-modeling with Gaussian priors\protect\label{sec:Theoretical results: Gaussian}}
We now establish results on regret consistency for \( g \)-modeling with a Gaussian prior family. Specifically, consider the Gaussian prior family defined as \(\mathcal{G} = \{\mathcal{N}(0, \gamma^{-1}): \gamma > 0 \}\), and assume that the true prior, \( G_0 \), belongs to this family, with \( G_0 = \mathcal{N}(0, \gamma_0^{-1}) \). The empirical Bayes (EB) estimate, \(\hat{\gamma}\), is then obtained by solving the corresponding sample moment equation provided in \eqref{eq:sample_moment_Gaussian_prior}. The \( g \)-modeling estimator for \(\theta_i\) is the posterior mean under the estimated prior:
\[
\hat{\delta}_i^{\text{EB}} = \frac{Z_i}{1 + \hat{\gamma} \sigma_i^2}.
\]
\subsubsection{Leave-one-out estimation}
We are able to derive a remarkably simple proof of regret consistency when the experiments are independent of each other, i.e., $(Z_i, \tau_i)$ are independent across $i$, and \(\gamma_0\) is estimated using a leave-one-out (LOO) methodology. Let \(\hat{\gamma}^{(-i)}\) denote the leave-one-out estimate of \(\gamma\), computed by excluding the \(i\)-th observation from the sample. The corresponding leave-one-out EB estimator is then denoted:
\[
\tilde{\delta}_i^{\text{EB}} = \frac{Z_i}{1 + \hat{\gamma}^{(-i)} \sigma_i^2}.
\]
To assess the performance of this estimator, we consider its regret ratio relative to the MLE estimator. Assume that experiments are indepdendent of each other,
\begin{align}
\label{eq:regret_ratio}
\frac{\mathcal{R}(\tilde{\bm{\delta}}^{\textrm{EB}},G_{0})}{\mathcal{R}(\bm{\delta}^{\textrm{mle}},G_{0})} & = \nicefrac{\sum_i \mathbb{E}_{G_{0}}\left[ \left\vert \frac{Z_{i}}{1 +\hat{\gamma}^{(-i)}\sigma _i^2 }-\frac{Z_{i}}{1+\gamma_{0}\sigma_i^2 }\right\vert ^{2}\right]}
{\sum_i \mathbb{E}_{G_{0}}\left[ \left\vert Z_i -\frac{Z_{i}}{1+\gamma_{0}\sigma_i^2 }\right\vert ^{2}\right]}
\nonumber \\
& =\nicefrac{\sum_i \mathbb{E}_{G_{0}}\left[\frac{Z_{i}^{2}\left(\hat{\gamma}^{(-i)}-\gamma_{0}\right)^{2} \sigma_i^4 }{\left(1 +\hat{\gamma}^{(-i)}\sigma _i^2 \right)^{2}\left(1+\gamma_{0}\sigma_i^2 \right)^{2}}\right]}
{\sum_i \mathbb{E}_{G_{0}}\left[\frac{Z_{i}^{2}\gamma_{0}^{2} \sigma_i^4}{\left(1+\gamma_{0} \sigma_i^2 \right)^{2}}\right]} \nonumber \\
& \le \nicefrac{\sum_i \mathbb{E}_{G_{0}}\left[\frac{Z_{i}^{2}\left(\hat{\gamma}^{(-i)}-\gamma_{0}\right)^{2} \sigma_i^4 }{\left(1+\gamma_{0}\sigma_i^2 \right)^{2}}\right]}
{\sum_i \mathbb{E}_{G_{0}}\left[\frac{Z_{i}^{2}\gamma_{0}^{2} \sigma_i^4}{\left(1+\gamma_{0} \sigma_i^2 \right)^{2}}\right]}
\le \sup_{1\le i \le n} \mathbb{E}_{G_{0}}\left[\left(\frac{\hat{\gamma}^{(-i)}-\gamma_{0}}{\gamma_{0}}\right)^{2}\right],
\end{align}
where the first equality is due to (\ref{eq:regret formula}), the inequality is due to $\hat{\gamma}^{(-i)} > 0$, and the last inequality follows from the fact $\hat{\gamma}^{(-i)}$ is independent of $Z_{i},\sigma_{i}$ under the assumption that the experiments are independent of each other.
Remarkably, (\ref{eq:regret_ratio}) is both algorithm-independent and non-asymptotic. Since $\hat{\gamma}$ solves a sample moment condition based on $n$ observations, standard regularity conditions ensure that for all $i$:
\begin{equation}
\label{eq:gamma_bound}
\underset{1 \leq i \leq n}{\sup} \mathbb{E}_{G_0}\left[\left(\frac{\hat{\gamma}^{(-i)} - \gamma_0}{\gamma_0}\right)^2\right] = O(n^{-1}).
\end{equation}
Substituting this result into (\ref{eq:regret_ratio}), we find:
\begin{equation}
\label{eq:regret_bound_gaussian}
\mathcal{R}(\tilde{\bm{\delta}}^{\textrm{EB}}, G_0) = O(n^{-1}) \cdot \mathcal{R}(\bm{\delta}^{\textrm{mle}}, G_0),
\end{equation}
implying that the regret of the leave-one-out EB estimator is a vanishingly small fraction of the MLE regret as $n \to \infty$.
In fact, in many practical scenarios, the regret of the MLE estimator,
$$
\mathcal{R}(\bm{\delta}^{\textrm{mle}}, G_0) = \frac{1}{n} \sum_i \mathbb{E}_{G_{0}}\left[\frac{Z_{i}^{2}\gamma_{0}^{2} \sigma_i^4}{\left(1+\gamma_{0} \sigma_i^2 \right)^{2}}\right] < \frac{1}{n} \sum_i \mathbb{E}_{G_{0}}[Z_i^2]
$$
is finite whenever $n^{-1} \sum_i \mathbb{E}_{G_0}[Z_i^2] = O(1)$. Under this additional moment condition, (\ref{eq:regret_bound_gaussian}) further establishes that the regret of $\tilde{\delta}_i^{\textrm{EB}}$ decays at an $O(n^{-1})$ rate.
\subsubsection{General results}
The assumption that experiments are independent may be violated in practice, and the computational cost of leave-one-out estimation can also be substantial. Still, even when these conditions fail, we can derive the following bound on the regret of the empirical Bayes (EB) estimator:
\[
\begin{aligned}
\mathcal{R}(\hat{\bm{\delta}}^{\textrm{EB}},G_{0})
& = \frac{1}{n} \sum_{i=1}^n \mathbb{E}_{G_{0}}\!\left[ \left| \frac{Z_{i}}{1 + \hat{\gamma}\sigma_i^2 } - \frac{Z_{i}}{1 + \gamma_0\sigma_i^2 } \right|^2 \right] \le \mathbb{E}_{G_{0}}\!\left[ \left( \frac{1}{n} \sum_{i=1}^n \frac{Z_i^2}{\sigma_i^4} \right) \left( \frac{1}{\hat{\gamma}} - \frac{1}{\gamma_0} \right)^2 \right].
\end{aligned}
\]
As a result, we can still show that $\mathcal{R}(\hat{\bm{\delta}}^{\textrm{EB}},G_{0}) = O(n^{-1})$,
albeit under stronger regularity conditions (we skip the formal statement for brevity).
\@startsection{subsection}{2}{\z@}
{0pt}
{0.4ex}
{\normalfont\normalsize\bfseries}{Regret consistency of NPMLE}\label{subsec:regret_consistency_NPMLE}
To study the regret consistency of NPMLE, we focus on the setting where the experiments are independent, so that \((Z_i,\tau_i)\) are independent across \(i\). Our treatment of the NPMLE follows \cite{jiang2009general} and \cite{jiang2020general}. Define
\[
\ell_n(G) = \frac{1}{n} \sum_i \ln f_{G,\sigma_i}(Z_i),
\]
with $\sigma_i^2 := \omega_i^2/\tau_i$ and let the NPMLE of \(G_0\) be given by
\begin{equation} \label{NPMLE}
\hat{G}_n \in \underset{G \in \mathcal{G}(\mathcal{A}_n)}{\operatorname*{arg\,max}}\; \ell_n(G),
\end{equation}
where \(\mathcal{G}(\mathcal{A}_n)\) is the collection of probability measures supported on a finite, possibly data dependent set \(\mathcal{A}_n \subseteq \mathbb{R}\). In applications, the finiteness of \(\mathcal{A}_n\) facilitates the use of convex optimization \citep{koenker2014convex} to compute the NPMLE.
Let \(\bm{\delta}^{NPEB}\) denote the plug-in estimator of the optimal Bayes rule obtained by substituting the NPMLE \(\hat{G}_n\) for the true prior \(G_0\). Concretely,
\[
\delta_i^{NPEB}
= \frac{\int \theta \,\frac{1}{\sigma_i}\, \varphi\!\left(\frac{Z_i - \theta}{\sigma_i}\right) d\hat{G}_n}
{\int \frac{1}{\sigma_i}\, \varphi\!\left(\frac{Z_i - \theta}{\sigma_i}\right) d\hat{G}_n}.
\]
The central theoretical claim of this section is that \(\bm{\delta}^{NPEB}\) achieves regret consistency. To establish this, we impose the assumptions stated below. Let \(\bm{Z} := (Z_1,\dots,Z_n)\) and \(\bm{\tau} := (\tau_1,\dots,\tau_n)\). We write \(\lesssim\) to indicate an inequality that holds up to a multiplicative constant independent of \(n\).
\begin{assumption}\label{boundedG0}
$\theta_1, \dots, \theta_n \sim _{iid} G_0$, where $G_0$ has bounded support on $[-R,R]$. Furthermore, $(Z_i,\tau_i)$ are independent across $i$.
\end{assumption}
\begin{assumption}\label{boundVar}
There exist $\omega_\ell >0$ and $\omega_u < \infty$ independent of $n$ such that $0 < \omega_\ell^2 \le \omega_i^2 \le \omega_u^2< \infty$ for all $i$. Additionally,
\[
\mathbb P\left(
\underline\tau\leq\tau_i\leq\bar\tau
\text{ for every }i=1,\ldots,n
\right)=1.
\]
Consequently, with $\sigma_\ell^2:=\frac{\omega_\ell^2}{\bar\tau}$ and $\sigma_u^2:=\frac{\omega_u^2}{\underline\tau}$, the random variances $\sigma_i^2:=\frac{\omega_i^2}{\tau_i}$ satisfy
\[
\mathbb P\left(
\sigma_\ell^2\leq\sigma_i^2\leq\sigma_u^2
\text{ for every }i=1,\ldots,n
\right)=1.
\]
\end{assumption}
\begin{assumption}\label{likfactor}
The likelihood of the data collected across $n$ experiments takes the form
\[
p(\mathcal{D} | \theta_1, \dots, \theta_n) = c(\mathcal{D} ) \cdot \prod_{i=1}^n \frac{1}{\sigma_i} \varphi\left(\frac{Z_i - \theta_i}{\sigma_i}\right),
\]
for some $c(\mathcal{D} )$ independent of $\theta_1, \dots, \theta_n$. Furthermore, there exist $\bar{c} < \infty$ such that
$$
\sup_{z_i, \tau_i \in [\underline{\tau}, \bar \tau]} \left\{ \sqrt{2\pi} \sigma_i e^{z_i^2/2\sigma_i^2} \cdot p(z_i,\tau_i | \theta_i = 0) \right\} \le \bar{c}.
$$
Here $p(z, \tau | \theta)$ denotes the density with respect to $m(z)\otimes\nu(\tau)$, where $m(\cdot)$ denotes the Lebesgue measure, and $\nu(\tau)$ is some probability measure over the support of $\tau$.
\end{assumption}
\begin{assumption}\label{tailbound}
Let $\bar Z_n = \underset{i}{\max} |Z_i| \vee 1$, then
$\mathbb{E}[\bar Z_n^4] \lesssim (\ln n)^2$ and $\mathbb{P}(\bar Z_n \geq M_n) \lesssim \frac{1}{n^2}$
with $M_n = \sqrt{\kappa \ln n}$ for some $\kappa>0$.
\end{assumption}
\begin{assumption}\label{likbound}
Let $q_n = \Big(\frac{e \sqrt{2\pi}}{n^2}\Big) \land 1$ and $\kappa_n = \frac{1}{n} \ln \frac{1}{q_n}$. For all sufficiently large $n$, the NPMLE $\hat G_n$ satisfies, almost surely,
\[
\underset{G \in \mathcal{G}(\mathbb{R})}{\sup} \ell_n(G) - \ell_n(\hat G_n) \leq \kappa_n.
\]
\end{assumption}
\begin{assumption}\label{support-geometry}
Let $r_i(z,\tau) = p(z,\tau | \theta_i = 0)$. For every $i$ and for $\nu$-almost every $\tau$, there is an open set
$\mathcal S_{i,\tau}\subseteq\mathbb R$ such that
\[
r_i(z,\tau)>0 \quad\text{for }z\in\mathcal S_{i,\tau},
\qquad
r_i(z,\tau)=0
\quad\text{for Lebesgue-a.e. }z\notin\mathcal S_{i,\tau}.
\]
The set $\{(z,\tau): z \in \mathcal{S}_{i,\tau}\}$ is measurable. There exists $\ell_*>0$, independent of $n,i,\tau$, such that every
bounded connected component $I$ of $\mathcal S_{i,\tau}$ satisfies
\[
|I|\ge\ell_*.
\]
\end{assumption}
\begin{assumption}\label{density-2}
Let $r_i(z,\tau) = p(z,\tau | \theta_i = 0)$. There is a fixed integer $m\ge 2$ and constants
$C_1,\ldots,C_m<\infty$, independent of $n$ and $i$, such that, for
$\nu$-almost every $\tau$, the function $z\mapsto \ln r_i(z,\tau)$ is
$m$-times continuously differentiable on each connected component $I$ of $\mathcal{S}_{i,\tau}$, and
\begin{equation}\label{eq:smoothness_assumption}
\left|\partial_z^l\ln r_i(z,\tau)\right|
\le C_l(1+|z|)^l, \qquad z \in I,
\qquad l=1,\ldots,m.
\end{equation}
\end{assumption}
Assumption \ref{boundedG0} stipulates that $G_0$ has compact support.\footnote{This can be relaxed to a tail assumption of $G_0$ at the cost of more burdensome notations for the theoretical results to follow.} Assumption \ref{boundVar} requires the stopping time $\tau_i$ in each experiment to be bounded and bounded away from zero almost surely, uniformly over \(i\), with deterministic bounds independent of \(n\). This condition is met by many algorithms. Moreover, each $\omega_i$ is finite and bounded away from zero, a standard requirement in the empirical Bayes literature; see, e.g., \cite{jiang2020general}, \cite{soloff2024multivariate}, and \cite{chen2022empirical}.
The first part of Assumption \ref{likfactor} reiterates the likelihood principle (\ref{likelihoodPrinciple:generalization}), which we have already established in Proposition \ref{Prop:likelihood principle}. The second part of Assumption \ref{likfactor} is more substantive. It requires that the conditional density of $(Z_i, \tau_i)$ exhibit Gaussian tails when $\theta_i = 0$. Under this condition, Lemma \ref{lem:density_dominance} in Section \ref{sec:Appendix:SB} of the online supplement demonstrates that the marginal density of $(Z_i, \tau_i)$ can be bounded above by a multiplicative constant times the working marginal density which is useful in establishing the result in Theorem \ref{thm:NPMLE-consistency}.
Lemma \ref{lemma:sufficient_condition_Asm3} in Appendix \ref{appendixC} demonstrates a bound on $\bar c$ when $\tau$ has finite support.
Assumption \ref{likbound} requires $\hat G_n$ to be an
approximate maximizer of the likelihood over the full class
$\mathcal G(\mathbb R)$.
It is used twice to establish the theoretical guarantee. It supplies the
true prior likelihood lower bound used to prove Proposition
\ref{HellingerAccuracy} in Section \ref{sec:Appendix:Hellinger} of the online supplement, and it also supplies the approximate generalized
MLE condition for Theorem 5 of \cite{jiang2020general}, both of which channels towards establishing Theorem \ref{thm:NPMLE-consistency}. Lemma \ref{lemma: qualityNPMLE} in the Appendix \ref{appendixC} constructs a data-dependent grid $\mathcal{A}_n$
using which the NPMLE satisfies Assumption \ref{likbound}.
Assumption \ref{support-geometry} places a mild restriction on the connected set over which the conditional density of $(Z_i, \tau_i)$ when $\theta_i=0$, denoted as $r_i(z,\tau)$, is strictly positive. If $r_i(z, \tau)$ is strictly positive in $z$ for $\nu$-almost every $\tau$, then we can take $\mathcal{S}_{i,\tau} = \mathbb{R}$. In this case there are no bounded connected components, so Assumption \ref{support-geometry} always holds. If instead $r_i(z,\tau)$ vanishes on part of the $z$-spce, the assumption prevents any of its bounded connected components from being arbitrarily short as $n, i,$ or $\tau$ varies. The lower-length conditon is needed to keep a lower-order term in the Gagliardo-Nirenberg interpolation inequality uniformly controlled when applied to each connected sub-component; for details, see the proof of Lemma \ref{GNinequality}.
Assumption \ref{density-2} concerns the smoothness properties of $r_i(z,\tau)$. In particular, it requires polynomial control of its first
$m$ log-derivatives wherever that density is positive. The integer $m$ then governs the rate at which the Bayes regret converges to 0.
Finally to show all these assumption can hold for adaptive experiments, as a demonstrating example, we verify all the Assumptions for group sequential trials with two stages in Appendix \ref{sec:two-stage-verification}.
Under the above assumptions, we show that the plug-in Bayes estimator \( \bm{\delta}^{NPEB} \) is regret-consistent.
\begin{thm} \label{thm:NPMLE-consistency}
Under Assumptions \ref{boundedG0}-\ref{density-2}, for some $m \geq 2$,
\[
\mathcal{R}(\bm \delta^{NPEB}, G_0) \lesssim (\ln n) \cdot\Big( \frac{(\ln n)^2}{n} \Big)^{(1-\frac{1}{m})}.
\]
\end{thm}
\subsubsection{Proof sketch}
By Proposition \ref{HellingerAccuracy} in Section \ref{sec:Appendix:Hellinger} of the online supplement, which extends the findings of \cite{jiang2009general} and \cite{jiang2020general} to the adaptive framework, we obtain
\begin{equation}\label{eq:Hellinger_rate}
\mathbb{P}\!\left(
\bar h^2(p_{\widehat G_n},p_{G_0})\geq t^2\epsilon_n^2
\right)
\leq3n^{-t^2} \quad \forall\ t\ge 1,
\end{equation}
where $\epsilon_n^2=C_{\epsilon}(\ln n)^2/n$.
Because the common factor $c_i$ in
$p_{G,i}=c_i f_{G,\sigma_i}$ drops out of the density ratios, Tweedie's
formula, together with the bound $\sigma_i \le \sigma_u$ for all $i$, implies
$$
\mathcal R(\bm\delta^{NPEB},G_0)
=
\mathbb{E}\!\left[
\frac1n\sum_i\sigma_i^4
\left\{
\partial_z\ln \frac{p_{G_0,i}}{p_{\widehat G_n,i}}
(Z_i,\tau_i)
\right\}^2
\right]
\lesssim
D_f\left(p_{G_0,i}, p_{\hat G_n,i} \right),
$$
where $D_f(p,q) = E_p[(\nabla \ln(p/q))^2]$ denotes the Fisher divergence between $p$ and $q$. As in \cite{jiang2009general}, \cite{jiang2020general}, and \cite{soloff2024multivariate}, the central task is to relate the Fisher divergence to the Hellinger distance. However, the techniques in the aforementioned papers cannot be applied verbatim here since they crucially rely on properties particular to Gaussian convolutions. In our setting, the marginal density takes the form
\[
p_{G,i}(z,\tau)
=c_i(z,\tau)f_{G,\sigma_i(\tau)}(z).
\]
While the common factor $c_i$ vanishes in the density ratio, it remains present in the Hellinger distance and thereby affects the smoothness properties of $p_{G,i}(z,\tau)$. To handle this, we employ a different and more general approach based on Gagliardo–Nirenberg interpolation inequalities.
Set $u_i=\sqrt{p_{G_0,i}}$ and $v_i=\sqrt{p_{G,i}}$ for some fixed $G$. Using $\nabla \ln p = 2\, (\nabla \sqrt{p}/\sqrt{p})$ and some straightforward algebra,
\begin{align*}
D_f\left(p_{G_0,i}, p_{G,i} \right)
& \lesssim \int\left(\partial_z u_i- \partial_z v_i\right)^{2}
+ \int\left(\partial_z\ln p_{G_0, i}\right)^{2}\left(u_i-v_i\right)^{2}\\
& \lesssim \int\left(\partial_z u_i- \partial_z v_i\right)^{2}
+\left\Vert \partial_z\ln p_{G_0, i}\right\Vert _{\infty}^2 h^{2}(p_{G_0,i},p_{G,i}).
\end{align*}
Assume for simplicity that $\mathcal{S}_{i,\tau}$ in Assumption \ref{support-geometry} coincides with $\mathbb{R}$. Then, the Gagliardo-Nirenberg interpolation inequality states that
\begin{equation}\label{eq:GN-inequality}
\left\Vert \nabla g\right\Vert _{2} \lesssim \left\Vert g\right\Vert _{2}^{1-1/m}\left\Vert \partial^{m}g\right\Vert _{2}^{1/m},
\end{equation}
for any function $g$ and $m \geq 2$. Applying this inequality yields
\begin{equation} \label{eq:Fisher_Hellinger_bound}
\begin{aligned}
D_f\left(p_{G_0,i}, p_{G,i} \right)
& \lesssim \left\Vert \partial^{m}(u_i-v_i)\right\Vert_{2}^{1/m} [ h^{2}(p_{G_0,i},p_{G,i})]^{1-1/m} \\
& \quad + \left\Vert \partial\ln p_{G_0, i}\right\Vert_{\infty}^2 h^{2}(p_{G_0,i},p_{G,i}).
\end{aligned}
\end{equation}
Under Assumption \ref{density-2}, we can verify that both $\left\Vert \partial^{m}(u_i-v_i)\right\Vert_{2}$ and $\left\Vert \partial_z\ln p_{G_0, i}\right\Vert _{\infty}$ are bounded up to some $\ln n$ factors.
Suppose for the moment that we were allowed to replace $G$ with $\hat{G}_n$ in (\ref{eq:Fisher_Hellinger_bound}). In that case, combining this substitution with (\ref{eq:Hellinger_rate}) would yield a proof of Theorem \ref{thm:NPMLE-consistency}.
However, since $\widehat G_n$ is a random quantity estimated from the very same data used to evaluate its score, the deterministic inequality (\ref{eq:Fisher_Hellinger_bound}) cannot be invoked directly. To address this, one must instead rely on covering-number arguments and truncation of $Z_i$ to a compact set, as in \citet[Lemma 10]{soloff2024multivariate}. The detailed steps are provided in the proof of Theorem \ref{thm:NPMLE-consistency} in Section \ref{sec:Appendix:S} of the online suppment.
\subsubsection{Discussion}
The exponent in Theorem \ref{thm:NPMLE-consistency} reflects the smoothness order $m$ appearing in Assumption \ref{density-2}. In the standard
Gaussian normal-means model with fixed exogenous scales, the marginal
density of $Z_i$ is $f_{G,\sigma_i}$. The Gaussian base density satisfies
Assumption \ref{density-2} for any fixed finite $m$. Hence, for every
fixed $m \ge 2$, Theorem \ref{thm:NPMLE-consistency} yields
\[
\mathcal R(\bm\delta^{NPEB},G_0)
\lesssim
n^{-1+1/m}(\ln n)^{3-2/m}.
\]
Therefore, for any fixed $\varepsilon>0$, choosing
$m$ sufficiently large
(but still fixed) leads to the rate
$O\!\left(n^{-1+\varepsilon}(\ln n)^3\right)$.
In the classical Gaussian setting, \cite{jiang2009general} and
\cite{soloff2024multivariate} obtain near-parametric NPMLE
regret upper bounds of order $n^{-1}$, up to logarithmic factors, under
compact-support or related tail assumptions. While the formal limit
$m\to\infty$ of our rate has the same $n^{-1}$ order (up to log factors), this limit does not
rigorously follow
from Theorem \ref{thm:NPMLE-consistency}, since the
constants
$\{C_l\}_{l=1}^m$ in Assumption \ref{density-2}, as well as those
entering the Gagliardo--Nirenberg interpolation inequalities, may depend
on $m$.
Obtaining a bound that is uniform over $m \in [2,\infty)$, and
clarifying whether the
rate for finite $m$ is optimal, is left to future
research.
\section{Extensions}
\@startsection{subsection}{2}{\z@}
{0pt}
{0.4ex}
{\normalfont\normalsize\bfseries}{Parametric models and local asymptotics}\protect \label{subsec:local asymptotics}
Thus far, our $g$-modeling framework has relied heavily on the assumption of Gaussian outcomes. However, as we demonstrate below, the Gaussian likelihood naturally emerges as an approximation to the true likelihood within a local asymptotic regime.
Consider the multi-arm sequential experiment setup described in Section \ref{subsec: Generalizing}, but now assume that the outcomes from each arm follow a parametric distribution:
\[
Y_{j,i} \sim P_{\mu_i},
\]
where \( P_{\mu} \) is a known parametric family and \( \mu_i \) is a scalar experiment-specific parameter. For instance, in A/B testing applications, outcomes are often binary, in which case \( Y_{j,i} \sim \text{Bernoulli}(\mu_i) \). As before, let \( \tau_i \) denote the proportion of times arm \( i \) was sampled, relative to the total sample size \( N \).
Following the standard local asymptotic framework, we assume that each \( \mu_i \) represents a local perturbation of a reference parameter \( \mu_0 \):
\[
\mu_i = \mu_0 + \frac{\theta_i}{\sqrt{N}}.
\]
Furthermore, as in the earlier sections, we also suppose that \( \theta_1, \dots, \theta_n \) are independent and identically distributed (i.i.d.) draws from a common prior \( G_0 \).
We emphasize that the reparametrization of \( \mu_i \) in terms of \( \theta_i \) serves primarily to provide a theoretical justification for our procedures. In practice, knowledge of \( \mu_0 \) is not required, as placing a prior on \( \mu \) is equivalent to placing a prior on \( \theta \).
To ensure a well-posed asymptotic analysis, we impose the standard assumption that the parametric family \( P_{\mu} \) is quadratic mean differentiable (QMD) around the reference point \( \mu_0 \). This condition guarantees the existence of a well-defined score function \( \psi(Y_{j,i}) \) and an associated inverse Fisher information given by $\omega_i^2 := \mathbb{E}_0[\psi^2]^{-1}$, where \( \mathbb{E}_0[\cdot] \) denotes expectation under \(P_0 := P_{\mu_0} \).
We define the (normalized) average score in experiment \( i \) as
\[
Z_i = \frac{\sqrt{N} \, \omega_i^2 }{\lfloor N \tau_i \rfloor} \sum_{j=1}^{\lfloor N \tau_i \rfloor} \psi(Y_{j,i}).
\]
Let \( p(\mathcal{D} \mid \theta_1, \dots, \theta_n) \) denote the likelihood of the data \( \mathcal{D} \) given the parameter vector \( (\theta_1, \dots, \theta_n) \). By Lemma 2 in \cite{adusumilli2021risk}, under the QMD property and provided that \( \max_i \tau_i \) is bounded, the likelihood ratio of \( (\theta_1, \dots, \theta_n) \) relative to the reference \( (0,\dots,0) \) is approximated as
\[
\ln \frac{p(\mathcal{D} \mid \theta_1, \dots, \theta_n)}{p(\mathcal{D} \mid 0, \dots, 0)} = \sum_i \left\{\frac{1}{\sigma_i^2} Z_i \theta_i - \frac{1}{2\sigma_i^2} \theta_i^2 \right\} + o_{P_0}(1),
\]
where \( \sigma_i^2 := \omega_i^2/\tau_i \).
Consequently, up to an asymptotically negligible error, the log-likelihood can be expressed as
\begin{align*}
\ln p(\mathcal{D} \mid \theta_1, \dots, \theta_n) &= \ln c(\mathcal{D}) + \sum_{i=1}^n \ln \left\{ \frac{1}{\sigma_i} \varphi \left(\frac{Z_i - \theta_i}{\sigma_i} \right) \right\} + o_{P_0}(1), \quad \text{where} \\
c(\mathcal{D}) &:= p(\mathcal{D} \mid 0, \dots, 0) \cdot \prod_{i=1}^n \sqrt{2\pi \sigma_i^2} \exp \left\{\frac{ Z_i^2}{2\sigma_i^2} \right\}.
\end{align*}
Thus, up to a proportionality constant that does not depend on \( (\theta_1, \dots, \theta_n) \), the log-likelihood of the data \( p(\mathcal{D} \mid \theta_1, \dots, \theta_n) \) is well-approximated by the working log-likelihood
\[
\sum_{i=1}^n \ln \left\{\frac{1}{\sigma_i} \varphi \left(\frac{Z_i - \theta_i}{\sigma_i} \right) \right\}.
\]
This indicates that our $g$-modeling framework remains fundamentally unchanged in the parametric setting, provided \( Z_i \) is redefined as the average score in each experiment \( i \).
\@startsection{subsection}{2}{\z@}
{0pt}
{0.4ex}
{\normalfont\normalsize\bfseries}{Alternative loss functions}
While our theoretical results have focused on the MSE loss function, our methodology easily generalizes to other losses. For instance, \cite{chen2022empirical} describes a loss function geared towards selecting $\theta_i$ that are larger than 0:
$$
L(\bm{\delta}, \bm{\theta}) = \frac{-1}{n} \sum_{i=1}^n \delta_i \theta_i.
$$
Here $\bm{\delta} := (\delta_1, \dots, \delta_n)$ and each $\delta_i \in \{0,1\}$ is a binary action, with $\delta_i = 1$ indicating selection. For a given prior $G$, the optimal Bayes decision is $\delta_i = \mathbb{I}\{\theta_{i,G} \ge 0\}$, where $\theta_{i,G}$ is the posterior mean of $\theta_i$ under that prior. The Empirical Bayes (EB) analogue of this strategy would be to replace $G$ with its EB estimate $\hat{G}$; in terms of our previous notation, this corresponds to using $\hat{\delta}_i^{\textrm{EB}}$ or $\hat{\delta}_i^{\textrm{NPMLE}}$ in place of $\theta_{i,G}$, depending on which $g$-modeling procedure is employed.
The selection problem could be generalized to top-$m$ selection, where the aim is to determine the best $m$ values of $\theta_i$. Again, the optimal Bayes decision, given a prior $G$, would be to select the top $m$ values based on the posterior means $\theta_{i,G}$. The EB strategy again simply replaces $G$ with $\hat{G}$.
For threshold selection, if $\widehat\theta_i$ and
$\theta_{i,G_0}$ denote the estimated and oracle posterior means,
the conditional regret of the plug-in decision $\mathbb{I}\{\hat{\theta}_i \ge 0\}$ is bounded by
\[
\frac1n\sum_i
|\theta_{i,G_0}|\,
1\{\operatorname{sign}(\widehat\theta_i)
\ne\operatorname{sign}(\theta_{i,G_0})\}
\leq
\frac1n\sum_i|\widehat\theta_i-\theta_{i,G_0}|.
\]
Taking expectations and invoking the Cauchy–Schwarz inequality shows that this quantity is upper bounded by the square root of the associated MSE regret. For top-$m$ selection, the reasoning in \cite{chen2022empirical} leads to an analogous result. Consequently, our regret guarantees for MSE loss, established in earlier sections, extend to these alternative loss functions as well.
\section{Simulations}
We now present simulation results evaluating the performance of our methods for treatment effect estimation using adaptively generated data from common algorithms. For simplicity, we focus on one-armed bandit experiments employing Thompson Sampling (TS) and the Upper Confidence Bound (UCB) algorithm, and we compare various shrinkage estimators. The maximum sample size of each experiment is set to be 50. Each experiment yields a pair of summary statistics \((Z_i, \sigma_i)\), corresponding to the sample mean and sample standard deviation.
In the first data-generating process (DGP), the prior distribution of \(\theta_1, \dots, \theta_n\) is assumed to be Gaussian: \(\theta_i \sim N(m_0, s_0^2)\) with \(m_0 = 0\) and \(s_0^2 = 1/4\). The outcome distribution for each experiment follows a Gaussian with mean \(\theta_i\) and unit variance. This setup favors linear shrinkage estimators. Table \ref{tab:gaussianMSE} reports the mean squared error (MSE) for estimating \(\theta_i\). The Oracle assumes perfect knowledge of \(G_0\) when constructing the Bayes rule, while NPMLE estimates \(G_0\) non-parametrically. Several intermediate linear shrinkage estimators are also considered:
\begin{itemize}
\item L-Marginal estimates the prior parameters \((m_0, s_0^2)\) using marginal moments of \((Z_i, \sigma_i)\). The linear shrinkage estimator takes the form:
\[
\hat{m}_0 + (Z_i - \hat{m}_0) \frac{\hat{s}_0^2}{\sigma_i^2 + \hat{s}_0^2}.
\]
This corresponds to the standard empirical Bayes approach (e.g., \cite{kane2008does}, \cite{chetty2014measuringa}), where \((m_0, s_0^2)\) is estimated from the marginal moments of \((Z_i, \sigma_i)\).\footnote{A commonly used estimator is \(\hat{m}_0 = \frac{1}{n} \sum_i Z_i\) and \(\hat{s}_0^2 = \widehat{\text{Var}}(Z_i) - \hat{\mathbb{E}}(\sigma_i^2)\), based on the identities \(\mathbb{E}[Z_i] = \mathbb{E}[\theta_i] = m_0\) and \(\text{Var}[Z_i] = \mathbb{E}[\text{Var}[Z_i | \theta_i]] + \text{Var}[\mathbb{E}[Z_i | \theta_i]]\).}
\item L-Posterior follows the same linear shrinkage formula but estimates \((m_0, s_0^2)\) via posterior moment matching.
\item L-LOO is similar to L-Posterior, except that prior parameter estimation is performed using a leave-one-out approach, as discussed earlier.
\end{itemize}
For comparison, we also report the performance of the MLE estimator for \(\theta_i\) as a benchmark.
Under the Gaussian prior, both L-Posterior and L-LOO closely track the performance of the Oracle estimator. NPMLE is also highly competitive, particularly for larger sample sizes. While its performance lags behind for small samples (\(n = 100\)), its performance catches up with the L-posterior as $n$ increases. In contrast, L-Marginal performs worse due to a mismatch between the working likelihood and the true marginal density of \((Z_i, \sigma_i)\). As a result, marginal moments of \((Z_i, \sigma_i)\) fail to provide consistent estimates of the prior parameters.
Table \ref{tab:momentGaussian} further illustrates this point by reporting the mean squared error (MSE) for estimating \((m_0, s_0^2)\) under the various methods. The columns labeled ``Marginal'' correspond to estimates based on marginal moments, which exhibit high MSE due to significant bias. In contrast, the ``Posterior'' approach, which estimates \((m_0, s_0^2)\) via posterior moment matching, and NPMLE, which reports the mean and variance of \(\hat{G}_n\), both achieve substantially lower MSE across all algorithms.
\begin{table}
\centering
\small
\caption{One armed bandit, Gaussian prior} \label{tab:gaussianMSE}
\begin{tabular}{rrrrrrrr}
\hline
n & Oracle & NPMLE & L-marginal & L-loo & L-posterior & MLE \\
\hline
\multicolumn{7}{c}{Thompson Sampling }\\
$100$ & 0.0560 & 0.0686 & 0.0646 & 0.0574 & 0.0572 & 0.1744 \\
$500$ & 0.0561 & 0.0605 & 0.0641 & 0.0564 & 0.0564 & 0.1790 \\
$1000$ & 0.0562 & 0.0586 & 0.0637 & 0.0564 & 0.0564 & 0.1770 \\
$5000$ & 0.0562 & 0.0570 & 0.0635 & 0.0563 & 0.0563 & 0.1775 \\
\hline
\multicolumn{7}{c}{UCB algorithm }\\
$100$ & 0.0605 & 0.0717 & 0.0720 & 0.0623 & 0.0622 & 0.1959 \\
$500$ & 0.0605 & 0.0647 & 0.0715 & 0.0608 & 0.0608 & 0.1976 \\
$1000$ & 0.0605 & 0.0628 & 0.0716 & 0.0606 & 0.0606 & 0.1986 \\
$5000$ & 0.0607 & 0.0613 & 0.0716 & 0.0607 & 0.0607 & 0.1983 \\
\hline
\end{tabular}\\
\begin{flushleft}
\setstretch{1.2}
\scriptsize{Notes: One-armed bandit experiments. $\theta_i \sim N(0, \tfrac{1}{4})$. Mean squared error for $\theta_i$ estimation. Results are based on 500 simulation repetitions.}
\end{flushleft}
\end{table}
\begin{table}
\centering
\small
\caption{One armed bandit, Gaussian prior} \label{tab:momentGaussian}
\begin{tabular}{rrrr|rrr}
\hline
n & Marginal & Posterior & NPMLE & Marginal & Posterior & NPMLE \\
\hline
\multicolumn{7}{c}{MSE of estimator for Prior Mean}\\
\hline
\multicolumn{1}{c}{}& \multicolumn{3}{c}{Thompson Sampling }& \multicolumn{3}{c}{UCB algorithm}\\
\hline
100 & 0.0222 & 0.0038 & 0.0039 & 0.0336 & 0.0036 & 0.0038 \\
500 & 0.0195 & 0.0007 & 0.0007 & 0.0297 & 0.0007 & 0.0008 \\
1000& 0.0185 & 0.0004 & 0.0004 & 0.0298 & 0.0004 & 0.0004 \\
5000 & 0.0181 & 0.0001 & 0.0001 & 0.0295 & 0.0001 & 0.0001 \\
\hline \hline
\multicolumn{7}{c}{MSE of estimator for Prior Variance}\\
\hline
\multicolumn{1}{c}{}& \multicolumn{3}{c}{Thompson Sampling }& \multicolumn{3}{c}{UCB algorithm}\\
\hline
100 & 0.0218 & 0.0023 & 0.0048 & 0.0203 & 0.0027 & 0.0048 \\
500 & 0.0190 & 0.0005 & 0.0009 & 0.0185 & 0.0005 & 0.0010 \\
1000 & 0.0179 & 0.0002 & 0.0005 & 0.0180 & 0.0003 & 0.0005 \\
5000 & 0.0176 & 0.0000 & 0.0001 & 0.0176 & 0.0001 & 0.0001 \\
\hline\hline
\end{tabular}\\
\begin{flushleft}
\setstretch{1.2}
\justify \scriptsize{Notes: One-arm bandit experiments with Thompson Sampling or UCB algorithm. Mean-squared error for the estimation of prior mean and variances. The true prior is $G_0 = N(0, 1/4)$. `Marginal' corresponds to learning the moments of $G_0$ based on $(Z_i,\sigma_i)$, `Posterior' corresponds to learning the moments of $G_0$ through posterior moment matching, `NPMLE' corresponds to the moments implied by NPMLE estimate. Results are based on 500 simulation repetitions.}
\end{flushleft}
\end{table}
For the second DGP, we consider a discrete $G_0$ with two atoms. This is a setting that naturally favors NPMLE. Methods assuming a Gaussian prior in this case fail to recover the correct posterior distribution and therefore cannot replicate the optimal Bayes rule. We should emphasize that Stein's result - that MLE is inadimissible - no longer holds under adaptive sampling. Nevertheless, we observe that the MLE still exhibits higher MSE compared to linear shrinkage methods. Among all the estimators, NPMLE clearly demonstrates the best performance and also tracks the Oracle very closely. It mimics the optimal Bayes rule and is able to learn $G_0$ with remarkable accuracy even in moderate sample sizes. The detailed results are reported in Tables \ref{tab:twoGroupMSE} and \ref{tab:momentTwo_point}.
\begin{table}
\centering
\small
\caption{One armed bandit, two-point prior} \label{tab:twoGroupMSE}
\begin{tabular}{rrrrrrr}
\hline
& Oracle & NPMLE & L-marginal & L-loo & L-posterior & MLE \\
\hline
\multicolumn{7}{c}{Thompson Sampling }\\
$n = 100$ & 0.0000 & 0.0139 & 0.1267 & 0.1227 & 0.1243 & 0.2096 \\
$n = 500$& 0.0000 & 0.0029 & 0.1222 & 0.1196 & 0.1199 & 0.2037 \\
$n = 1000$ & 0.0000 & 0.0016 & 0.1227 & 0.1202 & 0.1204 & 0.2047 \\
$n = 5000$ & 0.0000 & 0.0004 & 0.1228 & 0.1205 & 0.1205 & 0.2047 \\
\hline
\multicolumn{7}{c}{UCB algorithm }\\
$n = 100$ & 0.0000 & 0.0104 & 0.1039 & 0.1013 & 0.1023 & 0.1654 \\
$n = 500$ & 0.0000 & 0.0032 & 0.1049 & 0.1030 & 0.1033 & 0.1683 \\
$n = 1000$ & 0.0000 & 0.0018 & 0.1050 & 0.1032 & 0.1033 & 0.1684 \\
$n = 5000$ & 0.0000 & 0.0004 & 0.1046 & 0.1029 & 0.1029 & 0.1682 \\
\hline
\end{tabular}
\begin{flushleft}
\setstretch{1.2}
\justify \scriptsize{Notes: One-arm bandit experiments. $\theta_i \sim \frac{1}{2}\delta_{-1} + \frac{1}{2} \delta_3$. Mean squared error for $\theta_i$ estimation. Results are based on 500 simulation repetitions.}
\end{flushleft}
\end{table}
\begin{table}
\centering
\small
\caption{One armed bandit, two-point prior} \label{tab:momentTwo_point}
\begin{tabular}{rrrr|rrr}
\hline
n & Marginal & Posterior & NPMLE & Marginal & Posterior & NPMLE \\
\hline
\multicolumn{7}{c}{MSE of estimator for Prior Mean}\\
\hline
\multicolumn{1}{c}{}& \multicolumn{3}{c}{Thompson Sampling }& \multicolumn{3}{c}{UCB algorithm}\\
\hline
100 & 0.0662 & 0.0529 & 0.0451 & 0.0582 & 0.0476 & 0.0431 \\
500 & 0.0205 & 0.0102 & 0.0077 & 0.0212 & 0.0100 & 0.0083 \\
1000 & 0.0169 & 0.0064 & 0.0042 & 0.0167 & 0.0054 & 0.0041 \\
5000 & 0.0134 & 0.0029 & 0.0008 & 0.0122 & 0.0017 & 0.0009 \\
\hline \hline
\multicolumn{7}{c}{MSE of estimator for Prior Variance}\\
\hline
\multicolumn{1}{c}{}& \multicolumn{3}{c}{Thompson Sampling }& \multicolumn{3}{c}{UCB algorithm}\\
\hline
100 & 0.2578 & 0.1781 & 0.0366 & 0.1567 & 0.1181 & 0.0348 \\
500 & 0.2395 & 0.1579 & 0.0056 & 0.1614 & 0.1146 & 0.0061 \\
1000 & 0.2426 & 0.1600 & 0.0028 & 0.1617 & 0.1144 & 0.0032 \\
5000 & 0.2426 & 0.1601 & 0.0005 & 0.1642 & 0.1149 & 0.0006 \\
\hline
\end{tabular}
\begin{flushleft}
\setstretch{1.2}
\justify \scriptsize{Notes: One-armed bandit with Thompson Sampling or UCB algorithm. Mean-squared error for the estimation of prior mean and variances, when true prior is $G_0 = \frac{1}{2} \delta_{-1} + \frac{1}{2}\delta_3$. `Marginal' corresponds to learning the moments of $G_0$ based on $(Z_i,\sigma_i)$, `Posterior' corresponds to learning the moments of $G_0$ through posterior moment matching assuming normal prior, `NPMLE' corresponds to the moments implied by NPMLE. Results are based on 500 simulation repetitions. }
\end{flushleft}
\end{table}
\section{Empirical Illustration \protect \label{sec:Empirical illustration}}
In this section, we utilize the $g$-modeling framework to estimate the prior distribution and treatment effects using data from the ASOS digital experiments dataset. This dataset comprises results from $n = 61$ A/B tests conducted under adaptive stopping rules. In each experiment, the treatment and control groups were sampled in equal proportions. When analyzing these experiments using NPMLE, we make the implicit assumption that the experiments are independent of each other, so that Assumption 1 in Section \ref{subsec:regret_consistency_NPMLE} holds.
The relevant terminology for this analysis is defined in Section \ref{subsec:Motivating example}. In this context,
$$
Z_i = \sqrt{N} \left(\bar{Y}^{(1)}_i - \bar{Y}^{(0)}_i \right)
$$
represents the difference in sample means between the treatment and control groups, scaled by the mean sample size $\sqrt{N}$, while $\tau_i$ represents the number of times each arm was sampled divided by $N$. Also,
$$\sigma_i^2 = \frac{\textrm{Var}[Y_{j,i}^{(1)}] + \textrm{Var}[Y_{j,i}^{(0)}]}{\tau_i}.$$
A key distinction from the discussion in Section \ref{subsec:Motivating example} is that the outcomes \( Y_{j,i}^{(1)} \) and \( Y_{j,i}^{(0)} \) are binary, following a Bernoulli distribution, rather than the normal distribution assumed in that section. Nevertheless, as demonstrated in Section \ref{subsec:local asymptotics}, the working likelihood $\prod_i \frac{1}{\sigma_i} \varphi\left(\frac{Z_i - \theta_i}{\sigma_i}\right)$
remains a valid approximation to the true likelihood, up to a proportionality constant that does not depend on the parameters \( \theta_1, \dots, \theta_n\).
Our objective is to estimate the scaled experiment-specific treatment effects
\[
\theta_i = \sqrt{N} \cdot \mathbb{E}[Y_{j,i}^{(1)} - Y_{j,i}^{(0)}],
\]
along with the common prior distribution, $G_0$, from which these effects are generated. For this purpose, we implement the $g$-modeling procedures outlined in Section \ref{sec:Empirical Bayes}.
These procedures rely on knowledge of $\sigma_i^2$, which is a function of $\textrm{Var}[Y_{j,i}^{(1)}]$ and $\textrm{Var}[Y_{j,i}^{(0)}]$, respectively. Although these variances are unknown, they can be approximated using standard sample variance estimators. Notably, these estimators remain consistent in adaptive experimental settings, even though their rates of convergence are slower compared to classical experiments.
Figure \ref{fig: histogram} plots the histogram of $Z_i$, the raw estimates of scaled estimated treatment effects. Most of the treatment effect estimates are fairly small, even after scaling, although we do see 6 experiments with noticeably larger positive treatment effects. The stopping time across the 61 experiments also is heterogeneous and its relationship with $Z_i$ is plotted in the left panel of Figure \ref{fig: close}. We observe a slight positive correlation. Experiments with a largest $Z_i$ tends to have larger $\sigma_i$, which is evidence of earlier stopping. Using our proposed method, we estimate $G_0$ using the NPMLE with the prescribed working model for $Z_i$ and the estimator is presented in the right panel of Figure \ref{fig: close}. Interestingly, the four experiments with the largest $Z_i$ are all given their own mass points. The rest of the mass points cluster around zero.
The estimator of $\hat G_n$ drives the shrinkage pattern for the nonparametric EB estimator. Figure \ref{fig:shrinkage} shows the amount of shrinkage instructed by the nonparametric along with the linear empirical Bayes method. The two methods agree on experiments with more extreme $Z_i$'s. The four experiments associated with the largest $Z_i's$ are not shrunk much by either method. Two experiments with $Z_i$'s close to 5 are shrunk aggressively, more so under the nonparametric method to the extent that their relative rankings are also changed. A close inspection shows these two experiments are associated with extremely large variances. On the other hand, experiments with their effects close to zero are adjusted more towards zero under nonparametric EB, while the linear method aligns them closely on the 45 degree line, hence leaving the MLE unrevised.
\begin{figure}
\includegraphics[scale = 0.35]{histogram.pdf}
\caption{Histogram of $Z_i$ across 61 A/B tests}
\label{fig: histogram}
\end{figure}
\begin{figure}
\includegraphics[scale = 0.45]{close.pdf}
\caption{The left panel plots the $\log_{10}(\sigma_i)$ on the x-axis and $Z_i$ on the y-axis. The right panel plots the NPMLE $\hat G_n$ of $G_0$ using the working likelihood model for $Z_i$.}
\label{fig: close}
\end{figure}
\begin{figure}
\includegraphics[scale = 0.45]{shrinkage}
\caption{The left panel compared g-modeling based nonparametric empirical Bayes estimator with the MLE estimator for $\theta_i$ (the latter is the 45 degree line). The right panel compares the linear shrinkage estimator, using either marginal moments matching or posterior moments matching, with the MLE.}
\label{fig:shrinkage}
\end{figure}
\section{Conclusion}
Adaptive data collection need not make Empirical Bayes practice
algorithm-specific. As long as the sampling rules and stopping times depend on
the data only through the observed history, the likelihood factors into the
familiar Gaussian working likelihood and a term that is parameter free.
Standard $g$-modeling therefore recovers the common distribution
of treatment effects, and with it, the posterior estimates, without any
knowledge of the sampling or stopping rules. Our moment-matching
interpretation explains this robustness and, at the same time, why naive
$f$-modeling fails when adaptivity distorts the marginal distribution of the
observed sample means.
We establish identification and regret guarantees for both parametric and
nonparametric $g$-modeling, and our simulations and the ASOS application show
that the method remains effective when the sampling algorithms
are complex or unknown. Two key assumptions still remain:
that the experiments are comparable enough to
share a common prior and, for the theoretical analysis of NPMLE, that they are independent.
Relaxing these requirements is a natural next step. More broadly, we envision that our results
can turn collections
of adaptive experiments into a resource not only for improving past estimates
but also for designing better experimentation strategies for future experiments.
\bibliographystyle{IEEEtranSN}
\bibliography{EB_adaptive_experiments}