Amortized simulation-based frequentist inference for tractable and intractable likelihoods

A classical hypothesis test is designed to have power in distinguishing between the null and the alternative hypothesis. Consider the α-level hypothesis $\theta = \theta_0$, as in equation (1). This hypothesis is to be rejected if $\mathbb{P}(\lambda \gt \lambda_D \mid \theta_0) \lt \alpha$. The corollary is that θ₀ is not to be rejected if $\mathbb{P}(\lambda \gt \lambda_D \mid \theta_0) \unicode{x2A7E} \alpha$, that is, if $\mathbb{P}(\lambda \unicode{x2A7D} \lambda_D \mid \theta_0) \unicode{x2A7D} 1 - \alpha \equiv \tau$. The set of points θ₀ that have not been rejected at level α, and therefore remain as potentially viable hypotheses for the true value of θ, is by definition a confidence set R(D) at 100τ% confidence level (CL). The boundary of the confidence set R(D) is determined by the equation

$$\begin{align} \mathbb{P}\left(\lambda \unicode{x2A7D} \lambda_D \mid \theta_0\right) & = \tau. \end{align} \tag{ 5 }$$

By construction, the coverage probability of such sets is

$$\begin{equation} \mathbb{P}\left(\theta \in R\left({\cal D}\right) \mid \theta\right) \unicode{x2A7E} \tau. \end{equation} \tag{ 6 }$$

Note that $R(\mathcal{D})$ is a random set. Under repeated observations, the frequentist principle requires that the relative frequency with which these sets include the true value of θ never falls below the stated CL $\tau = 1 - \alpha$ regardless of the true value of θ. To the degree that the probabilities $\mathbb{P}(\lambda \gt \lambda_D \mid \theta)$ and $\mathbb{P}(\lambda \unicode{x2A7D} \lambda_D \mid \theta)$ are accurately modeled, the LF2I and ALFFI approaches yield random sets that satisfy the frequentist principle.

The LF2I and ALFFI approaches are applicable to any test statistic $\lambda(\mathcal{D}; \theta)$ that is monotonic in the following sense: the test statistic is constructed so that a particular direction in its 1-dimensional space corresponds to hypotheses that are increasingly disfavored. In ALFFI, we consider test statistics for which large values correspond to disfavored hypotheses, or equivalently, small values of the $\textit{p}{\textrm{-}}\textrm{value} = \mathbb{P}(\lambda \gt \lambda_D | \theta = \theta_0)$, where $\lambda_D \equiv \lambda(\mathcal{D} = D; \theta_0)$ is the observed value of the test statistic for the specified hypothesis. ALFFI (see appendix A) approximates the probability

$$\begin{align} \mathbb{P}\left( \lambda \unicode{x2A7D} \lambda_D \mid \theta\right) & = \int^{\lambda_D} \textrm{d} Y \int \textrm{d}{\cal D} \, \delta\left(Y - \lambda\left({\cal D}, \theta\right)\right) \, p\left({\cal D} \mid \theta\right) , \end{align} \tag{ 7 }$$

where $\delta (.)$ is the Dirac delta function [9] and $p({\cal D} | \theta)$ is a statistical model, which may or may not be tractable. However, as in the LF2I approach [3], it is assumed that one has access to a large collection $\mathbb{S}$ of simulated pairs $({\cal D}, \theta) \in \mathbb{S}$ where for every point θ sampled from any convenient prior π_θ a single instance of a data set ${\cal D}$ is simulated with the same characteristics, including sample size, as the real data D. In contrast to LF2I, in ALFFI a second collection of data sets $\mathbb{S}^{\prime}$ is created from $\mathbb{S}$ by randomizing the order of the data sets $\{{\cal D}\}$ to form new pairs $({\cal D}^{\prime}, \theta) \in \mathbb{S}^{\prime}$ in which the parameters and the data sets ${\cal D}^{\prime}$ are statistically independent. The data sets in $\mathbb{S}^{\prime}$ serve as instances of 'observed' data sets.

For every parameter θ in $\mathbb{S}$, we compute two values of the test statistic, namely, $\lambda({\cal D}, \theta)$ and $\lambda_D = \lambda({\cal D}^{\prime}, \theta)$, as well as the discrete variable Z which is unity if $\lambda({\cal D}, \theta) \unicode{x2A7D} \lambda({\cal D}^{\prime}, \theta)$ and zero otherwise. This procedure results in a large collection of triplets of size B, $\mathcal{T} = \{(Z_i, \lambda_{D,i}, \theta_i )\}_{i = 1}^{B}$, which constitute the training data. In LF2I, the observed test statistic under the null that $\theta = \theta_0$—the second component of the triplet—is computed using a fixed data set ${\cal D}^{\prime} = D$, namely, the one actually observed, while ALFFI uses the data sets ${\cal D}^{\prime} \in \mathbb{S}^{\prime}$.

The cdf $\mathbb{P}( \lambda \unicode{x2A7D} \lambda_D | \theta)$ is the expectation value $\mathbb{E}(Z | \lambda_D, \theta)$ of the discrete variable Z, a fact that suggests a straightforward way to approximate the cdf for a fixed data set D: histogram the parameter points θ, thereby yielding the histogram $\mathbb{H}_1$, and using the same bins as $\mathbb{H}_1$ histogram the parameter points again, but this time weighted by Z, yielding the histogram $\mathbb{H}_Z$. The ratio $\mathbb{H}_Z / \mathbb{H}_1$ provides a piece-wise-constant approximation of $\mathbb{E}(Z |\lambda_D, \theta)$.

Following LF2I, a smooth approximation of $\mathbb{E}(Z |\lambda_D, \theta)$ is created with a deep neural network (DNN) trained (that is, fitted) by minimizing the empirical risk

$$\begin{align} E\left(\boldsymbol{\omega}\right) & = \frac{1}{N} \sum_{i = 1}^N \mathcal{L}\left(t_i, f_i\right), \end{align} \tag{ 8 }$$

where $\mathcal{L}(t, f)$ is a loss function, $f(\textbf{x}; \boldsymbol{\omega})$ is a DNN with inputs x and free parameters ω , and t denotes known targets. The empirical risk or average loss, equation (8), is a Monte Carlo approximation of the risk functional

$$\begin{align} E\left[\,f\,\right] & = \int \int \mathcal{L}\left(t, f\right) \, p\left(\textbf{x}, t\right) \, \textrm{d}\textbf{x} \, \textrm{d}t, \end{align} \tag{ 9 }$$

where $p(\textbf{x}, t) = p(t | \textbf{x}) \, p(\textbf{x})$ is the (typically unknown) probability distribution of the data $\textbf{x}, t$. From the calculus of variations, the function f that minimizes equation (9) is the solution of

$$\begin{align} \int \frac{\partial \mathcal{L}}{\partial f} \, p\left(t | \textbf{x}\right) \, \textrm{d}t & = 0, \end{align} \tag{ 10 }$$

assuming that $p(\textbf{x}) \gt 0~\forall \, \textbf{x}$. In the examples below, and following LF2I, we use the quadratic loss $\mathcal{L}(t, f) = (t - f)^2$, which, using equation (10), leads to the well-known result [10, 11]

$$\begin{align} f\left(\textbf{x}; \boldsymbol{\omega}^*\right) & = \int t \, p\left(t | \textbf{x}\right) \, \textrm{d}t \equiv \mathbb{E}_t\left[t \mid \textbf{x}\right], \end{align} \tag{ 11 }$$

where $\boldsymbol{\omega}^*$ are the best-fit parameters of the neural network model. Setting $\textbf{x} = \{\lambda_D, \theta$} and the targets t = Z in equation (11) yields

$$\begin{align} f\left(\lambda_D, \theta; \boldsymbol{\omega}^*\right) & \approx \mathbb{E}\left[ Z \mid \lambda_D, \theta \right],\nonumber\\ & = \mathbb{P}\left( \lambda \unicode{x2A7D} \lambda_D \mid \theta\right), \end{align} \tag{ 12 }$$

that is, it yields the quantity that we wish to approximate.

One of the key virtues of the LF2I and ALFFI approaches, as is evident in the result in equation (12), is that the neural network f is conditioned on θ, which implies that it is independent of the prior π_θ from which the parameter points are sampled. The form of the prior affects only the accuracy of the approximation: the accuracy of the approximation will be greatest where the density of the prior is greatest.

In the late 1990 s, fits of cosmological models to Type 1a supernova data led to the conclusion that the expansion of the Universe is accelerating [12, 13]. The fits then, as now, were performed using tractable likelihoods, typically a multivariate normal. In this example, we fit a cosmological model to the Union 2.1 data compilation of the Supernova Cosmology Project [14] via maximum likelihood and also with ALFFI. The Union 2.1 data set comprises measured distance moduli, $x \pm \sigma$, and redshifts, z, for 580 Type 1a supernovae. Given the size of the data sample, it is expected that accurate confidence sets for the cosmological parameters can be constructed using standard methods such as maximum likelihood. ALFFI is therefore not needed for this problem; it is simply used to showcase the algorithm.

Our cosmological model is defined by the equation of state

$$\begin{align} {\cal P} & = - n a^n \Omega / 3, \end{align} \tag{ 13 }$$

where n is a free parameter and a(t), $\Omega(a)$, and ${\cal P}$ are the dimensionless universal scale factor, the dimensionless energy density, and the dimensionless pressure, respectively, and t is the time since the Big Bang. Together with the Hubble constant ${\cal H}_0$ this is a 2-parameter problem. (Details of the model are given in appendix B). If the small correlations between the 580 Typa 1a data points are neglected, the likelihood for these data is a diagonal multivariate normal. Maximizing this likelihood with respect to its parameters is equivalent to minimizing the function

$$\begin{equation} \chi^2 = \sum_{i = 1}^N \left( \frac{x_i - \mu\left(z_i, \theta\right)}{\sigma_i}\right)^2 , \end{equation} \tag{ 14 }$$

where $\mu(z, \theta)$, the distance modulus[15]—an astronomical measure of distance, is given by

$$\begin{equation} \mu\left(z, \theta\right) = 5 \log_{10}\left[\left(1 + z\right) \, \sin\left(\sqrt{-\Omega_K} \, u\left(z, \theta\right)\right)\, / \, \sqrt{-\Omega_K}\right] + 5 \log_{10}\left(c \, / \, {\cal H}_0 / 10^{-5} \textrm{Mpc}\right). \end{equation} \tag{ 15 }$$

The quantity c is the speed of light in vacuum in km/s, $\Omega_K$ is the curvature parameter, which we set to zero, and

$$\begin{equation} u\left(z, \theta\right) = \int_{1/\left(1+z\right)}^{1} \frac{da}{a^2\sqrt{\Omega\left(a\right)}},\nonumber\\ = 2^{\frac{1}{2 n}} \left[\gamma\left(\frac{1}{2 n}, \frac{1}{2}\right) - \gamma\left(\frac{1}{2 n}, \frac{\left(1 + z\right)^{-n}}{2}\right)\right] \sqrt{e} / \, n , \end{equation}$$

is a dimensionless function. When the model is fitted to the Union 2.1 data set by minimizing equation (14) an excellent fit is obtained, as shown in figure 1.

Zoom In Zoom Out Reset image size

We now apply ALFFI to the same problem using the test statistic

$$\begin{equation} \lambda\left(D, \theta\right) = \sqrt{\frac{\chi^2}{N} }, \end{equation} \tag{ 16 }$$

where χ² is defined in equation (14), and satisfies the requirement by ALFFI that large values of the test statistic cast doubt on the null hypothesis. The form of the test statistic is chosen so that it is ${\cal O}(1)$ as it is an input to the neural network. The boundary of the associated 100τ% confidence set is given by equation (5), which requires a good approximation to the cdf $\mathbb{P}(\lambda \unicode{x2A7D} \lambda_D | \theta)$. The latter is approximated in two ways: with histograms and with a deep neural network (DNN) as described in section 2.3. The 2D histograms have 10 bins in both dimensions n and ${\cal H}_0$ with the parameter n scaled down by a factor 10 and ${\cal H}_0$ scaled down by a factor 100, so that both input parameters are of $\mathcal{O}(1)$.

The DNN is 1781-parameter fully-connected feed-forward neural network, with 3 input features, $\boldsymbol{x} = \{\lambda_D, n, {\cal H}_0 \}$, 5 hidden layers with 20 nodes each, and a single output. We use a ReLU [16] activation function, the output node is a sigmoid that constrains the output to lie within the unit interval, and the DNN is trained as described in section 2 using PyTorch [17].

We use a batch size $K = 50 \ll N$ randomly sampled from $N = 250\,000$ data sets of 580 simulated distance moduli per data set, sampled at the same redshifts as the observed data. The Adam [18] optimizer with a fixed learning rate of 10⁻³ is used to train the network. As the training proceeds, the network with the smallest average loss is saved. The average loss is computed using a validation data set of size 5000 that is not used by the optimizer. We employ the following early stopping [19] criterion: the training stops if after 50 000 steps no network is found with a smaller validation loss than the saved network. The best network as defined by this training protocol is found after about 100 000 iterations. If one defines an epoch to be $N / K$ iterations, which for this example is 5000, then $100\,000$ iterations corresponds to 20 epochs of training. The simulated data sets of 580 distance moduli are sampled from 580 independent normal distributions using the standard deviations taken from the observed data.

Approximating $\mathbb{P}(\lambda \unicode{x2A7D} \lambda_D | \theta)$ using histograms and using ALFFI leads to the confidence sets shown in figure 2.

Zoom In Zoom Out Reset image size

The best-fit values of the cosmological parameters $\theta = \{n, {\cal H}_0 \}$ are taken to be the location of the minimum of $\mathbb{P}(\lambda \unicode{x2A7D} \lambda_D | \theta)$, which as indicated in figure 2 agrees with the values obtained from the likelihood fit.

In the LF2I approach, the coverage over the parameter space is checked by modeling the coverage probability with another neural network as a function of the parameters of the theoretical model. There are pros and cons to that approach. It is certainly convenient to have a functional approximation of the coverage probability as a function of the parameter space point because one can then estimate the coverage at any given point, not only at points for which there are sufficient simulated data. Unfortunately, however, as is true of most machine learning models, a reliable estimate of the accuracy of the trained machine learning model is not available. In ALFFI, the coverage is checked explicitly by direct enumeration at all points for which there are sufficient data. Since the problem consists of counting how often a particular statement is true, the problem is binomial; therefore, a reliable estimate of the accuracy of the coverage calculation is easy to compute. On the other hand, the coverage is available only at the parameter points for which there are sufficient simulated data. In this example, for a given parameter point, $T = 4,000$ sets of 580 Type1a supernovae data are simulated. For each data set, a test statistic, λ_D , is computed. If $\mathbb{P}(\lambda \unicode{x2A7D} \lambda_D | \theta) \unicode{x2A7D} \tau$ then, by definition, θ lies within the confidence set associated with λ_D . If S is the number of times this statement is true over the collection of T simulated data sets, then the coverage probability is $p \pm \sqrt{p (1 - p) / T}$, where $p = S / T$. For continuous probability distributions and exact confidence sets, the coverage probability p is exactly equal to the CL τ. Therefore, if the confidence sets produced by ALFFI are accurate then we should find $p \approx \tau$ given that we are making an approximation. This calculation is performed at 500 randomly sampled points within the 95% CL set associated with the observed Type 1a data, as shown in the left panel of figure 3.

Zoom In Zoom Out Reset image size

The right panel shows the coverage probabilities calculated for the 500 randomly sampled parameter points, shown in the left panel, compared with the desired $1-\alpha = \tau$ CLs, depicted by the solid horizontal lines. Since the coverage probabilities are greater or equal to the CLs, we conclude that the coverage of the confidence sets computed using ALFFI satisfy the coverage condition in equation  (6).

The cosmological example served to illustrate the ALFFI algorithm, but, as noted, ALFFI is not really needed because the cosmological data are numerous, the likelihood function is tractable, and parameter estimation via maximum likelihood works well. Our second example addresses the signal/background problem in high-energy physics (see, for example, [20]), which in astronomy is referred to as the On/Off problem [8]. We choose an example in which the likelihood is tractable but the data are sparse and, consequently, asymptotic methods may not be reliable [7]. We demonstrate that the ALFFI method yields valid confidence sets even when the regularity conditions that underpin Wilks' theorem [21] and its variants [6] are violated.

The signal/background problem in high-energy physics and astronomy is as follows. An observation is made, for given period of time, which consists of counting N events: typically, photons in astronomy and particle collisions in high-energy physics. The count is potentially a sum of counts from signal and background sources. A second independent observation is made for the same duration (in the simplest case) where by design the background has the same characteristics as in the first observation but no signal is present. The second observation yields a count M. It is generally assumed that the likelihood function for the data $D = \{N, M \}$ is the product of two Poisson distributions,

$$\begin{align} \mathcal{L} \left(D ; \mu, \nu\right) & = \frac{\left(\mu + \nu\right)^N \exp\left(-\left(\mu + \nu\right)\right)}{N!} \,\frac{\nu^M\exp\left(-\nu\right)}{M!}, \end{align} \tag{ 17 }$$

where µ and ν are the mean signal and background counts, respectively. In the signal/background problem the parameter of interest is µ, while ν is a nuisance parameter. We shall comment on how one might deal with such parameters in the discussion.

Our specific example is from the first experiment to search for neutron-antineutron oscillations using free neutrons [22], which took place in the 1980 s at the Institut Laue-Langevin (ILL) in Grenoble, France. Neutron-antineutron ($n\bar{n}$) oscillations are predicted by many proposed theories of physics beyond the Standard Model of particle physics [23]. For our purposes it suffices to note that if $n\bar{n}$ oscillations can occur then a pure neutron state, when observed at time $t \ll$ than the mean neutron lifetime, will be observed with probability

$$\begin{align} P_t & = \left(\frac{\epsilon^2}{\epsilon^2 + \Delta E^2} \right)\sin^2\left(\left(\epsilon^2 + \Delta E^2\right)^{1/2} t\right),\nonumber\\ & \approx \left(\frac{\epsilon^2}{\Delta E^2} \right)\sin^2\left( \Delta E \, t\right), \end{align} \tag{ 18 }$$

as an antineutron state, where $2 \Delta E$ is the difference in neutron and antineutron energies in external fields and $\epsilon = \tau_{n\bar{n}}^{-1}$ is the energy characteristic of whatever new physics is responsible for the oscillations; $\tau_{n\bar{n}}$ is referred to as the oscillation time. We use units in which $\hbar = 1$. The experimental conditions are such that $\Delta E \gg \epsilon$, which justifies the approximation in equation (18). Furthermore, in the Grenoble experiment the quasi-free condition $\Delta E t \ll 1$ could be realized, leading to a transition probability,

$$\begin{align} P_t & \approx \left(t \, / \, \tau_{n\bar{n}}\right)^2, \end{align} \tag{ 19 }$$

independent of the energy perturbation $\Delta E$ arising from the neutron and antineutron interactions with the ambient magnetic field. In the quasi-free condition N = 3 events were recorded in this experiment.

The background in the Grenoble experiment was directly measured by applying a magnetic field to suppress the transition probability P_t by making $\Delta E$ large enough. This condition yielded M = 7 events. The maximum likelihood estimate of the signal is $\hat{\mu} = N - M = -4$ events. However, since $\mu \unicode{x2A7E} 0$, we choose to take the best estimate of the signal in such experiments to be

$$\begin{equation} \hat{\mu} = \left\{ \begin{array}{ll} N-M & \text{if } \quad N\gt M \\ 0 & \quad \textrm{otherwise.} \end{array} \right. \end{equation} \tag{ 20 }$$

The sparsity of the data in the Grenoble experiment and our choice of signal estimate explicitly violate two of the regularity conditions for standard asymptotic results to hold [7]: the data should be sufficiently numerous and estimates must not lie on the boundary of the parameter space. The first condition is violated by the limited number of observations, and the second is violated by $\hat{\mu} = 0$, which lies on the boundary of the parameter space $\hat{\mu}\in [0,\infty )$. The violation of these regularity conditions, however, is not a problem for LF2I and ALFFI.

To construct confidence sets in the parameter space of $\theta = \{\mu, \nu \}$, we use the test statistic

$$\begin{equation} \lambda\left(D ; \theta\right) = -2 \log \left[ \frac{\mathcal{L}\left(D ; \mu, \nu\right)}{\mathcal{L}\left(D ; \hat{\mu}, \hat{\nu}\right)} \right], \end{equation} \tag{ 21 }$$

where $\hat{\mu}$ is given by equation (20) and $\hat{\nu}(\mu)$ by

$$\begin{equation} \hat{\nu} = \left\{ \begin{array}{ll} M & \text{if } \quad \hat{\mu} = N - M \\ \left(M+N\right)/2 & \quad \textrm{otherwise.} \end{array} \right. \end{equation} \tag{ 22 }$$

The cdf, $\mathbb{P}(\lambda \unicode{x2A7D} \lambda_D | \theta)$, was again approximated with a fully-connected feed-forward DNN with 3 input features $\mathbf{x} = \{\lambda_D, \mu, \nu \}$, 6 hidden layers with 12 nodes each, and a single output, estimating $\mathbb{E} \left[ Z\mid\lambda_D,\mu,\nu \right]$. The activation function at each hidden node is a PReLU [24], and the network was trained with the Adam optimizer with a fixed learning rate of $6 \times 10^{-4}$. The training set is composed of 10⁷ examples, which were used in batches of size $5\times 10^3$, for the duration of 10⁵ iterations, that is, for 50 epochs. A batch normalization [25] layer was added after every hidden layer, which was found to improve the results.

The DNN was used to compute the confidence sets shown in figure 4 and the associated coverage probabilities shown in figure 5. The fact that the coverage probabilities are close to the desired CLs confirms the accuracy of the confidence sets obtained with ALFFI.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the cosmology and signal/background examples, the likelihood functions are tractable. In our third example, from the field of epidemiology, the likelihood is intractable [26]. Therefore, this example is one for which LF2I and ALFFI are the most useful.

In this example, we fit the well-studied SIR epidemiological model to a classic data set, reproduced in figure 6, from a flu outbreak at an English Boarding School [27]. The SIR model has a nearly 100-year history, beginning with the work of Kermack and McKendrick in 1927 [28]. The simplest version of the model comprises coupled ordinary differential equations (ODEs) and assumes a population of individuals that is closed and well-mixed in which individuals fall into one of three classes or compartments: susceptible (S), infectious (I), and recovered or removed (R). The rate of change of susceptible individuals depends on the rate at which new infections arise as a result of contact between infectious and susceptible individuals. It is typically assumed that the contact rate (number of contacts per unit time) is proportional to the total population size N or constant. The assumption that contacts are proportional to the total population size leads to a mass action term $\beta S I$ for the number of new infections per unit time (called the incidence in epidemiology). Infectious (I) individuals can leave the infectious class through recovery (or death) and progress to the the Recovered (or Removed) class, R. It is often assumed that the number of recoveries/removals per unit time is linear in I: $\alpha I$. Under this assumption, the time spent in I is exponentially distributed with mean $1/\alpha$. The resulting system of ODEs is

Zoom In Zoom Out Reset image size

$$\begin{align} \frac{\textrm{d}S}{\textrm{d}t}& = -\beta SI,\nonumber\\ \frac{\textrm{d}I}{\textrm{d}t}& = -\alpha I + \beta SI,\nonumber\\ \frac{\textrm{d}R}{\textrm{d}t}& = \alpha I. \end{align} \tag{ 23 }$$

The qualitative dynamics of this model are governed by an epidemiological quantity called the basic reproduction number, R₀, which for this SIR model is $R_0 = \beta/\alpha$. If $R_0\gt 1$, the model exhibits a single outbreak, with I increasing to a peak, then approaching zero as time tends to infinity. Under this scenario, S will decay and approach a positive, constant value, meaning that the epidemic does not infect the whole population. If $R_0 \unicode{x2A7D} 1$, I will decay and approach zero as time tends to infinity. This model has played an invaluable role in understanding infectious disease dynamics, informing control strategies, and serving as a building block for more complex compartmental epidemiological modeling frameworks.

The SIR model is often fitted to data by minimizing a weighted least squares function,

$$\begin{align} F\left(\theta\right) = \sum_{n = 1}^N w_n \left(x_n - I_n\right)^2 , \end{align} \tag{ 24 }$$

where the data $D = \{x_1,\ldots,x_N \}$ are the number of infected individuals at the corresponding reporting times $t_1,\ldots,t_N$, $I_n = I(t_n, \theta)$ is the predicted mean infection count at time t_n , found by solving equations (23) for the given θ, and w_n are weights. Following example 1, we choose a test statistic that is transformed and scaled

$$\begin{align} \lambda\left(D, \theta\right) & = \sqrt{F\left(\theta\right) / N} \, / \, 50 , \end{align} \tag{ 25 }$$

so that the statistic is ${\cal O}(1)$. The weights are set to $w_n = I_n^{-1}$, which we hasten to add should not be taken to imply that the counts x_n are Poisson distributed. On the contrary, the counts x_n are correlated and their fluctuations are super-Poissonian.

We follow a similar training protocol as for example 1, except that the training sample size for this example is 750 000, a fully-connected DNN is used with 6 hidden layers and 25 nodes each, and the best model is found after about 350 000 iterations, that is, after 23 epochs. The confidence sets obtained are shown in figure 7 and the coverage probabilities are shown in figure 8. Again, we find that ALFFI produces accurate confidence sets.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size