Probabilistic Recalibration of Forecasts

TL;DR

This paper introduces a recalibration method for probabilistic forecasts using Gaussian processes, improving calibration and enabling consistent betting advantages, demonstrated through experiments with nonlinear circuits and climate forecasts.

Contribution

It proposes a novel Gaussian process-based recalibration scheme that enhances probabilistic forecast calibration and betting performance, even with dependent PIT values.

Findings

Recalibrated forecasts achieve uniform PIT distribution.

The scheme improves betting success against original forecasts.

Effective in both laboratory and climate forecast case studies.

Abstract

We present a scheme by which a probabilistic forecasting system whose predictions have poor probabilistic calibration may be recalibrated by incorporating past performance information to produce a new forecasting system that is demonstrably superior to the original, in that one may use it to consistently win wagers against someone using the original system. The scheme utilizes Gaussian process (GP) modeling to estimate a probability distribution over the Probability Integral Transform (PIT) of a scalar predictand. The GP density estimate gives closed-form access to information entropy measures associated with the estimated distribution, which allows prediction of winnings in wagers against the base forecasting system. A separate consequence of the procedure is that the recalibrated forecast has a uniform expected PIT distribution. A distinguishing feature of the procedure is that it is appropriate even if the PIT values are not i.i.d. The recalibration scheme is formulated in a framework that exploits the deep connections between information theory, forecasting, and betting. We demonstrate the effectiveness of the scheme in two case studies: a laboratory experiment with a nonlinear circuit and seasonal forecasts of the intensity of the El Ni\~no-Southern Oscillation phenomenon.

Tables

Table 1. Table 1: NMME models selected for this study, and respective ensemble sizes.

Model	Ensemble Size
COLA-RSMAS-CCSM3	6
COLA-RSMAS-CCSM4	10
GFDL-CM2p1-aer04	10
GFDL-CM2p5-FLOR-A06	12
GFDL-CM2p5-FLOR-B01	12

Equations

$$\tilde{F}(x;J,C)={\int }_{-\infty }^{x}d{x}^{\prime }p(x^{\prime };J,C).$$

$$\pi \left(X_n=x_n|{\mathcal{J}}_n=J_n,C\right)$$

$$\pi \left(F|C\right)$$

$$\pi \left(F|\mathcal{F},C\right)$$

$$=$$

$$p_1(x;J,C)=\pi \left(F=\tilde{F}(x;J,C)|\mathcal{F},C\right)\times p(x;J,C),$$

$$KL[f_2||f_1]=\int dx{f}_2(x){\log }_2\frac{f_2(x)}{f_1(x)},$$

$$\Delta s$$

$$E_{\Pi (F|C)|\mathcal{F}}\left\{E_{\mathcal{J}}\left[\Delta s\right]\right\}$$

$$\tilde{G}(f;C)$$

$$\Pi (G=g|\mathcal{J},C)$$

$$\Pi (G=g|C)$$

$$E_{\Pi (F|C)|\mathcal{F}}\left\{\Pi (G=g|C)\right\}$$

$$\Delta S_{True}\equiv E_{\mathcal{J}}\left\{\Delta s_{True}\right\}={\int }_0^{1}df\pi (F=f|C){\log }_2\pi (F=f|\mathcal{F},C),$$

$$\text{Ign}[p]$$

$$\text{Ign[p]}\approx -\frac{1}{N}\sum _{n=1}^N{\log }_2\frac{p(x_n;J_n,C)}{\pi (X=x_n|C)},$$

$$\Delta \text{Ign}[p_1,p]$$

$$\Delta I$$

$$-\Delta \text{Ign}[p_1,p]$$

$$\Delta S_{True}$$

$$\Delta S\equiv E_{\mathcal{J}}\left\{\Delta s\right\}={\int }_0^{1}df\Pi (F=f|C){\log }_2\pi \left(F=f|\mathcal{F},C\right),$$

$$\overline{\Delta S}\equiv E_{\Pi (F|C)|\mathcal{F}}\left\{\Delta S\right\}$$

$$\text{Var}\left(\Delta S\right)\equiv E_{\Pi (F|C)|\mathcal{F}}\left\{\Delta S^2\right\}-{\left(\overline{\Delta S}\right)}^2.$$

$$\overline{\Delta S}$$

$$\text{Var}(\Delta S)$$

$$\text{FAM}=\frac{\overline{\Delta S}}{\sqrt{\text{Var}\left(\Delta S\right)}},$$

$$EI\left[\pi (F|\mathcal{F},C)\right]\equiv E_{\Pi (F|C)|\mathcal{F}}\left\{{\int }_0^{1}df\Pi (F=f|C){\log }_2\frac{\Pi (F=f|C)}{\pi (F=f|\mathcal{F},C)}\right\}.$$

$$EI\left[\pi (F|\mathcal{F},C)\right]=\frac{1}{2\ln 2}{\int }_0^{1}df\pi (F=f|\mathcal{F},C)\times C(f,f),$$

$$\dot{x}$$

$$\dot{y}$$

$$\dot{z}$$

$$\ln \pi =-n+n\ln n-\frac{1}{2}\frac{1}{{\sigma }_{\mu }^2}(\mu -\bar{\mu }{)}^2+R_1(n,\mu ),$$

Full text

1 Introduction

A forecast, being an expression of uncertainty about the future, is necessarily a probabilistic affair. Probabilistic forecasts of events falling along a continuum—such as short-term weather forecasts [1, 2], medium-term seasonal rainfall [3, 4, 5], fluctuations of financial asset prices [6, 7] or electrical demand [8], rates of spread of infectious disease [9, 10], macroeconomic indicators [11, 12], wind power availability [13, 14], species endangerment and extinction [15], human population growth [16], and seismic activity [17, 18]—are of urgent interest to many kinds of decision-makers, and have occasioned much scientific literature across a broad range of fields.

Weather forecasting through numerical weather prediction (NWP) has substantially improved its performance over the past few decades in consequence of improvements in observational data, computational models, and computational power and currently is capable of providing, on average, reasonably robust weather forecasts for periods on the order of 10 days [19]. Unfortunately, compared to empirical forecasting schemes, these physics-based forecasting schemes are known to lose skill for periods longer than $\approx$10 days, and the obvious question arises whether it is possible to achieve skillful forecasting for periods longer than 10 days, and if so, whether there is an upper bound on such forecasting.

Given the difficulty of achieving the computing power and model fidelity required to abate model errors of NWP, the best way forward for the present may well be to attempt to achieve some kind of synthetic hybrid between NWP and empirical forecast methods in an effort to leverage the information content of the former to enhance the predictive power of the latter. Important approaches that have been attempted include applying some kind of statistical recalibration to the NWP simulation output such as Model Output Statistics [20] and then infer probability distributions for predictands from the corrected simulations by some smoothing procedure [21, 22]; leaving the simulations as they are and adapt the smoothing procedure itself to validation data [23]; or doing a bit of both [24, 25, 26]. A difficulty of such programs is that the smoothing procedure itself has statistical properties that are usually not under very good control, since they frequently take the form of highly simplified models such as Gaussian mixtures, which are not generally in well-motivated correspondence with the processes that relate the simulation output to the random predictand.

One common feature of the continuous forecast probability density functions (PDFs) produced from NWP output is that they are more often than not probabilistically miscalibrated—that is, the long-term frequencies of observations do not match stated probabilities of predictions (see [27], for example). In such cases, the interpretation of the PDFs requires caution, and the value of having a probabilistic forecast rather than a point forecast can be questionable, particularly given the risk of underestimating frequencies of extremes.

As discussed by Diebold, Hahn, and Tay ([28], hereafter “DHT”), the phenomenon of probabilistic miscalibration creates another opportunity for recalibration: direct recalibration of the forecast probability distributions. This possibility arises because whatever the methodology adopted to produce a forecast system, long enough use of that system leads to additional performance information—through comparison of a series of forecasts with their predictands—that can be incorporated into current forecasts to produce improved forecasts. Such information, which is commonly used to assess forecast system quality, was shown by DHT [28] to permit correction of future forecasts, assuming an i.i.d. restriction on predictands. More recently, a similar approach was devised in the context of deep learning by Kuleshov, Fenner, and Ermon ([29], hereafter “KFE”), who recast the usual machine learning activities of regression and classification as forecasting problems, and regarded artificial neural network outputs as discrete or continuous predictands, respectively. Under an i.i.d. restriction on outputs, KFE [29] obtain recalibration procedures that are equivalent to those of DHT [28], now cast as calibration procedures for classification and regression.

In this work we generalize the work described in [28, 29] in two respects. In the first place, we establish a mathematical framework for treating predictands without i.i.d. restrictions. In the process of doing so, we strongly emphasize the role of conditional information in forecast distributions, and make use of ideas from information theory to characterize the effect of miscalibration. Additionally, we replace the PDF estimation schemes suggested in [28] and the isotonic regression adopted by [29] with a Gaussian-process (GP) density estimation scheme, which allows us to also estimate – with quantified uncertainties – information entropy measures associated with the estimated pdfs. Using this technique we demonstrate that if a series of forecasts shows evidence of poor probabilistic calibration then we may use past forecast performance information to produce new current forecasts that have well-calibrated expected distributions, and that have greater expected logarithmic forecast skill score than the original forecasts, irrespective of whether the predictands are i.i.d. Furthermore, the expected performance improvement in logarithmic skill score is computable in advance, together with an uncertainty estimate.

We demonstrate the method using two case studies: a laboratory experiment with a nonlinear circuit and seasonal forecasts of the intensity of the El Niño-Southern Oscillation (ENSO) phenomenon.

2 Probabilistic Forecasts

A probabilistic forecast of a continuous scalar random variable $X$ is simply a probability distribution $P(X|I)$ over the value of the predictand $X$, which is to be observed at a later date. The distribution is conditioned on prior information $I$ such as current and past conditions, (often approximate) deterministic and probabilistic model structure, empirically determined training parameters, and simulation output. Such forecasts are often generated as a time series $P(X_{n}|I_{n},C)$, with $n\in\mathbb{Z}$ an index that labels time $t_{n}$, so that $n_{2}>n_{1}\implies t_{n_{2}}>t_{n_{1}}$. Here, $X_{n}$ is the random predictand at time $t_{n}$, $I_{n}$ represents information that varies with $n$, while $C$ represents static conditioning information that is constant for a particular forecasting system. Typically, the information $I_{n}$ is stochastic, and fluctuates randomly with $n$. Consequently, the distribution $P(X_{n}|I_{n},C)$ is itself a distribution-valued random variable [30, 1]. Note that implicit in the notation $P(X_{n}|I_{n},C)$ is the assumption that the particular realization $I_{n}=i$ completely determines the distribution of $X_{n}$ irrespective of $n$, so that $P(X_{n}|I_{n}=i,C)=P(X_{m}|I_{m}=i,C)$ for $m\neq n$. This will allow us to consistently drop time subscripts from expressions such as $P(X|I,C)$ in what follows.

As an example, in the case of weather prediction, $I_{n}$ might represent a discrete vector of weather observations at a finite number of weather stations over the course of several previous days, while $C$ might represent climatological data. Another example is provided by the empirical time-series modeling that underlies many analyses of financial and economic data, where the $I_{n}$ could be the last $M$ values of the time series $X_{n}$ and $C$ the parameters of an autoregressive moving-average (ARMA) time-series model [31].

In our methodological development we assume the system is approximately stationary, so that secular drifts due to external forcings are ignored. We also overlook annual-type periodicities, which are in principle tractable by adding cycle phase information to $I$. One may easily show that a unique distribution $P(X|I,C)$ always exists in principle. This follows simply from the existence of a unique joint distribution $P(X,I|C)$, which is ascertainable empirically from a sufficiently large archive of $(X_{n},I_{n})$ values. The forecast distribution $P(X|I,C)$ is then just $P(X,I|C)/P(I|C)$. This unique distribution is called the ideal forecast with respect to the information $I$ [30, 1].

Above and beyond empirical observation, often some kind of dynamical law exists from which $P(X|I,C)$ could in principle be inferred. In such cases, however, accurate inference of $P(X|I,C)$ from first principles is often impractical, because either the dynamical law is not known (as in the case of most time series in economics) or it is known imperfectly (as is the case with weather forecasting), or it is not feasibly computable even where it is well understood.

Weather forecasting furnishes an instructive example. The dynamical origin of the distribution $P(X|I,C)$ is intelligible in terms of the nonlinear physics of weather systems. However, while describable, this forecast distribution is in no way feasibly computable, because of limitations in model fidelity and in computational resources. Instead, limited-fidelity computational models [32, 33] are used to filter the information in $I$, incorporating techniques of data assimilation [34], evolving ensembles of states not chosen by a fair sampling of the distribution on the observation-constrained submanifold of the chaotic attractor, and in any case with too few ensemble members to be sufficiently informative about the distribution’s structure. Postprocessing of a training set of ensembles and corresponding validation values of X must be used to construct an approximation to $P\left(X|I,C\right)$ [24, 35].

Clearly, by the time this approximation has been constructed, it is no longer necessarily conditioned directly on $I$, but rather on some highly processed information $J[I]$. In ensemble NWP, $J[I]$ has both a deterministic aspect (the NWP simulations) and a stochastic aspect (the selection of the random ensemble of initial conditions to evolve). Quite generally, we can assume that some probabilistic model $J\sim P(J|I)$ describes the dependence of $J$ on $I$. If that mapping should happen to be deterministic, the probability distribution $P(J|I)$ would degenerate to a product of Dirac $\delta$-distributions. In general, the dimensionality of $J$ is not necessarily inferior to that of $I$—in the case of NWP, the simulations generate data over grids whose data mass far exceeds that of the input information. Invariably, however, the information content of $J$ is degraded in comparison with that of $I$, by the very approximations described above. This is merely the observation that $P(X|J[I],C)$ is expected to be—and generally is—inferior to the computationally infeasible $P(X|I,C)$ in quality measures such as calibration and sharpness (discussed below).

Another practical concern in obtaining $P(X|I,C)$ is the fact that the static conditioning information $C$ may not be exactly known and must be estimated from data. For example, even if a financial time series were known to be well approximated by a stationary autoregressive process, the process parameters would not in general be known, and would have to be fit from data. Similarly, in weather, climatological information would have to be fit from noisy data.

We will assume that forecasts are appropriately modeled by absolutely continuous distributions, which may therefore be represented by probability densities over $X$. The forecasting system converts the information $J_{n}=J[I_{n}]$ and $C$ into a published forecast $p(X_{n};J_{n},C)$, a density over $X_{n}$, at each time $n$. The data-generating process that is being forecast then generates a realization $X_{n}=x_{n}$. Suppose we generate $N$ forecasts and $N$ corresponding observations. The series of pairs $\left\{{\cal P}_{n}=\left(x_{n},p(\cdot;J_{n},C)\right),n=1,2,\ldots N\right\}$ is called the *Forecast-Observation Archive *(FOA) [36, 37]. The pairs ${\cal P}_{n}$ may be viewed as elements of a set called a prediction space, whose mathematical properties were analyzed in [1].

One of the most important tools for assessing the validity of published forecasts is the Probability Integral Transform, or PIT [38, 6, 39]. This is defined in terms of the cumulative probability distribution function $\tilde{F}(\cdot;J,C)$ associated with $p(\cdot;J,C)$,

[TABLE]

The PIT associated with the pair ${\cal P}_{n}=(X_{n}=x_{n},p(\cdot;J_{n},C))$ is simply the value $f_{n}\equiv\tilde{F}(x_{n};J_{n},C)$. The reason for the usefulness of the PIT is that if the density $p(\cdot;J_{n},C)$ correctly models the stochastic behavior of $X_{n}|J_{n},C$, then the random variable $F_{n}=\tilde{F}(X_{n};J_{n},C)$ must be uniformly distributed over the interval $\left[0,1\right]$ irrespective of $J_{n}$. If this is the case, we say that $p(\cdot;J,C)$ is *probabilistically calibrated *[30, 1]. Probabilistic calibration is a desirable feature in a published forecast because it means that the forecast is “honest” about the probabilities of its quantiles, since those probabilities correspond to long-term average frequencies. The property of being probabilistically calibrated may be checked, given a sufficiently large FOA, by histogramming the values of $F_{n}$ and inspecting the histogram for evidence of nonuniformity [38, 6, 39]. All ideal forecasts are probabilistically calibrated, although the reverse is not true—many different calibrated forecasts can easily be constructed, but only one is ideal with respect to the input information.

It is perhaps surprising to realize that calibration, while a desirable feature of a forecast system, is not sufficient to prefer one forecast system to another. As discussed in [39, 40], it is quite possibly for forecast systems yielding distributions that vary widely in precision to all be equally probabilistically calibrated. For example, a “climatological” forecast system that uses only long-term historical averages to make predictions and an idealized perfect NWP model that makes approximation-free use of current information $I_{n}$ to make ideal forecasts are both equally probabilistically calibrated from the point of view of PIT histogram uniformity. Clearly, the former provides forecasts that are vague compared with those of the latter, which are more informative and precise. The term sharpness was introduced by Bross and Bross [41] to characterize this distinction. It refers to the degree of concentration on small outcome sets of the published forecast density, and is sometimes expressed as distributional variance or as width of a central fixed-probability (e.g., 90%) interval [39]. Thus a sharper forecast is less vague in its predictions than is a less-sharp one, independently of the relative degree to which their respective PIT histograms are close to uniform.

Clearly, the difference between probabilistically calibrated forecasts of different sharpness —the difference between the climatological and the idealized NWP forecaster, for example—is purely in the information on which the forecasts are conditioned. States of more specific information lead to sharper forecasts. In the case of the idealized NWP forecaster, for example, a large increase in the number of available weather stations necessarily leads to sharper forecasts, while an increase in the measurement uncertainty of current weather conditions necessarily leads to less-sharp forecasts. The explicit highlighting of the relevant conditioning information is therefore essential to the discussion of probabilistic forecasting.

In the case of poorly calibrated forecasts, the source of the misspecification of the forecast distributions is necessarily to be sought in erroneous conditioning information, such as model errors that distort the information borne by the processed input data $J[I]$, model errors in constructing the published forecast distribution, or poor approximations encoded in the static conditioning information $C$. Forecast interpretation in the presence of misinformation is an important subject in decision support [42]. One may be confronted with cases of sharp but uncalibrated forecasts, that are (for example) biased, but which have smaller mean-square error than climatology. In these cases the incorrect conditioning information may not entirely be condemned, because such forecasts can have better predictive skill than climatology. One naturally wonders about the extent to which this partially correct information can be exploited to produce probabilistically calibrated forecast distributions.

In the next section, we will show that knowledge of the PIT histogram of a sufficiently large FOA can be used to correct a current published forecast $p(\cdot;J_{n},C)$ prior to the observation of the predictand $X_{n}$, to produce a new, updated forecast $p_{1}(\cdot;J_{n},C)$ that outperforms $p(\cdot;J_{n},C)$, in that it has a better expected logarithmic (“ignorance”) score (the ignorance score is defined in [35], and is further discussed below in §3.3). This probabilistic recalibration procedure allows us to better exploit the correct part of the conditioning information.

3 Probabilistic Recalibration

The essence of the probabilistic recalibration procedure is that the PIT histogram of a sufficiently large FOA can be subjected to an empirical fit, so that the underlying distribution may be inferred by regression. Assuming that current forecasts suffer from the same miscalibration as those in the FOA, the fit distribution may then be used to correct a current forecast distribution to produce a new forecast that outperforms the original by various objective measures, including the ignorance score. We now set out the procedure.

Since we have raised the issue of misinformation in connection with miscalibration, we adopt notation that distinguishes between distributions that are ideal—that correctly reflect their conditioning information, that is—and distributions that may be misinformed or poorly calibrated. In what follows, therefore, we will reserve the symbol $\pi$ and the notation $\pi\left(A|B\right)$ or $\pi(A=a|B)da$ for the probability density function of a random variable $A$ correctly conditioned on information $B$, so that $\pi(A|B)$ is ideal. Published forecast densities, which may be “misinformed” and hence incorrectly reflect the dependence on conditioning information, we denote by simple function notation such as $p(x;J)$.

We denote the input information to the $t=t_{n}$ published forecast by the random variable ${\cal J}_{n}$, whose realization is $J_{n}$. We assume the existence of a unique ideal distribution relative to $J,C$ with density $\pi\left(X|{\cal J}=J,C\right)$. Again, such a unique forecast clearly exists, by its relation to the empirically ascertainable joint distribution $\pi(X,{\cal J}|C)$.

A published forecast $p(\cdot;J_{n},C)$ and the corresponding observations $x_{n}$ give rise to a PIT value $f_{n}=\tilde{F}(x_{n},J_{n},C)$, which is a realization of a random variable $F_{n}\equiv\tilde{F}(X_{n},J_{n},C)$. The variable $F_{n}$ is simply a change of random variables from $X_{n}$, which implies an ideal density $\pi(F_{n}|{\cal J}=J_{n},C)$ for $F_{n}$ satisfying

[TABLE]

Equation (2) connects the published forecast $p(\cdot;J_{n},C)$ to $\pi\left(X_{n}|{\cal J}_{n}=J_{n},C\right)$, the unknown ideal forecast distribution density relative to $J_{n}$, by a pointwise multiplication with another unknown density $\pi\left(F_{n}|{\cal J}_{n}=J_{n},C\right)$. 111For the sake of simplicity, we assume that published forecasts $p(x;J,C)$ are always non-zero for all $x$. If this were not the case, an interval $[x_{1},x_{2}]$ over which $p(x;J,C)$ is zero would be mapped to a single point $F_{1}$ by the $x\rightarrow F$ change of variables. Should the ideal distribution density $\pi\left(X|\mathcal{J},C\right)$ happen to have nonzero probability mass over such an interval, that finite probability would be mapped to the single point $F_{1}$. It would then be necessary to represent this effect by an additive Dirac-$\delta$ distributional component in $\pi(F|\mathcal{J},C)$. Such a generalization would be cumbersome, and we avoid having to address it by specifying $p(x;J_{n},C)>0$ for all $x$.

If the observables $F_{n}$ are i.i.d., we may drop the dependence of $\pi\left(F_{n}|\mathcal{J}_{n},C\right)$ on $\mathcal{J}_{n}$ in Equation (2). In this case one may proceed straightforwardly estimating the time-independent distribution $\pi(F|C)$ (where $F$ may be any of the identically-distributed $F_{n}$) by regression on the FOA PIT data $\mathcal{F}\equiv\left\{F_{n}:\,n=1,2,\ldots,N\right\}$ as described in DHT [28]. Denoting this estimate by $\pi\left(F|{\cal F},C\right)$ and replacing $\pi\left(F_{n}|\mathcal{J}_{n},C\right)$ by $\pi\left(F_{n}|{\cal F},C\right)$ in Equation (2) results in forecasts with improved calibration properties. Equivalently, KFE [29] perform isotonic regression on what is, in effect, the i.i.d CDF of the $F_{n}$ (as opposed to their i.i.d. PDF), to obtain improved probabilistic calibration of deep learning classifiers and regressors. The theory developed in [28, 29] does not address the important general case of non-i.i.d. $F_{n}$, however, because the regression estimate $\pi\left(F|\mathcal{F},C\right)$ from the FOA averages over all conditioning data $\mathcal{J}$, and in this sense is a “climatological” distribution that is ignorant of current conditionining information $\mathcal{J}_{n}=J_{n}$.222Note that the i.i.d assumption in [28, 29] applies to the distribution of the $F_{n}$, and not to the forecast distribution of the $X_{n}$. The latter are not “climatological” in that they are conditioned on their individual $\mathcal{J}_{n}$. DFE [28] address the i.i.d. restriction by considering published and ideal forecasts belonging to different “scale-location” families with the same “scale-location” parameters, showing that this case does give rise to i.i.d. $F_{n}$. The generality of this restriction is problematic, however. As we discuss below in §3.6, and explicitly demonstrate in §4.2, it is not uncommon for time-series of $F_{n}=f_{n}$ realizations to be exhibit strong correlations. In such cases the i.i.d. assumption on the $F_{n}$ is simply not tenable.

Equation (2) is nonetheless the starting point of our probabilistic recalibration procedure: we will show that if in Equation (2) we replace the density $\pi\left(F|{\cal J},C\right)$ with the predictive distribution density $\pi\left(F|{\cal F},C\right)$ estimated by a Bayesian regression fit to the FOA PIT data ${\cal F}$, we will obtain a new forecast distribution which is not ideal, but which nonetheless improves on the logarithmic (“ignorance”) skill score of the published forecast irrespective of whether the $F$ data are i.i.d.

3.1 Bayesian PIT-Fit

We now perform the regression fit to the data ${\cal F}$. Diebold et al. [28] recommended either using a kernel estimator, or simply and directly the empirical PIT distribution, whereas Kuleshov et al. [29] recommended isotonic regression on the PIT CDF. Here we adopt a non-parametric procedure: Gaussian process measure estimation (GPME), described in detail in A. This is a more complicated procedure than previously used for this task, but it has benefits that will be described presently. As used here, GPME effectively fits functions from an infinite-dimensional function space to estimate the predictive density $\pi\left(F|{\cal F},C\right)$.

Our stationarity assumption implies that the FOA PIT data ${\cal F}$ may be viewed as a realization of a process that repeatedly samples a climatologically averaged distribution, corresponding to many different random realizations $J$ of the random variable ${\cal J}$ representing the input information. The climatology gives rise to a distribution $\pi\left({\cal J}|C\right)$, which we use to average $\pi\left(F|J,C\right)$, obtaining

[TABLE]

The distribution $\pi(F|C)$ is unknown and must be estimated by regression from the noisy FOA PIT data in ${\cal F}$. This density function estimate bears uncertainty represented by the Gaussian process posterior distribution over the density function $\pi(F|C)$ given ${\cal F}$. This uncertainty is purely epistemic, in contrast to the uncertainty consequent on the stochastic nature of ${\cal J}$, represented by the climatological distribution $\pi({\cal J}|C)$.

We will notationally represent the uncertainty in the determination of the ideal density function $\pi(X|{\cal J},C)$ using a density function-valued random variable $\Pi(X|{\cal J},C)$, whose realizations are possible densities $\pi(X|{\cal J},C)$. We refer to density function-valued random variables such as $\Pi(X|{\cal J},C)$ as imperfectly known distributions. The distribution $\Pi(X|{\cal J},C)$ is imperfectly known because our knowledge of it comes from a database of time series of pairs (${\cal J}$,$X$), from which the joint distribution $\pi(X,{\cal J}|C)$, and hence $\pi(X|{\cal J},C)$, could be estimated empirically, with considerable uncertainty.

One may use Equation (2) to define the imperfectly known distribution $\Pi(F=f|{\cal J}=J,C)=\Pi\left(X=x|{\cal J}=J,C\right)/p\left(x;J,C\right)$, where $\tilde{F}(x;J,C)=f$, whose realizations are possible densities $\pi(F|{\cal J},C)$.

The uncertainty in $\pi(F|C)$ is then represented by an imperfectly known distribution $\Pi(F|C)\equiv E_{{\cal J}}\left[\Pi(F|{\cal J},C)\right]$, whose realizations are possible density functions $\pi(F|C)$. The prior distribution over $\Pi(F|C)$ is described in GPME by a Gaussian process over $\log\Pi(F|C)$ with a chosen kernel $K(f_{1},f_{2})$ (here squared-exponential) and a constant mean function. The posterior distribution over $\Pi(F|C)$ given ${\cal F}$ is described in GPME by an updated Gaussian process over $\log\Pi(F|C)$, with a mean function $\lambda(f)$ given by Equation (56), and a covariance $C(f_{1},f_{2})$ given by Equation (57).

The predictive distribution $\pi(F|{\cal F},C)$ is the expectation of $\Pi(F|C)$ under this posterior distribution $\Pi(F|C)|{\cal F}$, that is

[TABLE]

Equation (4) expresses the operation by which $\pi\left(F|{\cal F},C\right)$ is obtained from the GPME posterior distribution. Equation (5) illustrates the fact that the predictive distribution incorporates both epistemic fit uncertainties and climatological averaging. This distribution, which is the key quantity enabling the probabilistic recalibration procedure, is given explicitly in Equation (65) and is the principal output of the GPME procedure.

3.2 The Recalibration Procedure

As adumbrated above, the recalibration procedure consists in replacing the published forecast density $p(x;J,C)$ with the recalibrated forecast density $p_{1}(x;J,C)$ given by the recalibration equation

[TABLE]

which is obtained from Equation (2) simply through the replacement of $\pi(F|{\cal J},C)$ by the predictive distribution ${\cal\pi}(F|{\cal F},C)$, estimated by the GPME procedure.

We may gauge how much has been gained by the recalibration procedure, by introducing the Kullback-Leibler divergence or relative entropy [43] of a density $f_{1}(x)$ relative to another density $f_{2}(x)$,

[TABLE]

which may be usefully viewed as a measure of the information embodied by $f_{2}$ relative to a prior state of information that is embodied by $f_{1}$. The divergence $KL[f_{2}\,||\,f_{1}]$ has the well-known property of being non-negative definite, and of being zero only if $f_{1}=f_{2}$ almost everywhere.

By setting $f_{2}$ to the imperfectly known ideal distribution $\Pi(X|{\cal J}=J,C)$ and setting $f_{1}$ alternatively to the published forecast $p(x;J,C)$ and to the recalibrated forecast $p_{1}(x;J,C)$, we may, by taking expectations over the climatology ${\cal J}$ and over the posterior GPME model of the density $\Pi(F|C)\,|\,{\cal F}$, assess whether the ideal distribution is closer in information to the recalibrated forecast than the original published forecast.

This leads to the following theorem:

Theorem 1**.**

(i) The distribution $\Pi\left(X=x|{\cal J}=J,C\right)$ is closer to the recalibrated forecast $p_{1}(x;J,C)$ in expected (over ${\cal J}$ and $\Pi(F|C)\,|\,{\cal F}$) relative entropy than it is to the published forecast $p(x;J,C)$ unless $p_{1}=p$ almost everywhere. (ii) The recalibrated forecast $p_{1}(x;J,C)$ is on average (over $\Pi(F|C)\,|\,{\cal F}$) probabilistically calibrated.

To prove (i), we first define the difference in the relative entropies,

[TABLE]

This quantity is a random variable in consequence of the stochastic nature of ${\cal J}$ and the epistemic uncertainty in $\Pi(F|{\cal J},C)\,|\,{\cal F}$. Taking the required expectations, we obtain

[TABLE]

with equality holding only when $\pi(F=f|{\cal F})=1$ almost everywhere, which is to say, when the original distribution was probabilistically calibrated. In this case, $p_{1}=p$.

To see (ii), observe that the cumulative distribution of $p_{1}(x;J,C)$ defines a change of variables from the random variable $F$ to a new random variable $G$ through the function $\tilde{G}(F;C)$ given by

[TABLE]

The function $\tilde{G}(\tilde{F}(x;J,C);C)$ is the PIT function of the recalibrated forecast distribution $p_{1}(x;J,C)$. In terms of the imperfectly known distribution $\Pi(F=f|{\cal J},C)$, the imperfectly known distribution $\Pi(G=g|{\cal J},C)$ is

[TABLE]

where $g=\tilde{G}(f;C)$. The imperfectly known PIT distribution for $p_{1}$ is obtained by averaging over ${\cal J}$:

[TABLE]

Averaging over the imperfectly known distribution $\Pi(F|C)\,|\,{\cal F}$ and using Equation (4), we find

[TABLE]

Hence, on average over $\Pi(F|C)\,|\,{\cal F}$, the PIT of the recalibrated forecast $p_{1}$ is uniform. $\square$.

Note a curious feature of the theorem: it does not use any fact about the GPME procedure, other than that it allows averaging over the posterior distribution $\Pi(F|C)\,|\,{\cal F}$. The statements of the theorem would be true given any such distribution, even one that is wildly incorrect about the likely shape of the true $\pi(F|C)$. If the posterior distribution over $\Pi(F|C)\,|\,{\cal F}$ did happen to be wildly wrong, however, the theorem, while still true, would no longer furnish the basis for a working recalibration procedure. The reason is that the actual data-generating process produces an observable value of $E_{{\cal J}}\left\{\Delta s_{True}\right\}$ given by

[TABLE]

where the subscript indicates that the true (unknown) distribution $\pi(F|C)$ is used in the average, and where the distinction between $\Delta s$ and $\Delta S$ is that the latter is averaged over $\mathcal{J}$. This expression differs from the expression in Equation (9) and is under no obligation to be non-negative definite. This expression can be expected to be strongly positive only when $\pi(F|{\cal F},C)$ approaches $\pi(F|C)$, which is to say when the fitting procedure results in a posterior distribution that is well concentrated near density functions that look a lot like $\pi(F|C)$. The success of the density-fitting regression procedure is therefore essential to the success of the recalibration procedure.

This observation applies equally to other styles of regression estimates for $\pi(F|{\cal F},C)$, including the kernel density and histogram estimators of DHT [28], and the isotonic regression of KFE [29]. So long as those procedures succeed in furnishing an estimate of $\pi(F|{\cal F},C)$ that is reasonably close to $\pi(F|C)$ and not too uncertain they too should produce recalibrated forecasts that are on average (over the uncertainty in those estimates) informationally closer to the ideal forecast than is the published forecast. The novel thing here is that we can now see that this ought to happen irrespective of whether the $F_{n}$ are i.i.d.. This is an important observation, since it was not previously clear to what extent correlations among the $F_{n}$ could be expected to damage or even vitiate recalibration. As a result, the work of KFE and DHT can be seen to have broader applicability than might otherwise have been believed to be the case.

The recalibrated forecast $p_{1}(x;J,C)$, while on-average probabilistically calibrated, is not in general the same distribution as the unique ideal distribution $\pi(X=x|{\cal J}=J,C)$, even when the size of the FOA is large and the uncertainty in $\Pi(F|C)$ is very small. With a sufficiently large database of pairs $(x,J)$, one could empirically estimate $\pi(X=x|{\cal J}=J,C)$ and discern its differences from $p_{1}(x;J,C)$. What the theorem establishes is that on average $p_{1}$ is a better approximation to the ideal forecast than $p$, in the sense that it is closer in information to the ideal forecast.

3.3 Scoring, Betting, and Information

A natural connection exists between the relative entropy difference $\Delta S$ and the ignorance score [35]. For continuous predictands $X$, the ignorance score $\textrm{Ign}[p]$ of a forecast distribution $p(x;J,C)$ is defined, using the climatological density $\pi(X|C)$ as a reference distribution, by the expression

[TABLE]

Because this final average is over the joint distribution on $X,{\cal J}$, the ignorance score may be estimated empirically from an FOA simply by the average

[TABLE]

whose expected value is the expression in Equation (15). Note that for i.i.d. predictands $x$, the expression on the RHS of Equation (16) is proportional to the log-likelihood for the model represented by $p(x;J,C)$, up to an additive constant. The resemblance is purely formal for non-i.i.d. predictands, however.

The difference in the ignorance scores of the published and recalibrated forecasts is

[TABLE]

where in the third line we have approximated $\pi(F|C)$ by $\pi(F|{\cal F},C)$ and in the last line we appealed to Theorem 1. We therefore expect that the recalibrated forecast $p_{1}$ will have a lower (i.e., better) ignorance score than the original published forecast. By the standards of the ignorance score, then, $p_{1}$ is an improvement on $p$.

It was pointed out in [35, 44] that the Ignorance Score has an interpretation as a tool for practical decision-making under uncertainty. The interpretation is couched in terms of a horse race, in which a bettor and an oddsmaker make optimal decisions about their choices with respect to the discrete possible outcomes of a horse-race. Kelly [45] described a strategy for a player allocating wealth to bets on $N$ outcomes $i=1,\ldots,N$ offering wealth multiplier odds $o_{i}$, assuming the player works with a forecast distribution $h_{i}$, $\sum_{i=1}^{N}h_{i}=1$. Kelly showed that the optimal strategy for the player is to allocate wealth $W_{i}$ to bet on outcome $i$ according to the rule $W_{i}=Wh_{i}$. By contrast, the optimal strategy for the pari-mutuel bookmaker with a forecast distribution $g_{i}$ is to set odds $o_{i}=1/g_{i}$. With these strategies for wagering and odds-setting, supposing the ideal forecast is $p_{i},$ the player’s wealth grows at the expected rate $2^{\Delta I}$, where

[TABLE]

is the difference between the entropy divergences of the two forecast distributions relative to the true distribution $p_{i}$, and hence also the negative Ignorance Score difference. This game can be regarded as a symmetric game between two forecasters who alternate the roles of bookmaker and bettor—so long as the current bookmaker with forecast $g_{i}$ sets odds in reciprocal proportion to $g_{i}$ and the bettor with forecast $h_{i}$ sets bets in proportion to $h_{i}$, Equation $($18) shows that it makes no difference which is which. A player with a forecast system that is consistently closer in information to the ideal distribution $p_{i}$ (and hence has a lower ignorance score) than the other player’s system has a positive expected wealth growth rate in this game.

This connection is attractive because it furnishes an example of decision-making under uncertainty that is improved by using forecasts with lower ignorance scores, and in particular by using the recalibration procedure described above. The view that underpins this work is that there ought to be some form of symmetric-rule game in which decisions made according to a higher-scoring forecasting system should consistently dominate those made under a lower-scoring one. This view of forecast quality differs from the standard “proper scoring” outlook [1, 46, 35, 47, 2], in which forecasts vie amongst one another for high scores using scoring rules designed to encourage forecaster honesty and furnish usable utility functions for estimation. It is in better correspondence with the Bayesian decision theory outlook on proper scoring [48, 46].

The Kelly horse race is not ideal for furnishing a fiducial symmetric-rule game for continuous probabilistic forecast assessment since it concerns itself with discrete outcomes and contains an appearance of asymmetry in the bettor-odds–maker model that it is based on. It can be adapted to continuous forecast distributions, for example by using fine quantiles (such as percentiles) as outcomes to be wagered on. A more natural symmetric-rule game can be described, however, by observing that we may recast the second line of Equation (17) as

[TABLE]

Observe that the right-hand side of Equation (19) may be interpreted as the expected winnings in a symmetric-rule game, which we call the entropy game. The rules of the game are as follows: two players compete, one using the published forecast distribution $p(\cdot;J,C)$, the other using the recalibrated distribution $p_{1}(\cdot;J,C)$. At the $n$-th turn of the game, values $J_{n}$,$x_{n}$ are observed from the data-generating distribution $\pi({\cal J},X|C)$. The players compute the log-ratio of their respective densities at $x_{n}$, $w=\log_{2}\left[p_{1}\left(x_{n};J_{n},C\right)/p(x_{n};J_{n},C)\right]$. If $w>0$, then the player using the published forecast $p$ pays the amount $w$ to the player using the recalibrated distribution $p_{1}$. Otherwise, the recalibrated player pays $|w|$ to the player using $p$. Equation (19) states that $\Delta S_{True}$ is the expected winnings per turn of the recalibrated player.

The entropy game may be viewed as the logarithm of the Kelly horse race, since the expected per turn winnings of the entropy game are in fact equal to the log (base 2) of the expected wealth increase rate per turn of a bettor at a horse racing track. Kelly-style bets on (say) percentiles of the original forecast distribution have an expected wealth growth rate whose logarithm is approximately equal to $\Delta S_{True}$. The entropy game is a useful alternative to the Kelly horse race because it is better adapted to continuous distributions and because its rules present a more symmetric appearance than does the bettor-and-bookie model of the horse race. It is our fiducial game for assessing the improvement in recalibrated forecasts, and we will make extensive use of it in what follows.

3.4 Predicting Recalibrated Forecast Performance

The entropy game winnings $\Delta S_{True}$ in Equation (19) are expressed in terms of the distribution $\pi(X,J|C)\propto\pi(X|J,C)$, which is known only approximately through the GPME fit. An important and useful feature of the GPME procedure is that it allows us to compute expected entropy game winnings per turn in advance of the game.

We first simplify $\Delta S_{True}$ by re-expressing it in PIT-space, using the fact that $\pi(X=x|J,C)dx=\pi(F=f|J,C)df$:

[TABLE]

Note that unlike the expression in Equation (9), $\Delta S_{True}$ is not non-negative definite but rather may in principle attain negative values if $\pi(F|{\cal F},C)$ is a poor estimate of $\pi(F|C)$.

We define the quantity

[TABLE]

which is a random variable (unlike $\Delta S_{True}$) in consequence of the use of the distribution-valued random variable $\Pi(F|C)$ for the average. This randomness expresses our epistemic uncertainty about the value of $\Delta S_{True}$ consequent on our uncertainty about the shape of the distribution $\pi(F|C)$. The probability distribution of $\Delta S$ is not directly computable by using GPME. However, we can calculate

[TABLE]

and

[TABLE]

These quantities represent the expectation and variance of $\Delta S$ with respect to the epistemic uncertainty contained in the GPME posterior distribution over the density $\Pi(F|C)\,|\,{\cal F}$. They constitute predictions of the outcome of rounds of the entropy game to be conducted out of sample with respect to the FOA training data that furnishes ${\cal F}$. Thus, not only can we verify the improvement in the recalibrated forecast using many rounds of the entropy game out of sample, but we can also predict in advance, with uncertainty bounds, what the average outcome of those games will be. The degree to which the predictions match the outcomes can be a useful gauge of model validity, as we will see below.

The expression for $\overline{\Delta S}$ is readily obtained (even without appealing to the GPME theory):

[TABLE]

That is, $\overline{\Delta S}$ is just $KL[\pi(F|{\cal F},C)\,||\,U(F)]$, where $U(F)$ is the uniform distribution on $[0,1]$. Unsurprisingly, we have $\overline{\Delta S}\geq 0$. Note that it is not the case in general that $\overline{\Delta S}=\Delta S_{True}$, since the weighted averages in Equations $($20) and (24) are different. The two quantities are “close” only to the extent that the distribution $\pi(F|{\cal F},C)$ approaches $\pi(F|C)$. The GPME model fit makes certain approximations (such as statistical independence of the data in ${\cal F}$, and the Laplace approximation for the log-density) and assumptions (such as kernel and mean function choice) that can in principle produce model errors that make $\overline{\Delta S}$ a flawed estimator for $\Delta S_{True}$. As we will see in the case studies below, GPME does a good modeling job in general, and such modeling errors can be acceptably small in practical cases.

Note also that despite the fact that $\overline{\Delta S}\geq 0$, it is not the case that the random variable $\Delta S\geq 0$, as can be seen from Equation (21). The random variable $\Delta S$ itself may certainly attain negative values for some realizations $\pi(F|C)$ of $\Pi(F|C)$, just as $\Delta S_{True}$ could in principle turn out to be negative if $\pi(F|{\cal F},C)$ badly mis-estimates $\pi(F|C)$. It would be reassuring to have some measure of how unlikely it is for $\Delta S<0$. This can be obtained from the variance of $\Delta S$.

We reproduce here the GPME expression for $\textrm{Var}(\Delta S)$, given in Equation (87):

[TABLE]

where $C(f_{1},f_{2})$ is the GPME posterior covariance over \ln\text{\Pi}(F|C). Using this formula, we may compute the uncertainty in the estimate $\overline{\Delta S}$ in terms of a straightforward two-dimensional quadrature on $[0,1]\times[0,1]$. In the Appendix, we show that in the asymptotic limit $N\rightarrow\infty$ of unlimited training data, $\textrm{Var}(\Delta S)\sim N^{-1}$. Thus the uncertainty in the estimate $\overline{\Delta S}$ scales asymptotically as ${\cal O}(N^{-1/2})$.

We may use this fact to introduce a forecast advantage measure (FAM):

[TABLE]

which asymptotically scales as ${\cal O}(N^{1/2})$. The FAM gives us a measure of a-priori confidence in the positivity of the out-of-sample entropy game winnings of the recalibrated forecast.

To summarize the results so far: Given an FOA and its associated PIT data ${\cal F}$, we have a recalibration procedure that allows us to improve current published forecasts $p(X;J,C$), replacing them with recalibrated forecasts $p_{1}(X;J,C)$. The recalibrated forecasts are expected to outperform the original published forecasts in ignorance score, and, equivalently, in entropy game performance. The extent of this superior performance—the average per round Entropy Game winnings $\Delta S_{True}$—may be estimated in advance using only the data in ${\cal F}$, and that estimate is attended by an uncertainty that is also computable using only the data in ${\cal F}$. Our confidence in the positivity of $\Delta S_{True}$—that is, in the superiority of the recalibrated forecast—can be expressed by the expression in Equation (26) for the FAM, which in the GPME theory grows with training dataset size $N$ as $N^{1/2}$. Hence we can reassure ourselves of the superior performance of the recalibrated forecast by accumulating a sufficiently large training set.

3.5 Fit Quality

As discussed at the end of §3.2, success of the density estimation procedure is crucial to the success of the recalibration procedure, since the latter success depends on the estimated distribution $\pi(F|{\cal F},C)$ not departing too much from the true distribution $\pi(F|C)$. The GPME theory furnishes a quantitative measure of the departure $\pi(F|{\cal F},C)$ from $\pi(F|C)$, through the quantity

[TABLE]

The quantity inside the expectation is the K-L divergence of the imperfectly known distribution $\Pi(F|C)$ from the estimated distribution $\pi(F|{\cal F},C)$. The GPME theory yields a closed-form expression for this fit quality measure. It is given in Equation (79), which we reproduce here:

[TABLE]

where $C(f_{1},f_{2})$ is the GPME posterior covariance over \ln\text{\Pi}(F|C).

Equation (28) expresses a very sensible result: the expected divergence between the estimated probability density and the true density is proportional to the average of the posterior variance weighted by the effective posterior probability density. As the quality of the fit improves, the variance decreases and takes $EI[\pi(F|{\cal F},C)]$ down with it. Thus, in some sense, $EI[\pi(F|{\cal F},C)]$ expresses fit quality. One must be cautious in interpretation, however, since no probabilistic interpretation (such as a “P-value”) attaches to $EI[\cdot]$, so it is difficult to say in an absolute sense how small a value of $EI[\cdot]$ is adequate. Moreover, relying on $EI[\cdot]$ for fit quality is, in effect, asking the model to report on its own success. The result will be conditioned by model assumptions (such as covariance kernel choice), and cannot directly detect the effects of model error on fit quality.

It is more useful to employ the fact, derived in the Appendix, that asymptotically, $EI[\pi(F|{\cal F},C)]\rightarrow B/2N$, where $B$ is the number of bins used in the GPME fit and $N$ is the number of datapoint in ${\cal F}$. The $N^{-1}$ scaling is checkable by varying the size of the training set ${\cal F}$. This behavior can be a useful diagnostic of model adequacy, as we will see below, since departures from this scaling at large $N$ may indicate model errors that are masked by noise at smaller values of $N$.

3.6 Thinning Data To Remove Correlations

The GPME theory has the considerable benefit that besides providing an efficient method to obtain the predictive distribution $\pi\left(F|{\cal F}\right)$, it also provides closed-form expressions for the desired entropy-related quantities $EI$, $\overline{\Delta S}$, $\textrm{Var}(\Delta S)$. It has one serious defect for our application, however: it assumes that a pointlike Poisson process governs the generation of the PIT values in ${\cal F}$. This assumption is more often false than true. In the case of weather, forecast cadences are generally more rapid than the characteristic times on which the dynamical system loses memory (hence the adage that the best predictor of tomorrow’s weather is today’s weather). This means that successive members of an FOA—both forecasts and observations—typically resemble each other more than do well-separated members. This effect manifests itself in correlations of successive PIT values in ${\cal F}$ – examples are displayed in the right panels of Figure 4 in §4.2. These correlations technically invalidate the point-like Poisson process assumption that underlies the GPME theory. Two bad consequences are that (1) the fit may be skewed, especially if the training set is not large; and (2) the a priori estimates of recalibrated forecast improvement over base forecast are not reliable, since they are based on an inaccurate statistical model.

To recover the utility of the GPME theory in such cases, one must thin the training dataset by a factor that may be inferred from the autocorrelation function of the PIT values in ${\cal F}$. This is a process analogous to the thinning of samples output by an MCMC chain [49, p. 149] and is necessary for the same reason: the samples obtained after appropriate thinning have good independence properties. They may therefore be appropriately modeled by a point Poisson process. The downside is that if data is not abundant, the thinning of the training set may be harmful to predictive performance.

One may, with some justice, ask what was the point of emphasizing the non-i.i.d. nature of the recalibration procedure, if an i.i.d. restriction is then re-introduced through the GPME fitting procedure. The answer is that GPME only imposes an i.i.d. restriction on the model training. The results in §3 on performance improvement of recalibrated forecasts are still valid for non-i.i.d. forecasts of future events, given an acceptable, statistically consistent regression estimate of $\pi\left(F|\mathcal{F},C\right)$ from the GPME procedure. There may quite possibly exist a generalization of GPME that takes proper account of non-i.i.d. behavior in $\mathcal{F}$. Locating such a procedure would be a promising avenue of future research, since the result would be a recalibration procedure that is entirely free of the i.i.d. restriction.

Note, however, that the GPME i.i.d. restriction on training data is only necessary to preserve the predictive performance of our recalibration procedure – that is, to be able to state in advance the expected improvement in forecast logarithmic skill. The restriction is not necessary to improve forecast skill by some (possibly difficult to predict) amount. As demonstrated by [28] and [29], several different styles of regression on the data $\mathcal{F}$ are capable of furnishing estimates of $\pi(F|\mathcal{F},C)$ that improve calibration. We can say that the present work advances the state of the art from the work of [28, 29] in that it is now clear that recalibration may be expected to work even in the case where the $F_{n}$ are not i.i.d., so long as some reasonable regression model for $\pi(F|\mathcal{F},C)$ is produced.

4 Verification

We now exhibit practical examples of the forecast recalibration procedure in two separate applications: a laboratory experiment with predictions of the output from a nonlinear circuit and a seasonal metereology example using ensemble forecasts of El Niño Southern Oscillations (ENSO) temperature fluctuations.

4.1 A Nonlinear Circuit

Our first application of the forecasting recalibration scheme is a laboratory experiment with predictions of the output from a nonlinear circuit. The circuit was first introduced in [50] and later discussed in more detail in [51]. References [51] and [52] discuss different aspects of its predictability properties. The circuit is designed to produce output voltages that mimic the Moore-Spiegel [53] three-dimensional system of ordinary differential equations:

[TABLE]

This system is a simplified model of a parcel of fluid moving vertically in a stratified fluid, with which it exchanges heat, while tethered by a harmonic force to a point [53]. The variable $z$ represents the height of the fluid element.

The circuit is set to operate at parameter values $R=10$, $\Gamma=3.6$, at which values the system exhibits chaotic behavior. Voltages $(V_{1},V_{2},V_{3})$ corresponding to the variables $(x,y,z)$ are measured at three points on the circuit. The ODE system of Equations (29) is scaled to endow its variables with the dimensions of voltage. The voltage $V_{3}$ corresponding to $z$ is our predictand. As noted in [51], the system of Equations (29) poorly predicts the behavior of the circuit due to model imperfection. An alternative prediction model is constructed using radial basis functions.

Probe voltages for the three voltage probes were collected over a duration of 14 hours at a sampling rate of 10 kH. A sample of 2000 points corresponding to the $z$ voltage probe was used to empirically estimate the climatological distribution $\rho(z)$ of $V_{3}$ (corresponding to $z$) (top-left panel of Figure 1) . Then 2,048 uncorrelated voltage states were sampled to furnish initial conditions from which forecasts could be initialized. Each of these states was used to create an ensemble of 127 forecasts by small Gaussian additive perturbations about the observed state and evolving the resulting states using a radial basis function model up to eight time steps ahead, that is, up to a forecast lead time of 0.8 ms. These 0.8 ms lead time forecasts were then converted to probabilistic forecasts for $V_{3}$ by kernel dressing and blending with climatology [23], wherein the forecast distribution density is expressed as a sum of kernels, each centered at the value of one of the 127 simulation values, and the result is linearly blended with the climatology $\rho(z)$. The kernels were chosen to be Gaussians, with equal widths chosen to minimize the ignorance score, and the linear blending parameter was also chosen to minimize the ignorance score [52].

The continuous ensemble forecasts thus generated were compared with the corresponding observed values of the $z$ voltage to obtain the PIT distribution shown in the top-right panel of Figure 1.

The 2,048 available observations and forecasts were divided into training and test sets, with training set sizes $N_{t}\in\{200,283,400,566,800,1131,1600\}$ (each about a factor of $\sqrt{2}$ larger than the previous value). In each case, the test set comprised all the remaining data. We carried out the recalibration procedure with each training set to compute the corresponding PIT posterior predictive density $\pi(F|{\cal F},C)$ and carried out entropy games over the corresponding test sets, recording the performance predictors ($EI[\pi(F|{\cal F},C)]$, $\overline{\Delta S}$, $\textrm{Var}(\Delta S)$, FAM) and the game outcomes.

The lower-left panel of Figure 1 shows the PIT fit from the training data (red line) superposed on the PIT histogram from the test data for the case $N_{t}=566$. The lower-right panel shows the PIT distribution of the recalibrated forecast for the same case. Comparing this figure with the one to its left we can see that the recalibration procedure successfully produced updated forecasts that are probabilistically calibrated.

The top-left panel of Figure 2 shows the run of FAM with $N_{t}$, displaying the expected $N_{t}^{1/2}$ trend. The top right panel of Figure 2 displays the run of $EI[\pi(F|{\cal F},C)]$ with $N_{t}$. The initial expected drop appears to level off at the highest values of $N_{t}$, possibly indicating some model inadequacy (for example, a poor choice of GP kernel) that reveals itself as the noise in the training histogram is suppressed by larger values of $N_{t}$.

The lower-left panel of Figure 2 displays Entropy Game winnings (red dots) together with the predicted winnings $\overline{\Delta S}$ (blue dots) and predicted uncertainty $\textrm{Var}(\Delta S)^{1/2}$ (error bars). Here again we see a tendency at the highest values of $N_{t}$ of the average winnings to depart from predictions—the actual winnings seem somewhat higher than predicted. Again, this discrepancy is possibly explainable in terms of inadequacies of the GP model used to estimate $\pi(f|{\cal F},C)$, which are perceptible only when the histogram noise abates at higher values of $N_{t}$. Nonetheless, the success of the model in predicting the entropy game winnings is gratifying, and the recalibrated model is clearly winning systematically against the ensemble forecasts.

The lower-right panel of Figure 2 shows a histogram of the outcomes of 1,482 rounds of entropy game for the case $N_{t}=566$, together with the empirical mean (green dashed line), predicted $\overline{\Delta S}$ (blue solid line) and predicted $1-\sigma$ interval (pink band), where $\sigma=\sqrt{\mathrm{Var}(\Delta S)}$ is computed using Equation (25). The distribution of outcomes is quite dispersed, comprising a sharp positive peak with a long tail of negative “bad busts.” The shape of the histogram is easily understood in terms of the GPME fit—the red line in the lower-left panel of Figure 1. In PIT space, this line approximates the actual PIT distribution, while the original forecast distribution is represented by a horizontal line at unit normalized frequency. The winnings cutoff at the right of the winnings histogram corresponds to the log of the maximum ratio between these two distributions, which coincides with the mode of the “actual” distribution. This is the reason that the mode of the winnings distribution is at the cutoff. By means of a second-order Taylor series at the mode of the “actual” PIT distribution, we can also show that the expected behavior in winnings space near the cutoff should approach $(w_{cutoff}-w)^{-1/2}$, which is (integrably) divergent. The long tail to the left is also intelligible, as it corresponds to the two tails near PIT values of 0 and 1 where the “actual” distribution is smallest. At these places, the fit distribution density has value of about $0.2$, so the tail should extend to values of $w$ near $w_{tail}=\log_{2}\left(0.2/1.0\right)=-2.3$. Note that these properties are to be expected of overdispersed base forecasts, because of the hump-shaped PIT distribution, and are not expected for, say, underdispersed predictions, where the PIT distribution looks like a pair of peaks near 0 and 1 with a valley in-between.

In Figure 3 we have displayed, for the case $N_{t}=1600$, the first nine ensemble forecasts (red dashed lines), recalibrated forecasts (blue solid lines), and observed $z$ voltages (green vertical lines). The plots show the overdispersion of the ensemble forecasts, manifest in the fact that the observations are too frequently near the median of the ensemble forecast. The recalibrated forecasts are sharper than the ensemble forecasts in this case. This would not be expected in general but is true here because of the overdispersion of the ensemble forecasts—the relative sharpness of the recalibrated forecasts restores the missing scatter in the PIT distribution. A noteworthy feature of this plot is that the recalibrated forecasts are less noisy than the ensemble forecasts, which are more prone to show their underlying discrete basis of ensemble-members that anchor the Gaussian mixture model of the continuous ensemble forecast. The transition from $p(x)$ to $p_{1}(x)$ appears to smooth out this noise somewhat.

In summary, the recalibration procedure is highly successful at improving the performance of the published ensemble forecasts of the nonlinear circuit. The average winnings of about 0.6 bits corresponds to a wealth amplification factor of $2^{0.6}=1.5$ per turn in a Kelly-style betting contest between the two forecasts—offering and wagering on odds on percentiles of the published ensemble forecast, say—which means that the ensemble forecaster would likely meet ruin in only a few rounds of betting.

4.2 El Niño Temperature Fluctuations

Our second application of the recalibration technique is to a seasonal forecasting problem. The seasonal forecast dataset used in this study is from the North American Multimodel Ensemble (NMME) project [54]. The NMME is a collection of global ensemble forecasts from coupled atmosphere-ocean models produced by operational and research centers in the United States and Canada. The NMME forecasts are generated in real time but also include a 30-year set of retropsective monthly forecasts (hindcasts) for assessing systematic biases in the models.

We examine the NMME model predictions of the intensity of El Niño-Southern Oscillation (ENSO) phenomenon. The ENSO state is often characterized by the NINO 3.4 index, which is the monthly mean sea surface temperature (SST) averaged over the equatorial Pacific region: 5S to 5N, and 170W to 120W. Tippet et al. [55] also use the NINO 3.4 index to assess the skill of the NMME models; that paper contains a useful description of all the models and their particular configurations for the hindcasts, and a discussion of the errors and known problems in the model forecasts. For this study, no data corrections have been made, and model climatologies have not been removed—the index values are based on real temperatures instead of temperature anomalies.

The NMME hindcast dataset is available at http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/. The NMME hindcasts are created monthly, have lead times that range from 1 to 12 months, and are validated with the observed NINO 3.4 index for the period January 1982–October 2017. The observed index values are derived from NOAA’s Optimum Interpolation Sea Surface Temperature data (OISST, version 2, [56]) which are available at http://www.cpc.ncep.noaa.gov/data/indices/sstoi.indices.

Our object in this study is to work with as long a stretch of data as possible and for that stretch of data to represent model output that is temporally as homogeneous as possible, because if the model composition were to fluctuate during the study, or be substantially different between training and test data sets, the recalibration procedure could not be expected to be effective. Of the 15 NMME models, only 6 were run daily during the entire period of the project, while others were retired at various stages of the project. Of those 6, 5 ran with the same ensemble size throughout, while the remaining model had an ensemble size that varied with sufficient frequency to create concern for the homogeneity of the sample. Consequently, we subsetted the hindcast data, choosing only the 5 models that were run consistently monthly for 35 years of hindcasts. These models and their respective ensemble sizes are displayed in Table 1.

To convert the forecast simulation ensembles to continuous forecast distributions, we chose the method of Bayesian model averaging (BMA) [24], as adapted to multi-model ensembles with exchangeable members by Fraley et al. [25]. Briefly, BMA models the forecast PDFs consequent from the ensemble predictions by using a mixture model, with mixture weights ascribed to different components. Ensemble members from a single model are “exchangeable” [25] in that none of them may be regarded as bearing better or worse information than their ensemble partners bear. They therefore are all assigned equal weights under the scheme. Ensemble members from different models are “nonexchangeable” and so have different weights. Following [25, 24], we use Gaussians for the mixture component PDFs, with Gaussian widths that are the same within each model ensemble but may differ from model to model. We first bias-correct the forecasts model by model using linear regression of the observations on the forecasts in the training set. Then we center the Gaussians on the bias-corrected ensemble forecast values and optimize the likelihood of the training data iteratively by the EM algorithm [57, 58], updating weights and Gaussian widths at each EM iteration [24, 25]. The converged weights and widths are used to create forecasts in the test set.

We have available a total of 430 monthly hindcast simulations. We consider lead times of 1–11 months for each hindcast. We train the BMA forecasting machinery on the first 36 hindcasts and use the machinery to create forecasts from the remaining 394 hindcasts, one for each lead time.

The ensemble forecast results are summarized in Figure 4. The three rows of the figure correspond to lead times of 4, 8, and 11 months. The left column depicts the PIT histograms, which show clear evidence that the BMA forecasts are biased to high values of NINO 3.4 index, despite the preliminary bias corrections. Clearly potential leverage exists here for the recalibration procedure to do its work. However, there is a fly in the ointment: the right column of Figure 4 displays the temporal autocorrelation functions of the PIT values, which are clearly significantly correlated out to 15-month lags and beyond. While not entirely surprising, this is a serious potential restraint on the effectiveness of the method, since with only 394 forecasts to work with, a thinning by a factor of 15 leaves hardly enough forecasts to form a training set, to say nothing of a test set.

We compromise, faute de mieux, on a thinning by a factor of 5—below this factor we find correlations unacceptably compromise the fit of $\pi(F|{\cal F},C)$, above it we have too few forecasts to work with. We choose a recalibration training set of 64 forecasts, leaving $394-5\times 64=74$ forecasts in the test set. For each lead time, we fit $\pi(F|{\cal F},C)$ to the corresponding PIT histogram and use it to compute $EI[\pi(F|{\cal F},C)]$, $\overline{\Delta S}$, $\textrm{Var}(\Delta S)$ and FAM, then to run 74 rounds of Entropy Game between the BMA forecasts and the recalibrated forecasts.

Figure 5 shows the result of the recalibration procedure on the PIT histograms for forecast leads of 4, 6, and 8 months. The 74 test forecasts are histogrammed into 20 bins. The recalibrated forecasts (right column) have improved probabilistic calibration over the published forecasts (left column). The effect is not as dramatic as for the circuit data of §4.1, in part because of the paucity of data and in part no doubt because of model inadequacy due to residual correlations in the data. Nonetheless the improvement in calibration is clear.

Figure 6 shows a few sample published and recalibrated forecasts for leads of 4, 6, and 8 months. Comparison with Figure 5 shows that the recalibration is attempting to correct the bias by shifting the mass of the probability distributions to lower values of the NINO 3.4 index.

In Figure 7, the top-left panel shows FAM as a function of forecast lead. FAM increases dramatically from lead=1 month to lead=2 months, then settles down between a value of 1 and 2, suggesting moderate confidence in the performance of the recalibrated forecast for months 2 and later. The top-right panel shows $EI\left[\pi(F|{\cal F},C)\right]$, which is fairly steady with a shallow peak near lead=7 months, which indicates that the quality of the fit of $\pi(F|{\cal F},C)$ to the PIT training data is fairly uniform across lead times.

The lower left panel shows entropy game average winnings (red dots) compared with performance measures $\overline{\Delta S}$ (blue dots) and $\mathrm{Var}(\Delta S)$ (blue errorbars). The performance measures do modestly well in predicting winnings, given the relatively modest size of the available training set and the correlations that still reside therein. The performance advantage of the recalibrated forecasts is striking, especially beginning around lead=4 months.

The lower-right panel shows a histogram of outcomes over 74 rounds of the entropy game in the 8-month forecast lead case. Again, the histogram structure is interpretable in terms of the PIT histogram fit in the middle-left panel of Figure 5, with the mode at $\log_{2}1.5\approx 0.6$ and a tail extending to $\log_{2}0.5=-1$.

The empirical average entropy game winnings are in the range 0.2–0.6 bits, corresponding to a range of per turn wealth multipliers of $2^{0.2}$–2^{0.6}$$=1.15–1.52 in a Kelly-style odds-making/betting game. Even despite the limitations set by the small amount of available training data, the recalibrated forecaster can expect to bankrupt the BMA forecaster after relatively few turns, especially if betting on leads of about 6 months or so.

5 Discussion

In this work we have gone to great length to emphasize the importance of weighing information—and of using quantitative measures of information—in reasoning about probabilistic forecasts of continuous variables. We showed that the information interpretation of calibration is useful because we may use it to build out of an arbitrary forecast system a related recalibrated system that is expected to be much better calibrated than the original system in those cases where the probabilistic calibration of the original system was noticeably poor.

We validated the forecast recalibration theory on two very different examples: (1) a nonlinear circuit whose output is forecast by iterating a radial basis function model constructed in delay space (See [51] for details) and a smoothing using kernel dressing and climatology blending, and (2) 30 years of monthly NINO 3.4 index observations, using forecasts generated from NWP ensembles smoothed by BMA. In each case, the nature of the observational data, the input elements to the forecast system, and the method of generating probabilistic forecasts were different. We emphasize that the recalibration procedure was successful in both cases, producing objectively superior forecasts (as measured by entropy game outcomes) without needing to care much about the inner nature of the forecasts that it improves upon.

In the ENSO study, the recalibrated forecasts were easily able to outperform the BMA forecasts, despite the modest training set size. This result is particularly striking in view of the fact that training was performed on simulations and data spanning about 27 years (after thinning), during which time some secular evolution of NINO 3.4 index dynamics certainly occurred due to carbon forcing, so that the test set comprising the remaining data necessarily represents somewhat different climatology from the training set. The success of the method suggests that while the climatology may evolve, the miscalibration of the forecast system may be more stable over time and hence may remain a reliable guide to recalibration.

We showed that the recalibrated forecasts have better ignorance scores than do the original published forecasts and can consistently win bets in the entropy game, a game that, while not fair (because the recalibrated forecast has more information than the original forecast, and hence the player that wields it has an edge), is not in any way biased toward one player or another by its rules. We also pointed out that while the entropy game is an abstract game, its expected winnings are directly related to the wealth multiplication factor of the player with the recalibrated forecast in Kelly-style odds-setting-and-betting games on outcomes such as percentiles of the original forecast.

For recalibration to work well, much depends on the power and flexibility of the modeling system used to fit the training set of PIT values; and for the performance of the recalibrated forecast to be predictable, the modeling system must give access to entropy measures of the distributions being estimated. The Gaussian process measure estimation scheme described in A, by providing a hyperparametric regression estimate of the PIT measure that yields estimates of Kullback-Leibler divergences from the true distribution, provides a highly satisfactory solution for this application.

This is far from saying that further development is unnecessary. In the first place, the studies presented in this work employed only the most basic and simple GP kernel—the squared-exponential—in modeling PIT distributions. In fact, there was evidence in §4.1, in the top-right panel of Figure 2, that at the largest training-set sizes the $EI[\pi(F|{\cal F},C)]$ plot deviates from the expected $N_{t}^{-1}$ behavior (see Equation 80), which could indicate a model defect that is masked by noise for smaller training sets. This sort of situation is probably not rare, so it would be worth investigating the effectiveness of more flexible covariances and possibly mixtures of such covariances.

Furthermore, the GPME model’s reliance on i.i.d. training data to create a regression model $\pi\left(F|\mathcal{F},C\right)$ qualifies the success of the recalibration procedure in removing the i.i.d. restrictions of the methods described in the DHT [28] and KFE [29] papers. Additionally, the necessity of thinning forecasts to create an approximately uncorrelated PIT training sample can be a daunting prospect in cases where the data is not sufficiently abundant to support adequate thinning. In the ENSO case an acceptable compromise was fortunately found. Nonetheless, thinning seems an undesirable nuisance imposed by a somewhat simplistic model—the log-Gaussian Cox process, which leads to simple closed-form expressions at the cost of requiring that the data be statistically independent. It would be interesting and useful to develop the modeling in a way as to account for correlations, instead of ignoring them, possibly by modeling the PIT distribution in two dimensions, PIT value and time, by using a two-dimensional Gaussian process. Other alternatives are certainly worth considering.

In this work, we have also argued from a perspective on forecast quality assessment that emphasizes the importance of decision support. Forecast skill scores can often be difficult to interpret in terms of decision support by someone wishing to ascertain the superiority of some forecast set over another; and since different choices of skill score do not agree on a unique sort order of forecast excellence, the process of preferring some forecast sets over others on the basis of skill has something of a beauty-contest air about it [37]. We have seen that one can rate the performance of forecast sets concretely, in terms of their ability to consistently win bets against other forecast sets. It would perhaps be well to emphasize this concrete interpretation of skill, since it would seem to translate more directly into actionable decision-making information.

Acknowledgements

The authors wish to thank both referees and this journal’s associate editor, for critiques that materially strengthened this article, and L. A. Smith for valuable discussions. This material was based upon work supported by the US Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, under Contract DE-AC02-06CH11347. Jennifer Adams acknowledges support from NSF (1338427), NOAA (NA14OAR4310160) and NASA (NNX14AM19G).

The submitted manuscript has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory (“Argonne”). Argonne, a U.S. Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-06CH11357. The U.S. Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan http://energy.gov/downloads/doe-public-access-plan.

Appendix A Gaussian Process Probability Measure Estimation

A large literature on probability density estimation exists, and a number of popular techniques including kernel density estimation (KDE) [59, Chapter 6], nearest-neighbor estimation [60, p. 257], and Gaussian mixture modeling [60, p. 259], as well as more sophisticated methods based on stochastic processes well adapted to density estimation, such as Dirichlet processes [61].

For this work, we choose a nonparametric approach based on Gaussian process (GP) modeling. GP modeling is a popular approach for modeling spatial, time series, and spatiotemporal data [62, 31] and has recently received a lucid introductory treatment in [63]. Some work applying GP modeling to density estimates from Poisson-process data has appeared recently [64]. We have selected and developed this technique because it is easy to implement and leads to readily computed closed-form expressions for the Shannon entropies that we require.

Our scheme builds on the same foundation as described in [64]: we model a log-Gaussian Cox process (LGCP), wherein a Gaussian process model is placed on the log-density of an inhomogeneous Poisson process, described further below. Our scheme has a few new features compared to the work described in [64]: We show how to approximately normalize the density distributions so that they are, in fact, approximately probability distributions; we exhibit closed-form expressions for Shannon entropies associated with fit uncertainties in the estimated densities; and we point out a curious—and to our knowledge, previously unrecognized—feature of this LGCP, which is that the Laplace method approximation to the Poisson likelihood yields a much better approximation to a Gaussian when the log-density is approximated than when the density is approximated directly.

A.1 Poisson Number Density Estimation

Suppose we have some i.i.d. sample points from some space. How do we estimate the distribution that gave rise to those points?

More precisely: Suppose we have an absolutely continuous finite measure $\mu$ over a set $\Gamma\subset\mathbf{R}^{D}$, representable by a density $\rho(\bm{x})$ so that. $d\mu(\bm{x})=\rho(\bm{x})d^{D}\bm{x}$. For any subset $b\subset\Gamma$, we interpret $\mu(b)$ as the mean of a Poisson distribution describing the number of events of some type that may occur in $b$. Evidently, $\mu(\Gamma)$ is the total expected number of events in $\Gamma$, and $\mu/\mu(\Gamma)$ is a probability measure over $\Gamma$, represented by a normalized probability density $\pi(\bm{x})=\rho(\bm{x})/\mu(\Gamma)$. We are given $N$ samples from this probability measure, denoted $\bm{x}_{k}$, $k=1,\ldots,N$. Suppose that $N\gg 1$. In this situation, which is common in many fields, one often would like a sensible way to estimate $\mu$, or, equivalently, $\rho(\bm{x})$.

The density $\rho(\bm{x})$ is known imperfectly, since it must be estimated from the data $\bm{x}_{k}$. Just as statistical estimation of a real-valued scalar quantity $v$ leads naturally to consideration of a real-valued scalar random variable $V$ whose possible realizations are values of $v$, the estimation of a density $\rho(\bm{x})$ from data leads naturally to consideration of a density-function-valued random variable $R(\bm{x})$ whose realizations are possible density functions $\rho(\bm{x})$. The simplest nontrivial theory of such stochastic functions is the theory of Gaussian processes [62, 63], which is what we exploit here.

We partition $\Gamma$ into bins, which will constitute measurable training sets. In order to avoid a coarse binning of the space that smears out spatial structure in $\rho(\bm{x})$, the bins should be geometrically small compared to the length scales in the measure. On the other hand, we would like to write down a likelihood for the data that somehow leverages the Gaussian nature of the GP model. But the Poisson likelihood, appropriate for this kind of data, is acceptably Gaussian in $\mu$ (the Poisson mean) only when $\mu$ is dispiritingly large, $\mu>15$ or so. Hence our requirement for small bins is, in principle, in tension with our requirement for populous bins.

Furthermore, a GP model of $\rho(\bm{x})$ will certainly not respect the positivity constraint $\rho(\bm{x})>0$. It might do so approximately, for certain choices of kernel function, in regions with abundant sample points, but we would like our model to apply correctly to sparsely sampled regions as well as crowded ones.

While these obstacles seem considerable, they are surmountable. The key element of the estimation procedure described below is log-Gaussian Cox process (LGCP), which models $\ln\rho$, rather than $\rho$ directly, so that $\rho$ automatically satisfies the positivity constraint. It turns out that by a stroke of good luck, choosing to model $\ln\rho$ as a GP (instead of $\rho$ directly) solves two problems at once. In the first place, such a model automatically satisfies the positivity constraint $\rho(\bm{x})>0$. More subtly, as a function of $\ln\rho$, the Laplace approximation to the Poisson likelihood approaches the Gaussian regime much more rapidly than it does as a function of $\rho$!

This is not difficult to demonstrate. Consider a single Poisson variate $n$ with mean $\mu=\exp(l)$. The Poisson likelihood for an observation of $n$ is $\pi=e^{-\mu}\mu^{n}$. The normal approximation is achieved by expanding $\ln\pi$ in a Taylor series about its maximum. As a function of $\mu$, this is to say

[TABLE]

where $\bar{\mu}=n$, $\sigma_{\mu}^{\,2}=n$, and where we denote the approximation residual by $R_{1}(n,\mu)$. As a function of $l$, we may similarly expand around the maximum, obtaining

[TABLE]

with $\bar{l}=\ln n$ and $\sigma_{l}^{\,2}=1/n$. The functions $R_{1}$ and $R_{2}$ are defined implicitly by Equations (30–31) and may be computed directly from those formulae. The result is plotted in Figure 8.

The figure shows values of the residuals computed at their respective $\pm\sigma$ point about their respective means, that is, $R_{1}(n,\bar{\mu}\pm\sigma_{\mu})$ and $R_{2}(n,\bar{l}\pm\sigma_{l})$. They show that at these characteristic values of the normal distribution to be approximated, the residual $R_{2}$ is considerably smaller than the corresponding $R_{1}$, especially at the $-\sigma$ points. In fact, the accuracy attained by the Gaussian approximation in $\mu$ at $n=14$ is exceeded at $n=3$ by the accuracy of the approximation in $l$. By $n=10$ or so, the accuracy of the normal approximation to $l$ exceeds that attained by the approximation to $\mu$ at $n>40$. This is reassuring because it suggests that the Laplace approximation error can be well-controlled even for small bin counts of order 5–10.

Now suppose the space $\Gamma$ has volume $\int_{\Gamma}d^{D}\bm{x}=\Omega$. We break up the space into $B$ disjoint bins $b_{\nu},\nu=1,\ldots,B$, satisfying $\bigcup_{\nu}b_{\nu}\subset\Gamma$, $b_{\nu}\bigcap b_{\beta}=\emptyset$ if $\nu\neq\beta$, $\int_{b_{\nu}}d^{D}\bm{x}\equiv v_{\nu}$. We associate each bin $b_{\nu}$ with a coordinate label $\bm{x}_{\nu}$, which is usually the location of the center of the bin. In what follows, we will assume that the bins are small in the sense that the density $\rho(\bm{x})$ is approximately constant over each bin.

Since the density is approximately constant over each bin, we may set $\int_{b_{\nu}}d^{D}\bm{x}\,\rho(\bm{x})=\rho(\bm{x}_{\nu})v_{\nu}.$ We want to model the log of $\rho$ rather than $\rho$ directly. To do so, we need a reference volume scale to nondimensionalize $\rho$. The volume $\Omega$ will do nicely. We therefore set

[TABLE]

from which follows

[TABLE]

defining the dimensionless volume element $\omega_{\nu}\equiv v_{\nu}/\Omega$.

Suppose we observe $n_{\nu}$ samples from the density in bin $b_{\nu}$. If we are given the density $\rho(\bm{x})$, we may compute the Poisson process likelihood $\mathcal{L}\left(\bm{X}|\bm{l}\right)$ of the data. Setting $\rho(\bm{x}_{\nu})\equiv\rho_{\nu}$ and $\bm{X}=\{\bm{x}_{k},k=1,\ldots,N\}$ for brevity, we have

[TABLE]

where we have appealed to the normal approximation discussed above, and we defined

[TABLE]

Clearly at this point that we have committed to having no empty bins, since any empty bin completely compromises the normal approximation that we have just introduced. In fact, as discussed above, we will be happy with bin sample counts in the 5–10 region.

Note that we have implicitly assumed here that the data are well described by a Poisson point process and, in particular, that the $\bm{x}_{k}$ are i.i.d. If this assumption is incorrect, then Equation (34) is not the correct expression for the likelihood. In some cases, when the $\bm{x}_{k}$ arise from a stationary time-series with nonzero correlations $<\bm{x}_{k}\bm{x}_{l}>=f(|k-l|)$, we may be able to thin out the data by a factor suggested by the shape of the autocorrelation function, so as to obtain an approximately uncorrelated data sample that may be correctly modeled by using Equation (34).

We will model the function $l(\bm{x})$ using a constant-mean Gaussian process,

[TABLE]

where the mean function $l_{0}(\bm{x})$ is in fact a constant $l_{0}$ and the covariance $K(\bm{x},\bm{x}^{\prime};\theta)$ is a positive-definite (as an integral kernel) function, parametrized by some hyperparameters denoted by $\theta$. For example, we might choose a stationary kernel with a scale hyperparameter $\sigma$ and an amplitude hyperparameter $A$, that is,

[TABLE]

in which case $\theta=(A,\sigma)$. The specific form of the covariance function is not needed here, and it could in general be chosen from the many known valid covariance forms, as seems suited to the type of measure being modeled [63, Chapter 4].

The choice of a parametrized mean level $l_{0}$ is made here because a zero-mean GP (the more usual choice) effectively makes a choice of amplitude scale for $\rho$ that is not selected by the data. One commonly obviates this kind of issue by mean-subtracting the data $\bm{l}_{1}$. However, many weighted means of the data in $\bm{l}_{1}$ could be chosen for this purpose, and it is not a priori clear that the usual unweighted mean $\bar{l}=B^{-1}\sum_{\nu=1}^{B}l_{\nu}$ is optimal among these. Setting the mean level as an adjustable parameter selects a certain weighted mean as the maximum-likelihood estimate (MLE) of $l_{0}$ and turns out to be computationally inexpensive, as we will see below.

The covariance matrix $\bm{Q}$ arising from the GP model is just the Gram matrix of the covariance function,

[TABLE]

where the indices $\nu,\nu^{\prime}$ range over the $B$ bins. The mean vector arising from the process is $\bar{\bm{l}}=l_{0}\bm{u}_{B}$, where $\bm{u}_{B}$ is the $B$-dimensional “one” vector, $\left[\bm{u}_{B}\right]_{\nu}=1$, $\nu=1,\ldots,B$, and $l_{0}$ is a parameter to be estimated, residing in the mean function rather than in the covariance kernel.

In terms of $\bm{Q}$ and $\bar{\bm{l}}$, the probability of a density $\rho$ represented by $B$-dimensional vector $\bm{l}$ is

[TABLE]

where we have symbolically collected in $I$ all conditioning information such as parameters $\theta,l_{0}$ and covariance kernel choices.

The Poisson process likelihood of Equation (34) can be marginalized over the GP distribution for $\rho$ of Equation (41) to produce the marginal likelihood:

[TABLE]

This is, unsurprisingly, the form for the marginal likelihood of a GP trained on noisy data $\bm{l}_{1}$ with noise covariance $\bm{D}$. We may therefore take over the standard result for the marginal likelihood [63, Eq. 2.30], adapted for the case of non-constant mean and heteroskedastic noise:

[TABLE]

where

[TABLE]

We may use Equation (43) to obtain an MLE of $l_{0}$ as a function of the kernel parameters $\theta$. We define the function $\tilde{l}_{0}(\theta)$ by

[TABLE]

Then $l_{0}=\tilde{l}_{0}(\theta)$ is the conditional MLE of $l_{0}$ given $\theta$. Defining $\bm{g\equiv}\left(\bm{Q}+\bm{D}\right)^{-1}\bm{u}_{B}$, we may write

[TABLE]

which shows that the MLE of $l_{0}$ is in fact a weighted average of $\bm{l}_{1}$, as asserted previously.

Combining Equations (45) and (44), we obtain

[TABLE]

which, by the Schwartz inequality, is non-negative definite and can attain a zero value only when $\bm{l}_{1}=\alpha\bm{u}_{B}$ for some scalar $\alpha$.

Combining Equation(47) with the negative log of Equation (43), we conclude that the MLE of $\theta,l_{0}$ are obtained by minimizing the objective function $S(\theta)$:

[TABLE]

[TABLE]

Suppose that we would like to estimate the density $\rho(\bm{y}_{a})$—or more to the point, the log-density $l_{a}^{(pred)}\equiv\ln\left(\rho(\bm{y}_{a})\Omega\right)$—at a set of points $\bm{y}_{a},$ $a=1,\ldots,P$. We may do this by following the standard methodological path of GP regression. We first extend the covariance matrix to an $(B+P)\times(B+P)$ matrix $\hat{\bm{Q}}$, writing

[TABLE]

with

[TABLE]

We further define $\bm{l}^{(pred)}=[l_{1}^{(pred)},\ldots,l_{P}^{(pred)}]^{T}$, and $P$-dimensional “one” vector $[\bm{u}_{P}]_{a}=1$, $a=1,\ldots,P$. Then we may take over the standard formula for predictions by a GP trained with noisy data [63, pp. 16–18], again adapted for nonzero mean and heteroskedastic noise. That is, $\bm{l}^{(pred)}\sim\mathcal{N}(\bm{\lambda}^{(pred)},\bm{C}^{(pred)})$, with

[TABLE]

and

[TABLE]

This is essentially a new, updated Gaussian process, with “trained” mean function

[TABLE]

and “trained” covariance function

[TABLE]

Returning to the higher-level “random function” view, we may summarize the story so far as follows. There is an unknown unnormalized density $\rho(\bm{x})$, from which a set of points $\bm{X}=\{\bm{x}_{k},k=1,\ldots,N\}$ is iid sampled. Since $\rho(\bm{x})$ is imperfectly known, we represent it by a function-valued random variable $R(\bm{x})$ and its scaled logarithm $l(\bm{x})=\ln(\Omega\rho(\bm{x}))$ by a function-valued random variable $L(\bm{x})$ to be estimated by GPME. Realizations of $L(\cdot)$ are possible log-density functions $l(\cdot)$. Similarly, realizations of the function-valued random variable $R(\cdot)=\Omega^{-1}\exp(L(\cdot))$ are possible density functions $\rho(\cdot)$.

The prior distribution over $L(\cdot)$ is a hierarchical model featuring a Gaussian process with constant mean function $l_{0}$ and a covariance function $K(\bm{x},\bm{y};\theta)$, as well as some prior distribution over $l_{0},\theta$ that we will not need to specify since we will proceed by maximizing the likelihood with respect to these parameters (when the parameter priors change slowly compared with the likelihood, this is approximately MAP estimation). We represent this prior distribution by the notation $L(\cdot)|I\sim GP[l_{0},K(\cdot,\cdot)]$. Training with the data $\bm{X}$ yields an updated posterior distribution for $L(\cdot)\,|\,(\bm{X},I)$, which is a Gaussian process with mean function $\lambda(\bm{x})$ and covariance function $C(\bm{x},\bm{y})$. Notationally, $L(\cdot)\,|\,(\bm{X},I)\sim GP\left[\lambda(\cdot),C(\cdot,\cdot)\right]$.

So far, we have modeled the imperfectly known log-density $L(\bm{x})$, rather than $R(\bm{x})$. This choice has consequences for the inferred Poisson process density, which is not, as one might naively assume, simply a constant times $\exp(\lambda(\bm{x}))$. In a small volume $v\equiv\Omega\omega$ about a location $\bm{x}$, the expected number of events $n$* given* the imperfectly known log-density $L(\bm{x})$ is

[TABLE]

Given the GP posterior predictive distribution $L(\cdot)\,|\,\bm{X}\sim GP[\lambda(\cdot),C(\cdot,\cdot)]$ , the effective number density $\rho_{E}(\bm{x})$ at $\bm{x}$ is given by

[TABLE]

We see that the log expected number of events is shifted with respect to the log-density $l(\bm{x})$ by the nonconstant factor $C(\bm{x},\bm{x})/2$.

A.2 Probability Density Estimation

The posterior predictive probability density $\pi(\bm{x}|\bm{X},I)$ that a future event should occur within a differential volume $d^{D}\bm{x}$ of $\bm{x}$ is the normalized version of $\rho_{E}(\bm{x})$:

[TABLE]

where the notation $\pi(\bm{x}|\bm{X},I)$ will be justified below. In the asymptotic limit $N\rightarrow\infty$, the normalization constant may be directly estimated from the data. We replace the integral by a sum over the training bins, in effect selecting the same prediction points as training bin centers:

[TABLE]

where

[TABLE]

and

[TABLE]

In the asymptotic limit, $\bm{D}\rightarrow 0$, and we have $\bm{\lambda}\approx\bm{l}_{1}$, $\bm{C}\approx\bm{0}$. Since $\left[\bm{l}_{1}\right]_{\nu}=\ln\frac{n_{\nu}}{\omega_{\nu}}$, we have

[TABLE]

which is an entirely unsurprising result.

Combining Equations (59), (60), and (61), we may write

[TABLE]

We have seen that GPME furnishes a tractable posterior distribution over $R(\bm{x})$, the imperfectly known Poisson number density that estimates $\rho(\bm{x})$. Now consider the normalized probability density $\pi(\bm{x})=\rho(\bm{x})/J[\rho]$, where $J[\rho]\equiv\int d^{D}\bm{x}\,\rho(\bm{x})$. We would like to use GPME to obtain a posterior distribution over $\Pi(\bm{x})$, the imperfectly known probability density (a probability-density-valued random variable) that estimates $\pi(\bm{x})$.

An unfortunate property of the the GPME scheme is that while it yields a tractable distribution over log number density $L(\bm{x})$, it does not yield a tractable distribution over probability density $\ln\Pi(\bm{x})$. The reason is that the transformation from $L(\bm{x})$ to $\ln\Pi(\bm{x})$ is nonlinear and does not transform the Gaussian distribution in $L$ into another tractable distribution over either $\Pi$ or $\ln\Pi$. The normalization factor $A$ discussed above pertains to the effective number density $\rho_{E}(\bm{x})$, which, according to the first line of Equation (59), is the expectation of the density $R(\bm{x})$ over the posterior distribution of the GP. This factor allows us to transition from the expected density to the posterior predictive probability density $\pi(\bm{x}|\bm{X},I)$. The factor $A$ is not, in general, the normalization appropriate to the imperfectly known density $R(\bm{x})$.

If we are satisfied with approximate normalization, however, then the factor $A$ is an appropriate normalization. The reason is that, as we now show, in the asymptotic limit $E_{L(\cdot)\,|\,\bm{X}}\left\{J\right\}=A\approx N$, $\mathrm{Var}_{L(\cdot)\,|\,\bm{X}}\left\{J\right\}\approx N$ so that $E_{L(\cdot)\,|\,\bm{X}}\left\{J\right\}/\sqrt{\mathrm{Var}_{L(\cdot)\,|\,\bm{X}}\left\{J\right\}}\approx N^{-1/2}$. Consequently, for large $N$, only a small error is committed by replacing $J[\rho]$ by $A$, and we may set

[TABLE]

This amounts to a constant offset of $\ln\Pi$ from $L$, so that the Gaussian distribution over $L(\bm{x})$ is simply mean-shifted by $-\ln(A\Omega)$ to produce the Gaussian distribution over $\ln\Pi$.

To show the required expectations, we again approximate the integral $J$, a random variable, by the sum over observed bins,

[TABLE]

so that

[TABLE]

Furthermore,

[TABLE]

Defining the $B$-dimensional vector $\bm{m}$ by

[TABLE]

we may write this as

[TABLE]

We then have

[TABLE]

In the asymptotic limit, by Equation (63) $\bm{C}\rightarrow\bm{D}$, which is diagonal and has small matrix elements $1/n_{\nu}$, so that

[TABLE]

Equations (68) and (73) are the required relations that allow us to approximate $\Pi(\bm{x})\approx A^{-1}L(\bm{x})$ in the asymptotic regime and hence approximate the posterior distribution over $\ln\Pi$ by a Gaussian process. In this light, we may cast the posterior predictive distribution $\pi(\bm{x}|\bm{X},I)$, defined in Equation (60), as

[TABLE]

This equation provides the justification for attaching the notation $\pi(\bm{x}|\bm{X},I)$ to the posterior predictive distribution.

A.3 Entropy Estimation

The practical output of the probability density estimation procedure is the posterior predictive probability density $\pi(\bm{x}|\bm{X},I)$. One might ask how far this is from the imperfectly known true probability density $\Pi(\bm{x})$. The Kullback-Leibler divergence between the two distributions is

[TABLE]

which is a random variable that measures departure of $\Pi(\bm{x})$ from $\pi(\bm{x}|\bm{X},I)$. We may calculate the expected divergence

[TABLE]

We have

[TABLE]

and

[TABLE]

Putting all this together, we find

[TABLE]

This is an intuitively reasonable result: the expected divergence between the estimated probability density and the true density is proportional the average of the posterior variance weighted by the posterior predictive density. As the quality of the fit improves, the variance decreases and takes $ES\left[\pi(\bm{x}|\bm{X},I)\right]$ down with it.

We may obtain an asymptotic estimate of $ES\left[\pi(\bm{x}|\bm{X},I)\right]$ by using the training bins as prediction points and approximating the integral by its finite Riemann sum, as we did above. Then,

[TABLE]

This asymptotic behavior is reassuring since its simple dependence on the average number of events per bin is in accordance with intution. Of course, how rapidly the asymptotic result becomes a reasonable approximation depends on how quickly the first term in Equation (63) eclipses the second term as $N\rightarrow\infty$, which is to say, on the choice of covariance kernel function $K(\bm{x}_{1},\bm{x}_{2})$, on the best-fit hyperparameters and, therefore, ultimately on the data and the distribution that gave rise to it.

Suppose someone has proposed a different probability density $p(\bm{x})$ as the source of the data. Can we tell whether $\pi(\bm{x}|\bm{X},I)$ is an improvement on $p(\bm{x})$?

Define

[TABLE]

The quantity $\Delta S$ is a random variable, in consequence of the uncertainty in the imperfectly known distribution $\Pi(\bm{x})$. Taking the expectation value of Equation (81), we find

[TABLE]

By the properties of the entropy, we have $\overline{\Delta S}\geq 0$, with equality holding only if $p(\bm{x})=\pi(\bm{x}|\bm{X},I)$ almost everywhere. This is not to say that $\pi(\bm{x}|\bm{X},I)$ is always superior to any other distribution $p(\bm{x})$, however (what if $p$ were, in fact, the ideal distribution $\pi(\bm{x}|I)$?). The quantity $\Delta S$ is uncertain and may in fact be negative; $\overline{\Delta S}$ is merely its expected value. The distribution for $\Delta S$ is too difficult to compute, but we may compute its variance straightforwardly:

[TABLE]

But

[TABLE]

where $\bm{\lambda}_{2}^{T}=[\lambda(\bm{x}_{1}),\lambda(\bm{x}_{2})]$, $\left[\bm{C}_{2}\right]_{ij}=C(\bm{x}_{i},\bm{x}_{j})$ with $i,j=1,2$, and $\bm{u}_{2}$ is the two-dimensional “one” vector. In other words,

[TABLE]

Substituting in Equation (83), we have

[TABLE]

It follows immediately that

[TABLE]

The asymptotic approximations for $\overline{\Delta S}$ and $\mathrm{Var}(\Delta S)$ are

[TABLE]

Since in general $n_{\nu}/N$ tends to a finite value in the limit, we see that $\mathrm{Var}(\Delta S)\sim\mathcal{O}(N^{-1})$ and tends to zero in the limit. On the other hand, if $p(\bm{x})$ is misspecified, one expects $\lim_{N\rightarrow\infty}\overline{\Delta S}$ to be a finite positive value. It follows that the GP measure estimate $\pi(\bm{x}|\bm{X},I)$ can achieve significant performance improvement over a misspecified distribution $p(\bm{x})$ in the limit of large data and that we would expect to be able to exploit this superior performance (in betting against the owner of $p(\bm{x})$, say), once $N$ is large enough that $\overline{\Delta S}/\sqrt{\mathrm{Var}(S)}\gg 1$.

If $p(\bm{x})$ is not misspecified—if it happens to be the true distribution $\pi(\bm{x})$—then we can set $n_{\nu}/N=\pi(\bm{x}_{\nu})(1+\epsilon)$, where $\epsilon\sim n_{\nu}^{-1/2}$ , in the limit of large $N$. We therefore have $\overline{\Delta S}\sim\mathcal{O}(N^{1/2})$ in this limit, so that asymptotically $\overline{\Delta S}/\sqrt{\mathrm{Var}(S)}$ tends to a constant.