A Comparison of Three Soft Computing Techniques, Bayesian Regression, Support Vector Regression, and Wavelet Regression, for Monthly Rainfall Forecast

2.1 Study Area

Assam is a northeastern state of India that is located between latitudes 24°8′N to 28°2′N and longitudes 89°42′E to 96°E. The location of the state on the India map is shown in Figure 1. The state is located at the foothills of the eastern Himalayas and receives an average annual rainfall of around 300 cm from the southwest monsoon. There are around 51 forests and subforests in the states because of its topographical and climatic features. The state experiences heavy rainfall and floods every year, which cause riverbank erosions and landslides in many parts of the state. Rainfall and floods also affect agriculture, which accounts for more than a third of the state domestic product and is the principal occupation of the rural people, who constitute nearly 90% of the total population [41].

Figure 1:

Location of Assam and Stations on India Map, and Distribution of Monthly Mean Rainfall over Assam.

The numbers represent the stations.

There is significant spatial and temporal variation in rainfall over the state due to topographic and geographic features. Southwest monsoon brings heavy rainfall during the June to August period and causes flood in many parts of the state. The state suffers huge monetary losses every year because of heavy rainfall and floods. The amount of rainfall received by some regions in monsoon months is >800 mm. During the non-monsoon period, the monthly mean rainfall remains much smaller. A box plot showing the annual mean rainfall received by 21 stations in Assam between 1901 and 2002 for every month is presented in Figure 2. Due to this high temporal variability, prediction of rainfall becomes very difficult.

Figure 2:

Box Plot Showing Annual Mean Rainfall at All Stations for Every Month Between 1901 and 2002.

2.2 Data

The monthly mean precipitation data of 21 stations spread over the whole state is used for this study. Data for precipitation were obtained from the India water portal website (http://www.indiawaterportal.org/met_data/). The monthly mean rainfall is available for 102 years from January 1901 to December 2002. The distribution of rainfall is shown in Figure 1. The eastern part of the state receives less rainfall, compared to other parts of the state. The rainfall statistics of 21 stations for 102 years data are presented in Table 1. The maximum recorded monthly rainfall is 1594.27 mm at Kamrup station, which occurred in July 1974. The mean monthly rainfall varied between 172 and 294 mm. The standard deviation varied between 155.35 mm at station no. 21 (Tinsukia) and 301.43 mm at station no. 11 (Kamrup). The skewness varied between 0.67 at Lakhimpur station and 0.92 at Hailakandi station. Skewness is positive at all the stations, which shows that the frequency distribution of rainfall is asymmetric and it lies on the right side of the mean. The distribution of the data is less scattered as the values of skewness become lower. Out of 102 years data, the first 60 years data (from January 1901 to December 1960) were used for calibrating the models and the remaining 42 years data (from January 1961 to December 2002) were used for validation.

3.1 Bayesian LR

In BR, the uncertainty in β is represented by a probability distribution, p(β). Bayes’ theorem is used to express the posterior distribution for β as the product of prior distribution and likelihood function:

The value of β is chosen, which maximizes the posterior distribution. We integrate over the distribution of β to make predictions. The integration over the distribution overcomes the problem of overfitting, as we are averaging over a large number of possible solutions. The maximization of the posterior distribution in Eq. (2) is equivalent to the minimization of the negative logarithm of the distribution. Different samples of the solution are generated using the mean and variance of the distribution. Solving the minimization problem, the posterior can be expressed in terms of mean and variance.

For this study, we have assumed a Gaussian prior on β (regression coefficients) with mean zero and covariance matrix Λ⁻¹ so that β~N(0, Λ⁻¹), where Λ is the inverse covariance or precision matrix that is used to make computation easier. The posterior distribution on β can be represented as

where

Posterior mean, μ_β=(X^TX+σ²Λ)⁻¹X^Ty,

Covariance matrix, Σ_β=σ²(X^TX+σ²Λ)⁻¹.

This model works in two steps. In the first step, the mean and covariance matrix for posterior distribution of β is calculated using Eq. (3). In the second step, β is approximated from this normal distribution by averaging over a large number of samples generated by using multivariate normal distribution:

3.2 Support Vector Regression

SVM is a method of supervised learning that is based on SLT. It was developed by Vapnik and his team, and is generally used for classification and regression analysis [33, 34]. The most important advantage of SVM is the principle of structural risk minimization, which is shown to be better than traditional empirical risk minimization [11, 24]. It tries to minimize the upper-bound generalization errors instead of minimizing the local training errors. It tries to find a decision rule that provides better generalization by selecting some subset of training data called support vectors. To take into account the non-linearity present in the data, the kernel functions are used for non-linearly mapping the input space to a higher-dimensional feature space.

In SVR, a loss function (L_ε(t, f(x))) is used, which describes the deviation between the function f(x_i) outcomes and the actual target values t_i, and allows some deviation between the target values t_i and the function values f(x_i). Vapnik’s ε-insensitive loss function is used, which has some algorithmic advantages like providing sparse decomposition and rendering the computation more amenable [35]. Figure 3 shows Vapnik’s ε-insensitive loss function. It is defined as

Figure 3:

Vapnik’s ε-Insensitive Loss Function.

The points that lie outside the ε-insensitive zone (outside red dotted lines) add to the cost. The deviations outside the ε-insensitive zone are penalized in a linear manner. The associated convex quadratic programming problem is

where ξ_i and ξ_i∗ are slack variables. These variables represent the deviations for the training data points lying outside the ε-insensitive zone. The slack variables have zero value for all the data points lying inside the ε-insensitive zone and increase progressively for the data points lying outside the zone. C is the regularization constant that assigns a penalty when a training error occurs. The current optimization problem can be easily optimized in the dual form.

A dual set of variables, α_i and α¯i, is introduced to build the Lagrange function from the primal objective function [Eq. (6)] and the corresponding constraints [Eq. (7)]. The dual formation of the optimization problem is

Subject to∑i = 1n(α_i−α¯i)=00≤α_i≤C,i=1, 2, …, n0≤α¯i≤C,i=1, 2, …, n.

After determining α_i and α¯i, the Karush-Kuhn-Tucker method is utilized to find the parameters w and b. Prediction in an LR function is expressed as

To make SVR deal with non-linear cases, a function ϕ is used to map the decision space variable (X) to a higher-dimensional space. The dual problem, after applying this transformation, becomes

where k(x_i, x_j)=(ϕ(x_i).ϕ(x_j)) is called the kernel function. Any function can be used as kernel function if it satisfies Mercer’s theorem [34]. Gaussian radial basis function (GRBF) is used for this study, which can be defined as K(x_i, x_j)=exp(−γ||x_i−x_j||²).

The hyperparameters such as the SVR constant ε, regularization constant C, and σ (for radial basis function kernels) influence the performance of SVR. Parameter C allows a trade-off between training error and model complexity. If a model is not sufficiently complex, then it may fail to capture the underlying trend of the data and hence cause underfitting. On the other hand, if the model is too complex, then it may capture the noise of the data and, hence, can suffer from overfitting.

The magnitude of training error increases for the small values of C and hence tends to underfit the training data. However, a higher value of C will tend to a behavior as that of a hard-margin SVM and hence overfits the training data [13]. The parameter ε influences the number of support vectors by deciding the thickness of the ε-insensitive zone. Hence, its value affects both generalization capability and the complexity of the approximation function. A lower value of ε leads to higher number of support vectors and increases the complexity. A higher value of ε leads to smaller number of support vectors and results in more flat estimates of the regression function. The performance of SVR is sensitive to these parameters; hence, it is important to find appropriate values for C and ε [14]. Determining the appropriate value of these parameters is often a heuristic trial-and-error process [25]. It is shown that an exponentially growing sequence of parameters works better [29].

3.3 Wavelet Regression

Wavelet analysis is an advanced signal processing technique that has been extensively used in communications, data compression, acoustics, signal processing, image processing, etc. It was introduced by Grossmann and Morlet [10] in the early 1980s. The signal or function is expressed as the combination of “little waves” called wavelets. Unlike Fourier transform, where the signal is divided into smooth sinusoids of unlimited duration, the wavelet transform splits the signal into wavelets of finite duration and zero mean [16]. Wavelets are localized both in frequency and time domains. This property of wavelets is called time-frequency localization. Wavelet analysis features the analysis of temporal variations at different time scales, and can work on non-stationary signals as well. Fourier transform assumes infinite length signals and concentrates mainly on their frequencies only while wavelet transform can be applied to any kind and any size of the signal.

The wavelet function ψ(t) is defined as ∫−∞+∞ψ(t)dt=0. The term ψ_a,b(t) is obtained by compression and expansion of ψ(t):

where R represents real numbers, ψ_a,b(t) is the successive wavelet, and a and b are the scale of frequency factor and time factor. When the term ψ_a,b(t) satisfies Eq. (12), for a finite energy signal [16, 26, 40], the successive wavelet transform of f(t) is

where ψ(t)¯ represents the complex conjugate function of ψ(t). f(t) is decomposed at a different resolution level (scale) in Eq. (13). In real applications, the successive wavelets are generally discrete. Let a=a0j, b=kb0a0j,a₀>1, b₀∈R, and k, j are integers. The DWT of function f(t) is defined as

The widely used values for parameter a₀ is 2 and for b₀ is 1 time step. When a₀ is 2 and b₀ is 1, Eq. (14) changes into binary wavelet transform [16] as

W_ψf(a, b) or W_ψf(j, k) reflects the frequency characteristics (a or j) and time characteristics (b or k) of the original time series. The parameter a decreases when the temporal resolution is high and the frequency resolution of the wavelet transform is low. On the other hand, a increases when the temporal resolution is low and the frequency resolution is high [36]. For discrete time series f(t) for different time steps (integral), the DWT is defined as [42]

where W_ψf(j, k) is the wavelet coefficient for the discrete wavelet of scale a=2^j and b=2^jk.

The input signal (monthly rainfall time series) is decomposed into a certain number of sub-time series components (Ds) using the Mallat DWT algorithm [21]. Every component has an important role in the original time series, and all have distinct behavior. Thus, the contribution of different subcomponents to the original time series varies from each other [36]. The decomposition process consists of a number of successive filtering steps. The original signal is first decomposed into approximation and detail components. It is broken down into many lower-resolution components. The details are low-scale high-frequency components of the signal, while the approximations are the high-scale, low-frequency components (Figure 4). The higher scales consist of the extended version of a wavelet, and the coefficients refer to the slowly-changing coarse features of low-frequency components. The lower scales present the condensed wavelet and follow the rapidly changing details (high-frequency components) of the signal [21].

Figure 4:

Signal Decomposition by DWT for Three Resolution Levels.

S_t denotes the original signal, A_i denotes the approximation components, and D_i denotes the detail components.

The WR model used in the present study is a combination of DWT and LR. DWT is utilized for preprocessing of the original monthly rainfall series before use in the LR model. The decomposed subcomponents (detail and approximation) are utilized as inputs for LR, keeping the original rainfall series as input. The rainfall time series is processed using DWT separately for both calibration and validation period.

3.4 Models

Using these three techniques and two different sets of input variables, six different models were developed. The details of these models are provided in Table 2. All the models were calibrated and validated for the same data, as specified in Section 2.2.

4.1 Comparison of BR, SVR, and WR

Four different statistics, namely mean absolute error (MAE), root mean square error (RMSE), correlation coefficient (R), and Nash-Sutcliffe efficiency coefficient (NS) [43], were used to evaluate the performance of the three techniques. These statistics are defined as follows:

1. The MAE between the observed and the predicted outputs can be defined as

   MAE=1n∑i = 1n|fi−yi|.

2. The RMSE between the observed and the predicted outcomes can be defined as

   RMSE=1n∑i = 1n(fi−yi)2.

3. R is defined as

   R=∑i = 1n(yi−y¯i)(fi−f¯i)[∑i = 1n(yi−y¯i)2·∑i = 1n(fi−f¯i)2]12.

4. NS is defined as

   NS=1−∑i = 1n(fi−yi)2∑i = 1n(yi−y¯i)2,

where f_i are the estimated values, f̅_i is the average of estimated values, y_i are the observed values, and y̅_i is the average of observed values.

Figure 6A–D shows the comparison of the performance of six models during the validation period. As observed from the figure, the WR2 model outperformed all other models for all 21 stations. The MAE for WR2 was observed to lie between 62.89 mm at station no. 19 and 119.813 mm at station no. 11. The RMSE for WR2 was found to lie between 93.06 mm (for station no. 15) and 185.79 mm (for station no. 11), while the RMSE for SVR2 model was observed between 106.15 and 214.08 mm. The average RMSEs for BR1, SVR1, WR1, BR2, SVR2, and WR2 were 172.84, 151.24, 142.58, 173.89, 145.42, and 127.06 mm, respectively. The WR2 model had the least RMSE compared to all other models.

Figure 6:

(A) MAE for Different Models at 21 Stations. (B) RMSE for Different Models at 21 Stations. (C) Correlation Coefficient (R) for Different Models at 21 Stations. (D) Nash-Sutcliffe Index (NS) for Different Models at 21 Stations.

Figure 6C shows the R for different models at all stations during the validation period. The value of R for WR2 was found to be between 0.794 at station no. 21 and 0.855 at station no. 8. The average values of R for BR1, SVR1, WR1, BR2, SVR2, and WR2 were 0.667, 0.702, 0.750, 0.643, 0.729, and 0.828, respectively. Figure 6D shows the Nash-Sutcliffe efficiency index (NSE) for different models at all stations during the validation period. The maximum and minimum values of NSE for WR2 were 0.694 at station 8 and 0.578 at station 21. The average values of NSE for BR1, SVR1, WR1, BR2, SVR2, and WR2 were 0.334, 0.491, 0.547, 0.327, 0.529, and 0.640, respectively. Considering R and NSE, it is clear that the WR2 model performed better than all other models.

Rainfall is one of the most complex and challenging hydrological processes to understand and to model because of the complexities involved with the associated atmospheric processes that generate rainfall. The rainfall time series is affected by many complex factors having different frequency components. Considering only one resolution component, the model is unable to capture the internal mechanism behind the rainfall process. Using DWT, the complex rainfall time series were decomposed into several simple time series, and the WR model was developed utilizing appropriate subcomponents. The performance of WR model was better than the other models due to the fact that it uses different resolution components for LR. The ineffective frequency components were removed from the rainfall time series, which effectively improved the performance of the LR model. WR is a simpler model compared to the non-linear SVR model; however, WR has shown better accuracy in rainfall prediction. This signifies the role of wavelet preprocessing of the rainfall time series, which de-noised the inputs and resulted in improved prediction ability of the model.

To assess the spatial variations in the performance of the forecasting models, the state was divided into four zones: lower Assam, middle Assam, upper Assam, and southern Assam (Figure 7). Table 4 presents the average RMSEs for four regions, obtained using different models. A comparison of model performances reveals that the upper Assam part has the highest accuracy (least RMSEs), whereas middle Assam has the least accuracy (highest RMSEs). The average standard deviations of monthly rainfall for upper, middle, lower, and southern Assam were 168, 254, 222, and 259 mm, respectively. The upper Assam region had the minimum average standard deviation, which could be the reason for the better performance of models in this region, as variability in station rainfalls for this region was least. Figure 7 shows the performance of the WR2 model (best model) in the form of a map made using inverse distance weighted interpolation of RMSEs at different stations. The average RMSEs for lower, middle, upper and southern Assam were 126.26, 154.33, 95.31, and 141.25 mm, respectively. Upper Assam had the lowest RMSE followed by lower, southern, and middle Assam, respectively. The WR2 model performed best for the stations in the upper Assam zone. Stations in the middle Assam had the highest RMSEs, which imply that the forecasting model performed badly at these stations. Of 21 stations, station no. 15 (Lakhimpur) had the least RMSE for the WR2 model, whereas station no. 11 (Kamrup) had the highest RMSE. It is already shown in Table 1 that the Kamrup station had the highest values of standard deviation and skewness.

Figure 7:

Map Showing the Distribution of RMSEs for the WR2 Model over Assam.