Appraisal of Statistical Predictability under Uncertain Inputs: SST to Rainfall

In stochastic approximation, $p(\Theta |\mathcal{D})$ is numerically approximated by using sampling algorithms. Markov Chain Monte Carlo (MCMC) methods, developed in the physics literature, are the most popular sampling method used in Bayesian inference because they are efficient in sampling high dimensional spaces (large $p$), easy to design, and applicable to large class of distributions (Bishop 2006). These methods sample by constructing a Markov chain whose stationary distribution (equilibrium distribution) is the required posterior distribution. Two MCMC methods used in this paper, namely the Metropolis-Hasting algorithm and Gibbs sampler, are briefly outlined.

The Metropolis-Hasting (MH) algorithm was developed by Metropolis et al. (1953) and later extended by Hastings (1970). The algorithm draws samples from $p(\Theta |\mathcal{D})$ by generating a first order Markov chain in which each state ${\Theta }^t$ depends only on its previous state ${\Theta }^{t-1}$ through a jumping distribution $J({\Theta }^t|{\Theta }^{t-1})$, also known as proposal or candidate generating distribution, as follows:

After sufficient burn-in period (say $\tau$ steps), the Markov chain approaches its stationary distribution and the samples generated $({\Theta }^{\tau +1},\dots ,{\Theta }^{\tau +N_s})$ are from $p(\Theta |\mathcal{D})$.

Gibbs sampling (Geman and Geman 1984) is a special case of MH algorithm in which $r=1$, i.e., candidate samples are always selected. The method is used when it is easy to draw samples from the univariate conditional distribution $p({\theta }_i|{\mathit{\theta }}_{\setminus i})$, where ${\theta }_i$ is the $i$th unknown variable and ${\mathit{\theta }}_{\setminus i}$ denotes ${\theta }_1,\dots ,{\theta }_p$ but with ${\theta }_i$ omitted. The following steps are involved:After sufficient burn-in-period, the samples generated will be from the required posterior distribution $p(\Theta |\mathcal{D})$.

MCMC methods, owing to increased computational power of modern computers, have revolutionized the field of Bayesian inference. These methods are guaranteed to converge given infinite processing time. However, because processing time is finite, various diagnostics have been developed to ascertain convergence. The most commonly used tools are trace-plots of the generated samples, Geweke test (Geweke 1992), Raftery-Lewis test (Raftery and Lewis 1992), Gelman and Rubin test (Gelman and Rubin 1992), and diagnostics on the basis of autocorrelation of the generated samples.

The aim in this study is to approximate $p(\Theta |\mathcal{D})$ by $Q(\Theta |\mathcal{D})$ such that $Q(\Theta |\mathcal{D})$ is similar to $p(\Theta |\mathcal{D})$ and at the same time has a closed form. The methods are based on the principles of calculus of variations and are therefore referred to as variational inference or variational Bayes (Bishop 2006) methods. In the machine learning literature, they are also known as ensemble methods (MacKay 2005). The variational method used in this study has its origin in mean-field theory (Parisi 1988) where it was developed to approximate the lowest-energy state in quantum mechanical calculations. In this method, the approximating distribution $Q(\Theta |\mathcal{D})$ is assumed to factorize asIn the remainder of the paper, to make the notation uncluttered, explicit dependence of approximating distribution on data will be dropped and the $i$th component of the distribution will be denoted by $q({\theta }_i)$.

To estimate $Q(\Theta )$, the Kullback-Leibler (KL) divergence (Kullback and Leibler 1951) between $Q(\Theta )$ and $p(\Theta |\mathcal{D})$ [Eq. (8)], which is a measure of distance between the two distributions, is minimizedThe KL-divergence can take only positive values with minimum value of zero only when $Q(\Theta )=p(\Theta |\mathcal{D})$. Substituting $Q(\Theta )$ from Eq. (3) in Eq. (4) and minimizing KL-divergence with respect to distributions $q({\theta }_i)$ yields the following estimates:where $q^{\star }({\theta }_i)$ = optimal value of the approximating distribution; $p(\mathcal{D},\Theta )$ = joint distribution over data and the unknown variables; and ${\mathbb{E}}_{j\ne i}(·)$ denotes expectation with respect to variables ${\theta }_1,\dots ,{\theta }_p$, except ${\theta }_i$. Eq. (5) provides the optimal solution for $Q(\Theta )$ subject to the factorization constraint [Eq. (3)]. However, the solution is not explicit because $q^{\star }({\theta }_i)$ depends on the expected value with respect to variables ${\mathit{\theta }}_{\setminus i}$ and should be estimated iteratively. Nevertheless, convergence is guaranteed because KL-divergence is convex with respect to each factor $q({\theta }_i)$.

An additional benefit of variational methods is that they can be used for model selection, i.e., selecting the best model from a set of candidate models. The candidate that corresponds to the minimum value of KL-divergence is the best model. Because KL-divergence cannot be estimated directly as it is a function of $p(\Theta |\mathcal{D})$, the following decomposition of the model evidence is utilized:where $\mathcal{L}(Q)=\int Q(\Theta )\ln [p(\mathcal{D},\Theta )/Q(\Theta )]d\Theta$. Since $p(\mathcal{D})$ is constant, minimizing KL-divergence is equivalent to maximizing $\mathcal{L}(Q)$. The quantity $\mathcal{L}(Q)$ is known as lower bound and can be easily estimated for given $Q(\Theta )$. The best model is the one that maximizes the lower bound. A secondary advantage of calculating $\mathcal{L}(Q)$ is realized in monitoring the convergence of the iterative estimation procedure and checking the software implementation by noting that each iteration should always increase the value of $\mathcal{L}(Q)$.

Let a univariate independent variable $x$ and dependent variable $y$ follow bivariate Gaussian distribution with mean $\mathit{\mu }=[\begin{array}{cc}{\mu }_x& {\mu }_y\end{array}{]}^T$ and covariance $\mathbf{\Sigma }$, i.e., $p(x,y|\mathit{\mu },\mathbf{\Sigma })\sim \mathcal{N}([\begin{array}{cc}x& y\end{array}{]}^T|\mathit{\mu },\mathbf{\Sigma })$. Further, consider $N$ i.i.d. samples from this distribution $[\begin{array}{cccc}(x_1,y_1),& (x_2,y_2),& \dots ,& (x_N,y_N)\end{array}]$ measured with heterogeneous errorswhere $[\begin{array}{cc}{\tilde{x}}_n& {\tilde{y}}_n\end{array}{]}^T$ = $n$th observed sample; and $[\begin{array}{cc}{\zeta }_n& {\xi }_n\end{array}{]}^T$ = associated measurement error vector. The measurement errors are assumed to be independent and follow zero-mean Gaussian distributionwhere $b_n^{-1}$ and $h_n^{-1}$ denote the variance of the measurement error in $x_n$ and $y_n$, respectively. They are also referred to as precision of the individual measurements. In BaNCorr, the true values ($x_n$ and $y_n$) are treated as latent variables. The covariance matrix is decomposed asIn Eq. (10), $\mathit{s}=[\begin{array}{cc}{\sigma }_x& {\sigma }_y\end{array}{]}^T$ = vector of standard deviations; $\mathrm{diag}(\mathit{s})$ = diagonal matrix with diagonal elements $\mathit{s}$; and $\mathit{R}$ = correlation matrix given byTo estimate correlation between $x$ and $y$, ${\rho }_{xy}$, prior distributions are assigned to all unknown parameters, which for the most general (also most typical) case include $\mathit{\mu }$, $\mathit{s}$, and ${\rho }_{xy}$. The posterior distribution is estimated by using Bayes rule as follows:where the likelihood of the data $p(\tilde{\mathit{x}},\tilde{\mathit{y}}|\mathit{\mu },\mathbf{\Sigma })$ is given byIn Eq. (13), owing to the assumption of zero-mean Gaussian measurement errors, the distributions $p({\tilde{x}}_n|x_n,b_n^{-1})$ and $p({\tilde{y}}_n|y_n,h_n^{-1})$ also follow Gaussian distribution with mean $x_n$ and $y_n$, and variance $b_n^{-1}$ and $h_n^{-1}$, respectively.

The prior distribution in Eq. (12), $p(\mathit{\mu },\mathit{s},{\rho }_{xy})$, can be selected in various ways. In BaNCorr, the selection was guided by a desire to have mathematical elegance and simplicity in incorporating prior information into the learning process. The priors for $\mathit{\mu }$, $\mathit{s}$, and ${\rho }_{xy}$ are assumed to be independent, i.e., $p(\mathit{\mu },\mathit{s},{\rho }_{xy})=p(\mathit{\mu })p(\mathit{s})p({\rho }_{xy})$. For typical hydro-climatic variables, it is difficult to elicit reliable prior information on how mean, variance, and correlation are related. The prior distributions for ${\mu }_x$ and ${\mu }_y$ are selected to be Gaussian distributions, i.e., $p({\mu }_x)\sim \mathcal{N}({\mu }_x|{\nu }_{x0},{\beta }_{x0}^{-1})$ and $p({\mu }_y)\sim \mathcal{N}({\mu }_y|{\nu }_{y0},{\beta }_{y0}^{-1})$; Gamma distributions are selected for precision (inverse of the variance) of the individual variables, i.e., $p({\sigma }_x^{-2})\sim \mathrm{Gamma}({\sigma }_x^{-2}|c_{x0},d_{x0})$, $p({\sigma }_y^{-2})\sim \mathrm{Gamma}({\sigma }_y^{-2}|c_{y0},d_{y0})$, and $p(\mathit{s})=p({\sigma }_x)p({\sigma }_y)$; and the prior for ${\rho }_{xy}$ is set to be uniform distribution, $p({\rho }_{xy})\sim U({\rho }_{xy}|{\rho }_L,{\rho }_U)$. Prior knowledge of the data can be incorporated into BaNCorr by adjusting the hyperparameters in the prior distribution. If no a priori knowledge is available, the prior distribution can be made noninformative (broad) by setting: ${\nu }_{x0}={\nu }_{y0}=0$; ${\beta }_{x0}={\beta }_{y0}={10}^{-3}$; $c_{x0}=d_{x0}=c_{y0}=d_{y0}={10}^{-3}$; and ${\rho }_L=-1$, ${\rho }_U=1$. Fig. 1 shows the schematic of the BaNCorr model by using directed acyclic graph (DAG). The notational convention used in a DAG is given in the appendix.

The posterior distribution for the model parameters, $\mathit{\mu }$, $\mathit{s}$, and ${\rho }_{xy}$, and latent variables, $x_n$ and $y_n$, cannot be estimated analytically, and therefore they were approximated by using the Gibbs sampler. The prior distribution over latent variables is implicitly decided on the basis of the prior distribution of model parameters as shown in Fig. 1. The Gibbs sampling procedure requires distribution of each unknown variable conditioned on the values of all remaining variables. For the BaNCorr model, the required distributions are given by (to make notations uncluttered, explicit dependence of the distributions on the observed data $\tilde{x}$ and $\tilde{y}$ has been dropped)whereA similar set of equations can be easily derived for the variable $y$. To complete the sampling procedure, a method is required to draw samples from conditional distribution of ${\rho }_{xy}$. Because this distribution does not have an analytical form, the samples were drawn by using the Griddy Gibbs sampler (Barnard et al. 2000; Ritter and Tanner 1992). The sampling strategy used in BaNCorr was adopted from Chevrier (2007).

The goal of PCA is to linearly map ${\mathit{x}}_n$ from a $p$ dimensional real to a $q$ dimensional real space so that the transformed data ${\mathit{z}}_n$ (${\mathit{z}}_n\in {\Re }^q$; $q<p$), for a given $q$, preserves maximum variance in ${\mathit{x}}_n$. In probabilistic formulation of PCA, ${\mathit{x}}_n$ is defined as linear transformation of ${\mathit{z}}_n$ by using the following linear Gaussian latent variable model:where ${\mathit{z}}_n$ = latent variables; ${\mathit{\epsilon }}_n$ = isotropic (symmetric) Gaussian noise ${\mathit{\epsilon }}_n\sim \mathcal{N}({\mathit{\epsilon }}_n|0,{\sigma }^2{\mathit{I}}_p)$; $\mathit{\mu }$ = $p$-dimensional vector; and $\mathit{W}$ = $p\times q$ matrix. For notational simplicity, variance of the error term (${\sigma }^2$) will be represented by $\tau$, where $\tau ={\sigma }^{-2}$. Tipping and Bishop (1999) proved that if the prior distribution over ${\mathit{z}}_n$ is zero-mean Gaussian distribution, $p({\mathit{z}}_n)\sim \mathcal{N}({\mathit{z}}_n|0,{\mathit{I}}_q)$, then the maximum likelihood solution of $\mathit{W}$ contains principal directions of data matrix $\mathit{X}$ in its columns (although with a rotational ambiguity that can be easily removed as described in Golub and van Loan (1996) and Ahn and Oh (2003).

Tripathi and Govindaraju (2008), developing on the work of Sanguinetti et al. (2005), extended probabilistic PCA to engage measurement uncertainties of observed data in the analysis and developed a Bayesian framework for estimating model parameter $\mathit{W}$. The algorithm was termed as Bayesian Noisy PCA (BaNPCA) and is shown schematically in Fig. 2(a). BaNPCA was demonstrated to have the following advantages over standard PCA: (1) possesses better data compressibility, (2) yields smoother principal vectors, (3) imputes missing values in the data, and (4) identifies appropriate number of PCs. However, the algorithm was not fully Bayesian because apart from $\mathit{W}$, the other model parameters and hyperparameters were estimated using maximum a posteriori (MAP) method that approximates a posterior distribution by a delta function at its mode. This approximation is justified when $N>>p$, a case seldom satisfied by hydro-climatic data sets including sea surface temperature (SST) data. In this study, a complete Bayesian treatment of the BaNPCA algorithm is proposed; the new algorithm is referred to as variational BaNPCA (VBaNPCA). Mathematically, the algorithm is given by Eqs. (7) - (17). An important difference between BaNPCA and VBaNPCA is the introduction of true data ${\mathit{x}}_n$ as a latent variable in VBaNPCA, which substantially reduces the mathematical complexity in estimating noise variance ${\sigma }^2$, a computational bottleneck in the earlier BaNPCA algorithm.

In VBaNPCA, the following hierarchical priors are assigned to model parameters and hyperparameterswhere the column vector ${\mathit{w}}_j$ denotes $j$th column of $\mathit{W}$. The priors are made broad or noninformative by setting $a_0=b_0=c_0=d_0={\beta }_0={10}^{-3}$. The joint distribution of the observed variables, latent variables, and model parameters for VBaNPCA model [schematically shown in Fig. 2(b)] is given byThe posterior distribution obtained by using the mean-field approximation is given in the appendix.

Let the regression relationship between independent variable ${\mathit{x}}_n$ and dependent variable $y_n$ be given bywhere parameter $w_0$ that allows for a fixed offset in the data is termed as “bias”; $\mathit{w}=[w_1,\dots ,w_M{]}^T$ = weight vector; $\epsilon$ = noise term assumed to have zero-mean Gaussian distribution with variance ${\sigma }_{\epsilon }^2$ (precision ${\tau }_{\epsilon }={\sigma }_{\epsilon }^{-2}$), i.e., $\epsilon \sim \mathcal{N}(\epsilon |0,{\sigma }_{\epsilon }^2)$; and ${\mathit{\varphi }}_n=[\varphi ({\mathit{x}}_1,{\mathit{x}}_n),\dots ,\varphi ({\mathit{x}}_M,{\mathit{x}}_n){]}^T$ = vector of basis functions. In RVM, basis functions are kernel functions with one kernel associated with each data point in the training data, i.e., $M=N$. A popular choice for kernel function is Gaussian or radial basis function (RBF) kernel given bywhere ${\sigma }_{\mathrm{kernel}}$ is known as width of the kernel function. In standard RVM [Tipping 2001; Fig. 3(a)] typically used in the hydrologic literature (Khalil et al. 2005), the complexity of the regression model [Eq. (23)] is governed by assigning the following hierarchical priors, also known as automatic relevance determination (ARD) priors, over bias and weight vectors:where ${\alpha }_j$ = hyperparameter with a noninformative uniform prior in log space as with ${\sigma }_{\epsilon }^2$. Ideally, model parameters should be estimated by determining their posterior distribution by using Bayes rule [Eq. (1)]. Because this is analytically intractable, posterior distribution is estimated only for bias and weight vectors whereas point estimates are obtained for other parameters by using the Type-II maximum likelihood approach that approximates posterior distributions by a delta function at their modes. This assumption is justified when the number of data points ($N$) is very large because in that case posterior distribution is tightly concentrated around its mean. However, in typical hydrologic applications, $N$ tends to be small, and this approximation is rarely justified.

A complete Bayesian treatment of RVM through variational inference was proposed by Bishop and Tipping (2000). In this study, their method is extended to incorporate measurement uncertainties of input data into RVM learning. The methodology developed is motivated by the works of Weigend et al. (1996) for state-space models and Quiñonero Candela (2004) for Gaussian processes.

As in BaNCorr and VBaNPCA, the true input ${\mathit{x}}_n$ is treated as a latent variable in VNRVM. The schematic of the VNRVM is shown by using a DAG in Fig. 3(b). The following prior distributions were assigned to model hyperparameters and noise variance:The priors were made flat or noninformative by setting $a_0=b_0=c_0=d_0={10}^{-3}$. The joint distribution of the model parameters and observed data is given bywhere $\mathit{\alpha }=[{\alpha }_1,\dots ,{\alpha }_M{]}^T$. The mean-field approximation to the posterior distribution of model parameters and latent variables is given in the appendix.

Gridded sea surface temperature data available at monthly time step from HadSST2 (Rayner et al. 2006) prepared at the Hadley Centre for Climate Prediction and Research, U. K., 〈http://hadobs.metoffice.com/hadsst2/〉, were used in the study. HadSST2 is an un-interpolated sea surface temperature data set based on in situ measurements of SST obtained primarily from ships and buoys. The climatological values were subtracted from the SST measurements to estimate anomalies. The resulting anomalies were gridded onto a $5°$ latitude by $5°$ longitude monthly grid by taking the winsorized mean. Bias correction was applied on the gridded anomalies to remove spurious trends caused by changes in SST measuring practices, following the procedure of Folland and Parker (1995).

HadSST2 accounts for uncertainty arising from (1) measurement and sampling errors (M&SE). Measurement errors are related to precision of the instruments used for measuring SST, whereas sampling errors are caused by estimating an area-averaged quantity from a finite number of noisy observations; and (2) bias attributable to different measurement practices over the observation period. The combined uncertainty was obtained by adding M&SE and errors resulting from bias correction, in quadrature, i.e., by taking the square root of the sum of the squares of errors. This step is justified under the assumption that M&SE and errors resulting from bias correction are independent.

The area weighted rainfall data for India (Parthasarathy et al. 1994) was extracted from the Indian Institute of Tropical Meteorology, Pune, website, 〈http://www.tropmet.res.in〉. The data extend from January 1871 to December 2006. The primary source of data is the India Meteorological Department. The time series of summer monsoon rainfall was derived by adding monthly rainfall values for monsoon months (JJAS). Because the uncertainties in rainfall are not reported, they were assumed to be noise free in the study.

Synthetic data were used to evaluate the performance of the proposed algorithms because evaluation requires true observations and results that are, in general, not available with real-world data. The synthetic data sets were either obtained from the literature or designed to highlight the advantages of the proposed algorithms. They will be described in the sections in which they are used.

One hundred samples ($N=100$) from two normally distributed random variables $x$ and $y$, with mean ${\mu }_x=0$ and ${\mu }_y=0$, standard deviation $s_x=1$ and $s_y=1$, and correlation ${\rho }_{xy}=0.6$ were generated and subjected to additive zero-mean Gaussian noise with variance sampled uniformly from the interval $[0,u]$ to yield $\tilde{x}$ and $\tilde{y}$. The data were generated for different values of $u$ and correlation coefficients were estimated; the process was repeated 100 times. Fig. 4 shows the correlation estimated by the standard method and median of the 5,000 samples generated from the posterior distribution of correlation by using the BaNCorr algorithm (sampled by using MCMC algorithm after burn-in-time of 5,000 steps). As expected, the correlations estimated from the standard method deviate from the true value, whereas those from BaNCorr algorithm are closer to true values, even in presence of large noise, thereby highlighting the advantages of engaging uncertainty information.

The BaNCorr algorithm is next compared with the algorithm proposed by Kulathinal et al. (2002) in the medical science literature, which also engages data uncertainty but yields maximum likelihood estimate of correlation, abbreviated in this paper as MLECorr. Fig. 5 shows the results from the two algorithms for the data set described in the previous paragraph but for different sample sizes and $u=0.5$. For large $N$, the estimates from BaNCorr and MLECorr are similar; however, for small $N$, MLECorr tends to give higher correlation with large variability whereas BaNCorr estimates are more robust. Similarly for a given sample size, MLECorr gives higher correlation for data having large measurement errors, as shown in Fig. 6. The estimate of covariance matrix tends to be unstable in MLECorr for small sample size and large noise but is stable in the BaNCorr algorithm owing to the shrinkage effect implicit in Bayesian methods (Yang and Berger 1994; Barnard et al. 2000).

The confidence interval (CI) for the estimated correlation given by the standard method depends only on the sample size. Intuitively, measurement errors should widen CI. Nevertheless, it is not reflected in the standard method owing to its inability to engage uncertainty information. For BaNCorr results, this effect is exhibited in the form of increased dispersion in the posterior distribution of the correlation. Figs. 7 and 8 present the posterior distribution of correlation for different sample sizes and data subjected to different noise levels. The posterior distribution is distinctively different from Gaussian distribution.

The standard and BaNCorr algorithms differ markedly in their procedures for selecting a feature that reflects their frequentist and Bayesian origins. In the standard method, a feature is selected if its correlation with the dependent variable is significant, more precisely, if the null hypothesis of no correlation cannot be rejected at a given significance level. In BaNCorr, a feature is selected if the posterior probability of correlation exceeds a threshold. Both procedures have relative merits, although the Bayesian method is more straightforward to understand and implement. Additionally, the Bayesian method allows incorporation of knowledge of the physical process into the inference procedure, whenever available. For instance, a priori knowledge about the positive correlation between dependent and independent variables could be used by appropriately adjusting the prior distribution.

Tripathi and Govindaraju (2008) presented and discussed the advantages of BaNPCA over the traditional PCA methods. The desired properties of BaNPCA are also shared by the variational Bayes noisy PCA (VBaNPCA) algorithm. Hence, for brevity, only two properties, namely imputation of missing values and identification of appropriate number of PCs, will be briefly presented.

Beckers and Rixen (2003) designed a data set for evaluating the performance of data imputing algorithms; a modified version was used by Tripathi and Govindaraju (2008). The data consist of a deterministic and a random component. The $ij$th element of the deterministic component $\mathit{X}$ is given aswhere $t_i=2\pi i/n$ and $l_j=2\pi j/p$, with time index $i=1,\dots ,n$, and space or location index $j=1,\dots ,p$. For the simulations, $n$ and $p$ were fixed to 100 and 50, respectively. The random component consists of two parts; the first part mimics zero-mean normally distributed observation errors whose variance varies log-linearly between 0.05 and 5 across space and fluctuates randomly across time. The second part is an additive zero-mean Gaussian noise with standard deviation of 0.05, and it remains constant for all data points. A hundred replicates of the data set were generated and 10% of the values were randomly removed from each set. The performance is quantified by comparing imputed values with true values. Fig. 9 shows the mean absolute error (MAE) in filling deleted values by VBaNPCA, BaNPCA, and two state-of-the-art data imputation algorithms, namely Regularized EM (RegEM; Schneider 2001), and data interpolating empirical orthogonal functions (DINEOF; Beckers and Rixen 2003). BaNPCA and VBaNPCA have smaller MAE compared to RegEM and DINEOF. The median MAE from BaNPCA and VBaNPCA are similar; however, VBaNPCA estimates have less scatter, suggesting improved performance.

To test the capability of VBaNPCA in identifying the optimal number of principal components, a synthetic data set with $N=100$ samples from a multivariate Gaussian distribution in $p=10$ dimensional space, with a variance of 5 in the first four dimensions and a variance of 1 in the remaining six dimensions, was generated and subjected to (1) randomly fluctuating additive zero-mean Gaussian noise with variance sampled from the range [0, 0.5]; and (2) fixed additive zero-mean Gaussian noise with variance 0.1. The principal vectors were obtained by using standard PCA (SPCA), noisy PCA (Sanguinetti et al. 2005), BaNPCA, and VBaNPCA, and are displayed by using Hinton diagrams in Fig. 10, wherein each column represents a principal vector (the Hinton diagram provides a qualitative display of the values in a vector matrix. Each value is represented by a square whose size is associated with its magnitude, and whose color indicates its sign). Except for standard PCA (SPCA), the three other algorithms use variance of the randomly fluctuating noise, mimicking measurement errors, for estimating principal vectors. The results indicate that both BaNPCA and VBaNPCA can identify the correct number of PCs suggesting the efficacy of automatic relevance determination (ARD) priors [Eq. (18)] used in them.

Two experiments were performed to assess the performance of variational noisy RVM (VNRVM) for linear regression. The results were compared with simulation extrapolation (SIMEX) algorithm (Carroll et al. 2006), recently introduced in the hydrologic literature by Chowdhury and Sharma (2007). Briefly, in the SIMEX algorithm, a relationship is developed between noise and regression coefficient by adding known white noise to the data, which is subsequently extrapolated, typically by using polynomial extrapolation, to estimate the regression coefficient for no noise case.

In the first experiment, 100 samples of independent variable $x$ were generated from $x\sim U(x|-10,10)$ and were subjected to zero-mean Gaussian error with standard deviation drawn randomly from the range [0,5] to yield $\tilde{x}$. The dependent variable $y$ was estimated as $y=wx+\epsilon$, where $\epsilon \sim \mathcal{N}(\epsilon |0,1)$ and $w$ was set to 10. The top panel in Fig. 11 shows the value of $w$ estimated by different algorithms from 100 repetitions of the experiment, by using box-plots. Both SIMEX and VNRVM performed better in estimating the regression coefficient than standard least square regression, which ignores uncertainty in data. The results from the SIMEX algorithm were sensitive to the choice of extrapolation method, which is reflected in the results obtained from linear, quadratic, and cubic extrapolation. The bias in estimating $w$ over 100 simulations appears to be similar for ${\mathrm{SIMEX}}_C$ and VNRVM; however, VNRVM has a notably smaller variance which can again be attributed to ARD prior that regularizes estimates of $w$.

The second experiment was similar to the first, except that a redundant variable, drawn randomly from the range $[-10,10]$, was introduced into the regression problem. The experiment mimics typical problems in hydrology, in which, in the absence of a priori knowledge, redundant variables are included in regression analysis. The middle and the bottom panels in Fig. 11 show the estimated coefficients by different methods for 100 simulations of the experiment. The performance of all the methods in estimating the coefficient for the relevant variable degrades from the first experiment as reflected in slightly higher bias and increased variance. For the redundant variable, unlike SIMEX, VNRVM forces its coefficient to zero, thus effectively removing the variable from the analysis. The experiment highlights the feature-selection capabilities of linear VNRVM.

To test VNRVM for nonlinear regression, the data from the sinc function $\mathrm{sinc}(x)=\mathrm{sinc}(x)/x$ were used. The sinc function is a benchmark data set for evaluating kernel methods (Tipping 2001). The input variables were randomly generated from the range $x\in [-10,10]$ and subjected to an additive zero-mean Gaussian noise with standard deviation sampled from the range $[\begin{array}{c}0,8]\end{array}$. To the output variable $y=\mathrm{sinc}(x)$, a zero-mean Gaussian noise of standard deviation 0.1, was added. Fig. 12 shows the sinc function and the noisy data set sampled from it. RVM and VNRVM were used to learn the function from the noisy data by using an RBF kernel [Eq. (24)] with ${\sigma }_{\mathrm{kernel}}=2$. The learned functions are shown as dotted lines in Fig. 12. It is evident that the function learned by RVM is overly smooth, whereas that by VNRVM is closer to the true function.

The value for ${\sigma }_{\mathrm{kernel}}$ was assumed in the results reported for sinc data. In practice, however, it is not known and should be estimated as part of the learning process. Tripathi and Govindaraju (2007) demonstrated the importance of ${\sigma }_{\mathrm{kernel}}$ in RVM and provided a set of procedures for selecting ${\sigma }_{\mathrm{kernel}}$ for hydrologic applications. However, their methods are not directly applicable for the case with noisy input data. The use of the variational principle in VNRVM provides a handy way for selecting ${\sigma }_{\mathrm{kernel}}$. The model that maximizes the variational lower bound, $\mathcal{L}(Q)$, corresponds to the best model. Fig. 13 shows the plot of variational lower bound versus ${\sigma }_{\mathrm{kernel}}$ for the sinc data. The plot suggests ${\sigma }_{\mathrm{kernel}}\sim 2$ to be optimal, and the results reported in Fig. 12 support this choice. The selected value of ${\sigma }_{\mathrm{kernel}}$ is very close to the optimal value suggested by Tripathi and Govindaraju (2007).

Seasonal prediction of all Indian summer monsoon rainfall (AISMR) has been a topic of scientific studies for more than a century, and it is still pursued because of associated socioeconomic benefits. The literature abounds with the theories and models, both physical and empirical, but the goal of reliable prediction remains elusive (Shukla 2007; Gadgil et al. 2005).

In the early 1980s, Charney and Shukla (1981) proposed that seasonal monsoon rainfall is governed by the slowly varying surface boundary conditions of SST, soil wetness, vegetation, and snow cover. This theory, which is among the most influential theories on predictability of AISMR, forms the basis of the model developed in this study. The predictors for the model were selected from global SST data. Tripathi and Govindaraju (2008) studied the effect of engaging uncertainties in SST data on prediction of AISMR, but their focus was on predicting extreme states of AISMR by using robust optimization. In this study, predictability of seasonal mean monsoon rainfall is investigated by using the Bayesian tools developed in this study.

HadSST2 data from January 1901 to December 2006 were used in developing AISMR prediction models. The grid points with more than 50% of missing records were not considered in the analysis. The data were partitioned into five nonoverlapping sectors, (1) tropical Pacific sector (30°S–30°N), (2) North Pacific sector (north of 30°N), (3) tropical Atlantic sector (30°S–30°N), (4) North Atlantic sector (north of 30°N), and (5) Indian Ocean sector (north of 30°S), to reduce the influence of strong ENSO signals on variability of SST data outside of the tropical Pacific Ocean. The partioning of data allows for investigation of the influence of each oceanic sector on AISMR separately, and it also brings down the computational cost of processing global SST data. The monthly SST anomalies were averaged to get seasonal anomalies for four seasons (DJF, MAM, JJA, and SON). Because the measurement errors in monthly anomalies were independent, the error variance in the seasonal anomaly was estimated by taking the average of error variance in monthly anomalies. Seasonal SST tendencies, defined as a change in SST anomalies from the previous season, were estimated and the associated error variances were calculated by adding variances in the anomalies of the current and the previous seasons. The term “SST data” refers to both SST anomalies and SST tendencies. Following Sahai et al. (2008), the predictors for mean monsoon rainfall were identified from the SST data, one season lag to three years lag from the monsoon season. The data were partitioned into a training set (1901–1985) and a test set (1986–2006), and only the former set was used for identifying predictors.

To develop a prediction model, SST data for each season in an oceanic sector were processed by using PCA; the PCs obtained were correlated with AISMR to select predictors, which were then used by a regression model to learn a functional relationship between them and AISMR. Table 1 shows correlation and relative root mean square error (RRMSE) between observed and predicted values during training and testing periods. The algorithms in the first row, namely standard PCA (SPCA), Pearson product-moment correlation (Pearson), and RVM, do not use information on SST uncertainty, which is used by the algorithms in subsequent rows, with the last row using it in all the three steps. The VBaNPCA can automatically identify optimum number of PCs, whereas no objective way exists of identifying the number of PCs in SPCA. Therefore, to facilitate comparison of different algorithms, only the first three PCs of SST data from SPCA and VBaNPCA were used.

A PC was selected as a potential predictor if its Pearson product-moment correlation with AISMR was greater or smaller than zero at 10% significance level, or if the posterior probability of the absolute value of the correlation to be greater than zero was more than 0.9, i.e., $p(|r|>0)>0.9$ in which $p(r)$ represents posterior distribution of correlation obtained from BaNCorr algorithm. The regression relationships were learned by using linear versions of RVM and VNRVM algorithms. These algorithms can identify relevant predictors while learning the regression relationship. The first step of deriving PCs reduced the dimensionality of predictors from a few thousand to a few hundred. The screening of PCs by correlation futher reduced the number of potential predictors to approximately 50. Finally, the RVM identified approximately 10 out of those 50 features to be relevant for learning the regression relationship. Most of the finally selected predictors were from Tropical Pacific and India Ocean sectors.

The results presented in Table 1 indicate that the models engaging data uncertainty information perform better than the models that do not use uncertainty information for both training and testing phases. The VBaNPCA algorithm appears to significantly affect the prediction performance, whereas other two methods, BaNCorr and VNRVM, only marginally alter the prediction performance. Fig. 14 shows the observed and predicted time series during testing phase of the two models—one that does not engage uncertainty (SPCA-Pearson-RVM; NoUnc) and another that incorporates uncertainty in all three steps of model development (BaNPCA-BaNCorr-VNRVM; WithUnc). The gray shaded region represents one standard error of prediction. The results indicate that NoUnc model is poor in capturing extreme events, a behavior that was also noted in sinc data example (Fig. 12). The WithUnc model can predict extreme years of 1987, 1988, and 1994; however, it could not successfully predict weak monsoon years of 2002 and 2004. Furthermore, normal monsoon year of 1997 is wrongly predicted to be a weak monsoon year by the model.

The model that utilizes SST measurement uncertainties results in a slightly larger standard error of prediction (Fig. 15). The standard error in NoUnc is attributable to model uncertainty and uncertainties in the estimate of model parameters. In addition to these two components, measurement uncertainties in the SST data also contribute toward the prediction uncertainty in WithUnc. Among the different components of prediction uncertainty, the model uncertainty was observed to be largest. The model uncertainty was 53.5 mm for NoUnc and 46.4 mm for WithUnc. Engaging data uncertainty information marginally reduces the model uncertainty, but it is still very large, and therefore, the effects of data uncertainty are not conspicuous in the overall estimates of prediction uncertainty.

A small sample size and relatively large pool of potential predictors make assessment of prediction skill difficult. Traditionally, data are divided into a training and a validation set; the model is trained on the former set, and the performance of the trained model on the latter set gives an estimate of the prediction skill. This estimate is not reliable if the sample size of validation set is small. For such cases, cross-validation (CV; Wilks 2006) gives better estimate of prediction skill. In a typical $k$-fold CV, the data are divided into $k$ disjoint subsets, the model is trained on all subsets except one, and the prediction is made on the subset left out. The procedure is repeated $k$ times, each time selecting a different subset for validation, such that each data point becomes a part of the validation set once. If the predictors for the model are selected from the same data used for validation, it induces selection bias, a subtle point that is sometimes ignored not only in hydrology but even in statistical and machine learning literature (Ambroise and McLachlan 2002). The problem of selection bias becomes acute when the set of potential predictors is very large compared to sample size. In such cases, a reliable estimate of prediction skill requires two-level cross-validation (Wessels et al. 2005). In this procedure, at the first level, data are divided into $G$ disjoint sets, with $G-1$ subsets forming training data and the remaining set forming validation data. Training data undergoes an inner loop operation. In the inner loop, the $k$-fold CV is used to select predictors. A model is then trained, for the selected predictors, on the training data. The performance of the trained model is then evaluated on the validation data. The process of splitting the data into a training and a validation set is repeated $L$ times. Average performance over $L$ cycles gives an estimate of the prediction skill. Table 2 presents the CV results for NoUnc and WithUnc models by using the naïve $k$-fold CV with $k=N$ (leave-one-out CV), and two-fold CV with $G$, $k$ and $L$ set to 3, 5, and 100, respectively. It is evident from the table that the presence of selection bias in the naïve $k$-fold CV gives an overly optimistic estimate of prediction skill. The two-fold CV gives a more conservative estimate of prediction skill than that obtained from the testing data set in Table 1.

Predicting extreme monsoon states is the most challenging aspect of AISMR prediction. Table 3 presents the performance of NoUnc and WithUnc models for predicting below normal, normal, and above normal states of AISMR in training, testing, and cross-validation. The years with rainfall within 5% of the long term AISMR mean are termed as normal years. The below (above) normal years are years with rainfall 5% below (above) the long term AISMR mean. This table reinforces the earlier conclusion that engaging unceratinty improves the ability to predict extreme years. Neverthless, the WithUnc model could predict only approximately one-third of the extreme years, which is disappointingly low.

The following convention is used in the directed acyclic graphs (DAGs) shown in Figs. 1–3: (1) a circle or a node denotes a random number, observed variables are given by shaded circles, whereas the nonshaded circles represent the model parameters or hidden variables, (2) solid elliptical dots represent deterministic parameters, (3) an arrow denotes a distribution of the variable on its head conditioned on the variable on its tail, for e.g., an arrow connecting $x$ to $y$ represents $p(y|x)$, and (4) a plate marked $n=1,\dots ,N$ denotes that $N$ nodes exist of the type of which only one example is given inside the plate.

The mean-field approximation of the posterior distribution is given by $Q(\mathit{X},\mathit{Z},\mathit{W},\mathit{\mu },\mathit{\alpha },\tau )=q(\mathit{X})q(\mathit{Z})q(\mathit{W})q(\mathit{\alpha })q(\mathit{\mu })q(\tau )$, and by using Eq. (5), the optimal solution for the individual components takes the following closed form:where column vector $\acute{{\mathit{w}}_i}$ denotes $i$th row of $\mathit{W}$, andThe expected values required in the equations were evaluated by using the following properties: For a Gaussian distribution, $p(\mathit{x})\sim \mathcal{N}(\mathit{x}|\mathit{\mu },{\mathit{P}}^{-1})$and for a Gamma distribution, $p(x)\sim \mathrm{Gamma}(x|a,b)$

The mean-field approximation to the posterior distribution of VNRVM parameters and latent variables is given byThe optimal solution for $q(\mathit{X})$ depends on the choice of kernel function and will be discussed shortly. Optimal solution for other components of $Q$ were analytically estimated aswhereand the required moments are estimated by using Eqs. (30) - (31). For linear VNRVM, i.e., ${\mathit{\varphi }}_n={\mathit{x}}_n$, the optimal solution for $q(\mathit{X})$ has a closed form given bywhereFor RBF kernel, $q^{\star }(\mathit{X})$ is given bywhere ${\mathbb{E}}_{\setminus \mathit{X}}$ denotes expectation with respect to all the variables in the model except $\mathit{X}$. The distribution $q^{\star }(\mathit{X})$ does not resemble any standard distribution, and hence required moments for VNRVM have to be approximated. An easy choice is to approximate it by a Gaussian distribution by using either variational or Laplace methods. However, because $q^{\star }(\mathit{X})$ is often multimodal, Gaussian approximation gives a poor estimate for higher order moments. Stochastic approximation by using the MH algorithm yielded better estimates for moments and is therefore recommended for use in VNRVM. The jumping distribution was selected to be $\mathcal{N}({\mathit{x}}_n|{\tilde{\mathit{x}}}_n,{\mathit{B}}_n^{-1})$. Thus, for nonlinear case, VNRVM combines variational and MCMC methods; the usefulness of such a combination has been demonstrated in the context of other statistical models (de Freitas et al. 2001).