Bayesian Surprise-Enhanced Random Forests for Rare Anomaly Classification in Environmental Sensor Data

1. Introduction

Advances in information technology have enabled the deployment of sensors to monitor environmental phenomena, such as air and reservoir water quality. For example, the Taiwanese government has installed environmental sensors to measure total phosphorus, chlorophyll a, and transparency to warn of critical river water quality. The government shares historical sensor measurements with residents through websites. Flagging anomalous sensor readings is governed by specific rules published on these websites.

Nevertheless, recent severe climatic changes or extreme events may cause environmental sensor data to undergo distribution shifts. New observations and old sensor measurements follow different probability distributions. Thus, it is difficult to apply a single probability distribution to model the sensor data and flag anomalous sensor measurements. It may therefore be necessary to recalibrate frequently existing rules to accommodate new anomalies.

Machine learning algorithms can help mitigate the burden of frequent recalibration of existing rules. In practice, training tree-based models for accurately processing large volumes of data has become feasible and computationally efficient. Such models may be trained to separate normal and anomalous sensor data, and human experts next inspect the resulting predictions to decide whether to accept or override the detected anomalies. However, a challenge may be encountered in implementing such an idea. Since anomalous data are naturally infrequent, tree-based models often fail to achieve balanced classification performance between normal and anomalous sensor data [1,2]. They may output high weighted classification metrics, but anomaly-class recall values remain low, implying that excessive false negatives have been generated.

This study proposes that integrating Bayesian surprise (e.g., [3]), defined using relative predictive surprise ($RPS$), with a random forest classifier may balance weighted classification metrics and the class-specific classification accuracy. $RPS$ values are first computed, standardized, and then added as anomaly-informed features for real air quality, reservoir water quality, and river water quality datasets. Before augmenting with these $RPS$ values as new features, the resulting predictive surprise values are compared to the Rousseeuw–Croux estimators [4], the Tukey interquartile range [5], and a simple percentile-based rule to ensure that they correctly flag anomalies on complementary cumulative probability curves. A random forest classifier is next employed to classify the augmented datasets. Based on empirical results, this study investigates whether the corresponding weighted classification metrics and the class-specific classification accuracy are simultaneously improved after adding the $RPS$ values as new features.

This study anticipates that this research will fill three research gaps to make the corresponding contributions:

• The first research gap is the lack of a computationally simple anomaly measure for environmental sensor data, which usually exhibit heavy tails. Calculating the $RPS$ requires fewer steps than the computation of Rousseeuw–Croux estimators.
• The second research gap concerns the performance of tree-based models in classifying environmental sensor data. Since anomalous environmental sensor measurements are naturally infrequent, tree-based models tend to focus on improving weighted classification metrics and may cease to optimize the classification accuracy of anomalous classes once the weighted accuracy is maximized [1,2]. This study aims to balance weighted classification metrics and the class-specific classification accuracy by aggregating $RPS$ values as new anomaly-informed features.
• The third research gap concerns the robustness of anomaly thresholds for heavy-tailed sensor data. When existing anomaly measures (for example, Rousseeuw–Croux estimators and percentile rules) are applied to different heavy-tailed distributions, the thresholds for separating normal and abnormal sensor readings may need to be re-tuned for more accurate separation. In contrast, relative predictive surprise is robust to different tail behaviors. This study adopts a fixed anomaly-detection criterion based on $RPS$ across all heavy-tailed environmental sensor datasets.
• The final contribution of this study is enhanced explainability. Without an anomaly measure, a tree-based model is like a ’black box’ in classifying normal and anomalous sensor data. It provides only limited explainability for its classification results. This study integrates a tree-based model with $RPS$ values as anomaly-informed features to provide distribution-level information about how unexpected each sensor reading is, thereby improving the explainability of the resulting classifier.

The remainder of this article is organized into five sections. Section 2 presents a literature review relevant to this study. Section 3 elucidates Bayesian surprise and provides a brief review of the random forest classifier. Section 4 reports experimental results that test the combination of relative predictive surprise values and the random forest classifier on three real datasets of air quality, reservoir water quality, and river water quality. Based on these results, Section 5 offers a discussion. Section 6 presents the conclusions of this study.

2. Literature Review

Because environmental sensor data commonly exhibit heavy-tailed behaviors, the first part of this literature review focuses on prior studies that proposed anomaly measures for heavy-tailed data. The second part reviews published studies on data augmentation methods for heavy-tailed data, since this study calculates $RPS$ values from real environmental sensor datasets and uses them as additional anomaly-informed features.

2.1. Anomaly Measure for Heavy-Tailed Data

A brief review of the literature on anomaly measures for heavy-tailed data identifies two important research gaps. First, a reliable threshold for distinguishing normal observations from anomalous ones has not yet been established. If a reliable threshold is unavailable, augmenting with anomaly-informed features provided based this threshold may not be confidential. Second, the computation of existing anomaly measures often remains unnecessarily complex. In published studies, anomaly measures for heavy-tailed data have generally been grouped into three main categories:

The first approach is based on extreme value theory and tail modeling [7,8,9,10,11]. Rather than modeling the entire heavy-tailed distribution to distinguish normal observations from anomalies, this approach focuses on the tail behavior and fits the tail using a specific probability distribution, such as the generalized Pareto distribution. Anomalies are then identified on the basis of the fitted tail model. In general, the choice of tail model depends on the heaviness and shape of the tail. Representative studies suggest that the Hill estimator [7] is one of the most widely used methods for characterizing tail heaviness. To address finite-sample issues, a previous study [8] proposed a trimmed version of the Hill estimator. Building on this idea, Bhattacharya et al. [9] developed an automatic, data-driven procedure for flagging observations that deviate from the central and intermediate structure of the data. Verardi and Croux [10] extended Tukey’s interquartile range method [5] by incorporating the medcouple to account for skewness in both univariate and multivariate settings, thereby producing an adjusted boxplot for outlier detection. Li et al. [11] estimated sample skewness and constructed empirical cumulative distribution functions to identify tail regions, which were then used to aggregate tail probabilities and define an outlier score.

The second approach comprises probability- and kernel-based methods for modeling heavy-tailed distributions. In this approach, researchers first derive parametric or nonparametric expressions to describe a heavy-tailed distribution, and then identify anomalies based on the resulting probability models. Horváth et al. [12] employed copula functions to estimate the joint distribution of a multivariate heavy-tailed random vector and constructed an anomaly score from the resulting copula-based estimation. A copula enables separate modeling of the marginal distributions and their dependence structure, capturing only the dependency among variables in a multivariate distribution. Bourigault and Mandic [13] combined generalized hyperbolic processes with kernel density estimation and one-class support vector machines to develop anomaly detection frameworks, where the anomaly score is defined using the fitted probability density function. Tortora et al. [14] proposed a model-based clustering approach based on mixtures of multiple scaled contaminated normal distributions to effectively accommodate skewness, heavy tails, or their combination in the classification of heavy-tailed data. This framework allows (1) flexible cluster shapes in terms of skewness and kurtosis and (2) component-wise and directional outlier detection.

The third approach focuses on deep learning models for representing complex, high-dimensional heavy-tailed data distributions. However, most of these models were originally designed for detecting anomalous images. A key distinguishing feature among the models lies in their learning mechanisms. For example, Zenati et al. [15] investigated generative adversarial networks and proposed adversarially learned anomaly detection, in which images with large reconstruction errors are identified as anomalous. Dionelis [16] also explored generative adversarial networks and adopted a similar reconstruction-based criterion for anomaly detection; however, their learning objective is to reconstruct the tail region of a heavy-tailed data distribution. Yang et al. [17] developed a transformer-based reconstruction module for image anomaly detection.

2.2. Data Augmentation Methods for Heavy-Tailed Data

This study proposes that augmenting heavy-tailed data with $RPS$ values may improve the accuracy of classifying these data using a random forest classifier. A concise review of prior studies on data augmentation methods for heavy-tailed data reveals two research gaps. First, existing data augmentation methods for heavy-tailed data are largely restricted to specific application domains (for example, image classification), and their implementation can be inefficient because they typically require tuning many hyperparameters. For instance, a recent study [6] developed a distance-based scheme to enhance the performance of a random forest classifier; however, training this classifier must be performed twice to tune certain hyperparameters.

Other relevant studies in image classification are also available. Ma et al. [18] observed a substantial imbalance between the numbers of normal and abnormal images in an image classification problem. To address this issue, they proposed a four-phase training scheme for a conventional neural network. The second phase of this scheme incorporates a feature-uncertainty representation based on a Gaussian distribution, which provides probabilistic features of the decision boundary separating normal and abnormal images. Introducing this probabilistic feature improves the image classification performance of a convolutional neural network.

Song et al. [19] also investigated image classification problems. They developed a contrastive visual–language model that consists of visual and semantic components. The visual components were created using images, whereas the semantic features are generated using the text sequences associated with the images. Experimental results demonstrated that the proposed contrastive visual–language model improves the F1 scores for classifying normal and abnormal images.

Delaigle and Hall [20] investigated the effect of heavy tails on cardiomyopathy microarray data and Affymetrix spike-in data. They calculated statistical features to conduct their analysis. The purpose of this calculation is to recover influential genes that would otherwise be incorrectly discarded due to the presence of outliers.

3. Methods

Suppose ${\mathbf{X}}_1,{\mathbf{X}}_2\dots ,{\mathbf{X}}_N$ represent a collection of N environmental sensor measurements. These ${\mathbf{X}}_1,{\mathbf{X}}_2\dots ,{\mathbf{X}}_N$ may be acquired asynchronously. Each ${\mathbf{X}}_i$ is a an n-dimensional tuple, $\left(x_{i1},x_{i2}\dots ,x_{in}\right)(i=1,2\dots ,N)$. Each $x_{ij}(j=1,2\dots ,n)$ may be either normal or anomalous.

3.1. Bayesian Surprise

Suppose a null hypothesis ${\mathbf{H}}_0$ states that the collected $x_{ij}(i=1,2\dots ,N;j=1-n)$ sensor measurements follow a specific probability distribution, such as the Student’s t distribution. A new observation $x_{obs}$ may be a newly observed value for the j-th sensor variable. It may be an outlier of $x_{1j},x_{2j}\dots ,x_{Nj}$ if its resulting p value is below a predetermined significance level $\alpha$ (e.g., $\alpha =0.05$). This outlier is a potential anomaly. The p-value is calculated by [3]

$$p=Pr\left(t_{ij}\ge t_{obs}\right)\le \alpha$$

(1)

where $t_{ij}={\displaystyle \frac{x_{ij}-{\mu }_j}{s_j}}$ and $t_{obs}={\displaystyle \frac{x_{obs}-{\mu }_j}{s_j}}$, ${\mu }_j$ is the mean of referenced $x_{1j},x_{2j}\dots ,x_{Nj}$ sensor measurements, s is their standard deviation and Pr denotes the probability. This Pr is equal to ${\displaystyle \frac{\#(t_{ij}=t)}{N}}(i=1,2\dots ,N)$ in which $\#(t_{ij}=t)$ denotes the number of observations equals to t. However, the distribution of $x_{1j},x_{2j}\dots ,x_{Nj}$ may exhibit heavy tails. The p value-based anomaly detection may therefore yield excessive false positives. To address this issue, scholars proposed alternative approaches (e.g., [21]) that better characterize the rarity of an observation $x_{obs}$ with respect to the reference measurements $x_{1j},x_{2j}\dots ,x_{Nj}$. In this study, considering that reference measurements $x_{1j},x_{2j}\dots ,x_{Nj}$ follow an unknown distribution, the Bayesian surprise, defined through the $RPS$ [3], is adopted to assess whether a $t_{obs}$ is rare. For discrete sensor measurements, the $RPS$ is defined by [3,21,22]:

$$RPS={\displaystyle \frac{Pr\left(t_{obs}\right)}{sup\left[Pr\left(t_{ij}\right)\right]}}={\displaystyle \frac{Pr\left(t_{obs}\right)}{\left[{Pr}_{max}\left(t_{ij}\right)\right]}}$$

(2)

where sup is the supremum, the subscript $max$ denotes the maximum, $i=1,2\dots ,N$, and $j=1,2\dots ,n$. The theoretical foundation of Equation (2) is related to the Kullback-Leibler divergence [23]. Unlike conventional approaches that seek observations inconsistent with a specific probability distribution, Equation (2) is used in this study to identify sensor readings that are relatively rare in the reference sensor data.

The surprise index can be used to determine whether $t_{obs}$ is anomalous [21,22]. We can derive:

$$SI={\displaystyle \frac{{\sum }_{i=1}^N{Pr}^2\left(t_{ij}\right)}{Pr\left(t_{obs}\right)}}\le {\displaystyle \frac{{\sum }_{i=1}^{N}Pr\left(t_{ij}\right){Pr}_{max}\left(t_{ij}\right)}{Pr\left(t_{obs}\right)}}={\displaystyle \frac{{Pr}_{max}\left(t_{ij}\right){\sum }_{i=1}^{N}Pr\left(t_{ij}\right)}{Pr\left(t_{obs}\right)}}={\displaystyle \frac{{Pr}_{max}\left(t_{ij}\right)}{Pr\left(t_{obs}\right)}}$$

(3)

where $SI$ represents the surprise index. Equation (3) indicates that the maximum of the $SI$ is equal to the reciprocal of the $RPS$. This relationship may serve as an important basis for determining the threshold at which an $RPS$ value implies an anomaly. Table 1[24] presents the referenced rules for flagging anomalies based on $SI$ and $RPS$ [24].

3.2. Random Forest Classifier

Suppose ${\mathbf{R}}_1,{\mathbf{R}}_2\dots ,{\mathbf{R}}_N$ represent the N tuples consisting of raw sensor data and their $RPS$ vslues in which ${\mathbf{R}}_i=\left(x_{i1},x_{i2}\dots ,x_{in},r_{i1},r_{i2}\dots ,r_{in}\right)$ and $r_{ij}=RPS\left(x_{ij}\right)(i=1,2\dots ,N;j=1,2\dots ,n)$ and Y is the class label. To avoid data leakage and possible confusion with the training-test split used in random forest classification, it must be emphasized that the $RPS$ value of a new observation $x_{obs}$ (or $t_{obs}$) is computed based on the N raw reference observations $x_{ij},x_{2j}\dots ,x_{Nj}$.

A random forest is a tree-based ensemble learning algorithm in which each tree is constructed using a random subset of the available information [25]. In classification tasks, it learns a prediction function $f\left(\mathbf{R}\right)$ from the training observations ${\mathbf{R}}_{i}i=1,2\dots ,N$ to predict the class label Y. The learning objective is defined through a loss function $Loss\left(Y,f\left(\mathbf{R}\right)\right)$, which quantifies the mismatch between the predicted and actual labels. For classification under the zero-one loss criterion, the optimal prediction function is obtained by minimizing the expected loss $E_{RY}\left[Loss\left(Y,f\left(\mathbf{R}\right)\right)\right]$ where the expectation $E_{RY}$ is taken over the joint distribution of $\mathbf{R}$ and Y. A commonly used loss function is [26]:

$$Loss\left[Y,f\left(\mathbf{R}\right)\right]=I\left[Y\ne f\left(\mathbf{R}\right)\right]=\left\{\begin{array}{cc}0& \mathrm{if}Y=f\left(\mathbf{R}\right)\hfill \\ 1& \mathrm{otherwise}\hfill \end{array}\right.$$

(4)

where I is the indicator function. If $\mathfrak{Y}$ is the set containing all possible Y values, minimizing the $E_{RY}\left[Loss\left(Y,f\left(r\right)\right)\right]$ for the zero-one loss results in

$$f\left(\mathbf{R}\right)=\underset{Y\in \mathfrak{Y}}{argmax}Pr\left(Y=y|\mathbf{R}\right)$$

(5)

In ensemble learning, the prediction function f is constructed from a collection of base learners $h_1(\mathbf{R},{\Theta }_1),h_2(\mathbf{R},{\Theta }_2)\dots ,h_K(\mathbf{R},{\Theta }_K)$ in which ${\Theta }_k(k=1,2\dots ,K)$ denotes the random variables used to build the k-th learner. These random variables are assumed to be independent, and the resulting base learners are aggregated to form the final prediction function f:

$$f\left(\mathbf{R}\right)=\underset{Y\in \mathfrak{Y}}{argmax}{\displaystyle \sum _{k=1}^K}I\left[Y-h_k\left(\mathbf{R},{\Theta }_k\right)\right]$$

(6)

The learners $h_1(\mathbf{R},{\Theta }_1),h_2(\mathbf{R},{\Theta }_2)\dots ,h_K(\mathbf{R},{\Theta }_K)$ build the tree classifier by recursively splitting the training data ${\mathbf{R}}_i(i=1,2\dots ,N)$ data into a splitting criterion that measures the purity of a leaf node. In this study, the $Gini$ index is chosen as the splitting criterion. It is defined as

$$Gini=1-{\displaystyle \sum _{i=1}^C}{P}_i^2$$

(7)

where C is the number of classes and $P_i$ denotes the proportion of observations belonging to class i in the node. Data partitioning based on Equation (7) is carried out recursively until a stopping criterion is satisfied. For classification of a new observation, the observation is propagated through each tree in the forest until it reaches a terminal node. Unlike a single decision tree model, random forest does not require pruning during tree growth. The absence of pruning simplifies model construction and is one of the practical advantages of random forest.

Algorithm 1 [26] presents the procedure for constructing a random forest classifier. First, bootstrap samples are generated from the observations ${\mathbf{R}}_1,{\mathbf{R}}_2\dots ,{\mathbf{R}}_N$ data. This step returns each bootstrap sample to the ${\mathbf{R}}_1,{\mathbf{R}}_2\dots ,{\mathbf{R}}_N$ by sampling with replacement. Each bootstrap sample is then used to grow a classification tree until a stopping criterion is met. During this process, the sample is recursively split to determine the best partition at each node.

Algorithm 1 Random forest classifier [26]

Input:  Training observations ${\left\{({\mathbf{R}}_i,Y_i)\right\}}_{i=1}^N$, where ${\mathbf{R}}_i=(x_{i1},x_{i2},\dots ,x_{in},r_{i1},r_{i2},\dots ,r_{in})$ contains raw sensor measurements and $RPS$-based anomaly-informed features; number of trees K; number of randomly selected features m Output:  A trained random forest classifier 1: for$k=1$ to K do 2:     Draw a bootstrap sample from ${\left\{({\mathbf{R}}_i,Y_i)\right\}}_{i=1}^N$ with replacement 3:     Initialize the root node using all observations in the bootstrap sample 4:     while the stopping criterion is not satisfied do 5:         Randomly select m features from all available features 6:         Determine the best split among the selected features according to the $G_{ini}$ index 7:         Split the current node into child nodes according to the best split 8:     end while 9:     Grow an un-pruned classification tree $h_k(\mathbf{R},{\Theta }_k)$ 10: end for 11: Aggregate the predictions of ${\left\{h_k(\mathbf{R},{\Theta }_k)\right\}}_{k=1}^K$ by majority voting 12: For a new observation, compute its $RPS$ value using the reference observations 13: return the resulting random forest classifier

In addition to the $Gini$ index, several hyperparameters must be specified when implementing a random forest classifier, including:

• Maximum depth of the tree;
• Number of trees in the forest;
• Minimum number of samples required to split an internal node; and
• Minimum number of samples required to be at a leaf node.

Section 4.4 further examines the effects of these four hyperparameters on the classification accuracy of a random forest classifier with and without the inclusion of $RPS$ values as anomaly-informed features.

4. Applications

To improve the quality of life, the Taiwanese government has deployed environmental sensors to monitor environmental conditions, including air quality, reservoir water quality, and river water quality. It has also constructed websites to share historical sensor measurements. This study uses three real datasets collected from these websites to evaluate the performance of the random forest classifier and $RPS$.

4.1. Air Quality

The first real dataset contains 24,007 records of SO₂, CO, NO₂, PM2.5, and air quality index (AQI) measurements [27]. These measurements were collected daily in 2024 from 60 monitoring stations across Taiwan. The units of SO₂ and NO₂ are ppb, whereas that of CO is ppm. The unit of PM2.5 is $\mu g/m^3$. Figure 1 depicts the distributions of SO₂, CO, NO₂, PM2.5, and AQI. Table 2 presents their summary statistics, where $IQR$ denotes the Tukey interquartile range, and ${\gamma }_{Hill}$ represents the trimmed Hill index [8].

As mentioned earlier, the trimmed Hill index [8] was used to characterize the heaviness of a tail. For the j-th sensor measurement, suppose ${\chi }_{\left(1\right)}\le {\chi }_{\left(2\right)}\le \dots \le {\chi }_{\left(N\right)}$ denote the ordered $x_{1j},x_{2j},\dots ,x_{Nj}(j=1,2\dots ,n)$ sensor data. The trimmed Hill index is defined as [7,8]

$${\gamma }_{Hill}={\displaystyle \frac{1}{k-k_0}}{\displaystyle \sum _{i=k_0+1}^k}ln{\displaystyle \frac{{\chi }_{(N-i+1)}}{{\chi }_{(N-k)}}}$$

(8)

where $1\le k\le N$, and $k_0$ denotes the number of the largest extremes trimmed from the estimation, with $0\le k_0<k$. Based on the published article [8], Table 3 presents our interpretation of the resulting ${\gamma }_{Hill}$ values. These values are used as practical reference thresholds rather than universal decision rules. In practice, $k_0$ and ${\gamma }_{Hill}$ can be computed using standard spreadsheet software.

As expected, Figure 1 and Table 2 suggest that the distributions of SO₂, CO, NO₂, PM2.5, and AQI exhibit substantial upper-tail heaviness. Evidence includes the large kurtosis values and high ${\gamma }_{Hill}$ values in Table 2, as well as the pronounced upper-tail spread shown in Figure 1. Under such distributional characteristics, applying the standard Z-score to anomaly detection may yield incorrect detections. Extreme $x_{ij}(i=1,2\dots ,N;j=1,2\dots ,n)$ can inflate their mean and standard deviation, thereby masking potential anomalies.

Furthermore, Equation (9) presents the Pearson correlations among SO₂, CO, NO₂, and PM2.5. In particular, CO and NO₂ exhibit a strong positive correlation. Equation (9) suggests that these two pollutants may share common emission sources, potentially related to fossil-fuel combustion.

$$\begin{array}{cc}& \begin{array}{cccc}\mathrm{S}{\mathrm{O}}_2& \mathrm{C}\mathrm{O}& \mathrm{N}{\mathrm{O}}_2& \mathrm{PM}2.5\end{array}\\ \begin{array}{c}\mathrm{S}{\mathrm{O}}_2\\ \mathrm{C}\mathrm{O}\\ \mathrm{N}{\mathrm{O}}_2\\ \mathrm{PM}2.5\end{array}& \left[\begin{array}{cccc}1& 0.208& 0.266& 0.284\\ 0.208& 1& 0.735& 0.592\\ 0.266& 0.735& 1& 0.434\\ 0.284& 0.592& 0.434& 1\end{array}\right]\end{array}$$

(9)

Next, this study combines a log-log plot and dumbbell plots to compare the ranges of anomalies identified by the $RPS$ and four baseline metrics. Since the ${\gamma }_{Hill}$ values reported in Table 2 suggest substantial upper-tail heaviness, the Tukey interquartile range ($IQR$), the Rousseeuw-Croux estimators $Z_{S_n}$ and $Z_{Q_n}$, and a simple percentile-based rule are adopted as baselines. According to the $IQR$, a sensor measurement $x_{ij}(i=1,2,\dots ,N;j=1,2,\dots ,n)$ is considered anomalous if it satisfies

$$x_{ij}<Q_1-1.5IQR\mathrm{or}{x}_{ij}>Q_3+1.5IQR$$

(10)

in which $Q_1$ and $Q_3$ are the first and third quartiles; respectively. Meanwhile, Rousseeuw-Croux estimators $Z_{Q_n}$ and $Z_{S_n}$ [4] are defined as

$$Z_{Q_n}=\left|{\displaystyle \frac{x_{ij}-median\left(x_{ij}\right)}{Q_n}}\right|$$

(11)

$$Z_{S_n}=\left|{\displaystyle \frac{x_{ij}-median\left(x_{ij}\right)}{S_n}}\right|$$

(12)

where, for each k, $S_n=1.1926{median}_i\left({median}_j\left|x_{ik}-x_{jk}\right|\right)$ and $Q_n=2.219{\left\{\left|x_{ik}-x_{jk}\right|:i<j\right\}}_{\left(L\right)}$. Here, $\left|·\right|$ is the absolute value, and ${\{·\}}_{\left(L\right)}$ denotes the L-th order statistic obtained by sorting all pairwise distances $\left|x_{ik}-x_{jk}\right|$ for $i<j$ in an ascending order and selecting the L-th element. Moreover, $L={\displaystyle \frac{floor\left(\frac{N}{2}+1\right)!}{2!floor\left(\frac{N}{2}-1\right)!}}$ where floor denotes the floor function, ! represents the factorial, $i,j=1,2,\dots ,N$, and $k=1,2,\dots ,n$.

A log-log plot is employed to separate the regions associated with normal and anomalous observations on the complementary cumulative probability curve. In this plot, anomalous observations tend to occupy an approximately linear segment. As shown in Figure 2(a)-Figure 2(d), the slope between adjacent black points is calculated to detect this approximately linear segment. This segment is shaded in gray in Figure 2(a)-Figure 2(d). The 90th percentile is further marked as a reference threshold to verify the placement of the gray region. Blue, orange, red, and green dumbbells in Figure 2(a)-Figure 2(d) represent the anomaly ranges identified by $RPS$, $IQR$, $Z_{S_n}$, and $Z_{Q_n}$, respectively, and the legend in each panel lists the corresponding criteria.

Figure 2(a)-Figure 2(d) demonstrate that Equation (2) is more effective than Equations (10)-(12) in detecting anomalous pollutant concentrations. In Figure 2(a)-Figure 2(c), the left boundaries of the gray regions closely align with the 90th percentiles, while in Figure 2(d) the discrepancy between the left boundary and the 90th percentile is negligible. These results indicate that the gray regions provide sufficiently accurate estimates of the ranges occupied by anomalous pollutant concentrations in all four plots. The blue dumbbells are consistently located close to the gray regions, supporting the effectiveness of Equation (2). In contrast, although the orange dumbbell in Figure 2(a) is also close to the gray region, the orange dumbbells in Figure 2(b)-Figure 2(d) are shorter than the widths of the corresponding gray regions, suggesting that Equation (10) yields an excessive number of false negatives. Moreover, lowering the thresholds from $Z_{S_n}>3.5$ and $Z_{Q_n}>3.5$ to $Z_{S_n}>3$ and $Z_{Q_n}>3$ could improve the agreement between the red and green dumbbells and the gray regions. This observation suggests that the original thresholds, $Z_{S_n}>3.5$ and $Z_{Q_n}>3.5$ may be overly conservative for anomaly detection.

Next, this study applies a random forest classifier to classify the collected air quality data. Before training the random forest classifier, standardizing the $RPS$ values of these air quality data is implemented; thus, they can be incorporated as anomaly-informed features. For the SO₂ concentration, its standardized $RPS$ is calculated by

$${\overline{RPS}}_{\mathrm{S}{\mathrm{O}}_2}={\displaystyle \frac{RP{S}_{\mathrm{S}{\mathrm{O}}_2}-\mathrm{Median}\left(RP{S}_{\mathrm{S}{\mathrm{O}}_2}\right)}{IO{R}_{\mathrm{S}{\mathrm{O}}_2}}}$$

(13)

where the overline denotes the standarization, the $RP{S}_{\begin{array}{c}\hfill \mathrm{S}{\mathrm{O}}_2\end{array}}$ represents the $RPS$ value derived from the SO₂ concentration, $\mathrm{Median}(·)$ and $IQ{R}_{\begin{array}{c}\hfill \mathrm{S}{\mathrm{O}}_2\end{array}}$ are its median and $IQR$; respectively. Applying the same procedures yields ${\overline{RPS}}_{\begin{array}{c}\hfill \mathrm{C}\mathrm{O}\end{array}}$, ${\overline{RPS}}_{\begin{array}{c}\hfill \mathrm{N}{\mathrm{O}}_2\end{array}}$, and ${\overline{RPS}}_{P{M}_{2.5}}$.

In Taiwan, the AQI is used to indicate whether air quality may be unhealthy. Based on the AQI regulation [27], the official six-category AQI is consolidated into three health-oriented categories. The class of good to acceptable contains the official categories of good and ordinary air quality are combined into a class of good to acceptable, whereas the class of unhealthy is composed of the categories of unhealthy, very healthy, and harmful. Accordingly, this study defines the class label Y in the form as

$$Y=\left\{\begin{array}{c}0\mathrm{AQI}\le 100(\mathrm{Good}\mathrm{to}\mathrm{acceptable})\hfill \\ 1101\le \mathrm{AQI}\le 150(\mathrm{Unhealthy}\mathrm{for}\mathrm{sensitive}\mathrm{acceptable}\mathrm{groups})\hfill \\ 2\mathrm{AQI}\ge 151\left(\mathrm{Unhealthy}\right)\hfill \end{array}\right.$$

(14)

where Equation (14) yields 22,534 instances for $Y=0$, 1,367 instances for $Y=1$, and 106 instances for $Y=2$, indicating substantial class imbalance.

The random forest classifier is subsequently trained on the collected air quality data enriched with standardized $RPS$-based features. Within this framework, the standardized $RPS$ variables do not determine or alter the class labels Y. Instead, they function as auxiliary predictors that characterize the abnormality intensity of individual pollutants. Accordingly, the proposed framework is more appropriately viewed as an AQI-category classification approach with anomaly-informed feature augmentation, rather than an anomaly-based relabeling scheme.

Unless otherwise specified, the random forest classifier is implemented with the scikit-learn package. Default hyperparameter settings are employed, except that the maximum depth of the tree is set to 8 and the class weight is set to ’balanced’. The latter setting is used to address class imbalance automatically. Specifically, scikit-learn assigns class weights inversely proportional to the empirical class frequencies in the training data. Thus, minority-class observations receive greater weights during model fitting, whereas majority-class observations receive lower weights, reducing the model’s tendency to favor majority classes.

The collected air quality data are randomly partitioned into training and testing subsets, with 75% of the instances used for model training and the remaining 25% reserved for performance evaluation. Because both the train-test split and random forest training involve randomness, results from a single run of the random forest classifier may be unstable. Therefore, it is independently trained using different random seeds. Nineteen random seeds are selected from the integer range of 42 to 65.

Across 19 runs, the standard deviation of each macro- or weighted-average classification metric, such as the macro F1-score, is below 5% of its corresponding mean, indicating that the performance estimates are sufficiently stable. For each run, the confusion matrix is row-normalized by dividing every entry in a row by the sum of that row. Thus, each row represents the distribution of predicted classes conditional on the true class.

Equation (15) and Table 4 present the mean and standard deviation of the normalized confusion matrices and classification metrics obtained across 19 repeated runs. Each classification metric is reported as mean ± standard deviation. Because Equation (14) describes a multi-class classification problem, the AUC is computed using a one-vs-rest approach in which each class is treated as the positive class in turn, and all remaining classes are treated as the negative class. The resulting class-specific AUC values are then aggregated to obtain an overall multi-class AUC value. For binary classification, the accuracy, precision, recall, specificity, and F1 scores are defined by ${\displaystyle \frac{\mathrm{True}\mathrm{positive}+\mathrm{True}\mathrm{negative}}{\mathrm{All}\mathrm{samples}}}$, ${\displaystyle \frac{\mathrm{True}\mathrm{positive}}{\mathrm{True}\mathrm{positive}+\mathrm{False}\mathrm{positive}}}$, ${\displaystyle \frac{\mathrm{True}\mathrm{positive}}{\mathrm{True}\mathrm{positive}+\mathrm{False}\mathrm{negative}}}$, ${\displaystyle \frac{\mathrm{True}\mathrm{negative}}{\mathrm{True}\mathrm{negative}+\mathrm{False}\mathrm{positive}}}$, and ${\displaystyle \frac{2\times \mathrm{Precision}\times \mathrm{Recall}}{\mathrm{Precision}+\mathrm{Recall}}}$; respectively. In the multi-class setting, these metrics are calculated for each class using the one-vs-rest framework and further macro- or weighted-averaged.

$$\begin{array}{cc}\hfill & \begin{array}{cc}& \begin{array}{ccc}\mathrm{AQI}\le 100& 100<\mathrm{AQI}\le 150& \mathrm{AQI}\ge 151\end{array}\\ \begin{array}{c}\mathrm{AQI}\le 100\\ 100<\mathrm{AQI}\le 150\\ \mathrm{AQI}\ge 151\end{array}& \left[\begin{array}{ccc}0.843\pm 0.0072& 0.101\pm 0.011& 0.056\pm 0.0084\\ 0.232\pm 0.022& 0.597\pm 0.022& 0.171\pm 0.0270\\ 0.169\pm 0.083& 0.426\pm 0.129& 0.405\pm 0.0500\end{array}\right]\end{array}(\mathrm{Raw}\mathrm{data})\hfill \\ \hfill & \begin{array}{cc}& \begin{array}{ccc}\mathrm{AQI}\le 100& 100<\mathrm{AQI}\le 150& \mathrm{AQI}\ge 151\end{array}\\ \begin{array}{c}\mathrm{AQI}\le 100\\ 100<\mathrm{AQI}\le 150\\ \mathrm{AQI}\ge 151\end{array}& \left[\begin{array}{ccc}0.848\pm 0.0108& 0.097\pm 0.016& 0.055\pm 0.0094\\ 0.224\pm 0.0250& 0.603\pm 0.0450& 0.173\pm 0.028\\ 0.127\pm 0.0360& 0.348\pm 0.0660& 0.525\pm 0.035\end{array}\right]\end{array}\begin{array}{c}(\mathrm{Raw}\mathrm{data}\\ \mathrm{and}RPS\mathrm{values})\end{array}\hfill \end{array}$$

(15)

Accordingly, Table 4 and Equation (15) show that augmenting raw air quality data with standardized $RPS$-based features improves the performance of the random forest classifier on the collected air quality dataset. The improvement is reflected simultaneously in the macro-, weighted-average classification metrics reported in Table 4, and in the class-specific recall values shown in Equation (15). The improved recall values may imply that the classification of 133($\approx (0.88-0.843)\times 22,534+(0.603-0.597)\times 1,367+(0.525-0.405)\times 106$) instances has been corrected after the augmentation. Since the air quality was monitored daily, these 133 instances denote the measurements of air quality over 133 days.

Specifically, the augmented random forest classifier achieves modest but consistent improvements in macro-averaged classification performance, while maintaining competitive weighted-average performance. More importantly, the recall for the most severe pollution class, AQI≥151, increases from 0.405 ± 0.05 to 0.525 ± 0.035. This result indicates that incorporating $RPS$-based features substantially improves the model’s ability to identify observations with anomalous AQI values, which are typically rarer and more difficult to classify in a heavy-tailed dataset.

Surprisingly, the augmented features also improve the stability of predictions for the class $Y=2$. Without $RPS$-based feature augmentation, the standard deviation of the recall for $Y=2$ (AQI ≥ 151) is 12.3% of the corresponding mean, whereas it decreases to 6.71% of the corresponding mean after augmentation. This finding suggests that the $RPS$-based representation provides more robust discriminatory information for identifying severe air-pollution events.

However, the improvement for $Y=2$ is accompanied by a reduction in the probability of classifying these observations as the intermediate class, $100<AQI\le 150$, from 0.426 ± 0.129 to 0.348 ± 0.066. This shift is consistent with the improved diagonal value for the severe-pollution class and indicates that the augmented model better distinguishes the two higher-AQI categories. In contrast, the classification performance of the low-AQI and intermediate-AQI classes remains broadly stable, with only small changes in their diagonal entries.

4.2. Reservoir Water Quality

In Taiwan, reservoirs serve as a primary source of drinking water. Therefore, continuously monitoring water quality is essential for protecting public health. The second real dataset comprises 6,518 valid monthly water-quality records collected from 66 reservoirs between 2016 and 2025 [28]. The dataset includes conductivity, dissolved oxygen, turbidity, ammonia nitrogen (NH₃−N), and Carlson’s Trophic State Index (CTSI). Dissolved oxygen is reported as a percentage, whereas conductivity, turbidity, and NH₃−N are measured in $\mu \mathrm{S}/\mathrm{cm}$, NTU, and mg/L, respectively. Figure 3 presents the distributions of these variables, and Table 5 summarizes their descriptive statistics.

The resulting ${\gamma }_{Hill}$ values in Table 5 indicate five heavy-tailed distributions. Applying Equation (2) to anomaly detection for conductivity, dissolved oxygen, and turbidity introduces an additional methodological challenge. Specifically, anomalous observations may occur in both the upper and lower tails of the distributions; therefore, two approximately linear segments may need to be identified on the corresponding complementary cumulative probability curves. Such bidirectional anomalies are physically meaningful in water-quality monitoring, where dissolved oxygen and turbidity are commonly monitored indicators. For example, exceptionally high dissolved oxygen concentrations may be associated with water-pollution events, whereas abnormally low dissolved oxygen concentrations may indicate eutrophication. Eutrophic water may be unsuitable for drinking-water use.

Equation (16) presents the Pearson correlation of conductivity, dissolved oxygen, turbidity, and ammonia nitrogen (NH₃−N). Inconsistent with the findings in Section 4.1, no pair of variables in this subsection exhibits a strong correlation.

$$\begin{array}{cc}& \begin{array}{ccc}\begin{array}{c}\mathrm{Conduc}\text{-}\\ \mathrm{tivity}\end{array}& \begin{array}{c}\mathrm{Dissolved}\\ \mathrm{oxygen}\end{array}& \begin{array}{c}\mathrm{Turbi}\text{-}\\ \mathrm{dity}\end{array}{\mathrm{NH}}_3\text{-}\mathrm{N}\end{array}\\ \begin{array}{c}\mathrm{Conductivity}\\ \mathrm{Dissolved}\mathrm{oxygen}\\ \mathrm{Turbidity}\\ {\mathrm{NH}}_3\text{-}\mathrm{N}\end{array}& \left[\begin{array}{cccc}1& 0.034& 0.132& 0.169\\ 0.034& 1& -0.031& -0.230\\ 0.132& -0.031& 1& 0.083\\ 0.169& -0.230& 0.083& 1\end{array}\right]\end{array}$$

(16)

Figure 4(a)-Figure 4(d) compare the anomalous ranges of current conductivity, dissolved oxygen, turbidity, and NH₃−N data identified using the $RPS$, $IQR$, $Z_{S_n}$, $Z_{Q_n}$, and the 90th, 10th, and 5th percentile rules. As in Figure 2(a)-Figure 2(d), the black points represent the complementary cumulative probability curves derived from the reservoir water-quality observations. Blue, orange, red, and green dumbbells indicate the ranges of anomalies identified by the $RPS$, $IQR$, $Z_{S_n}$, $Z_{Q_n}$, respectively. Dashed lines denote the locations of 90th, 10th, and 5th percentiles. Gray regions represent approximately linear segments detected by calculating the slopes between successive black points.

Different from Section 4.1, conductivity, turbidity, and dissolved oxygen may exhibit extremely low values in this section. Therefore, after identifying the gray regions, dashed vertical lines are added to Figure 4(a)-Figure 4(c) to indicate the locations of the 5th or 10th percentiles. The selected percentile threshold, either 5th or 10th, is determined according to the extent of the corresponding gray region.

The need to use the 5th, 10th, and 90th percentile thresholds to verify the gray regions in Figure 4(a)-Figure 4(d) illustrates the limitations of applying percentile rules to anomaly detection in an ad hoc manner. The selection of percentile thresholds must be guided by the physical interpretation and distributional properties of each water-quality variable. In particular, the potential occurrence of unusually low or high values should be considered before defining one-sided or two-sided percentile-based anomaly criteria.

In contrast, calculating the $IQR$, $Z_{S_n}$, and $Z_{Q_n}$ metrics fails to identify extremely low values of conductivity and turbidity in Figure 4(a) and Figure 4(c). Moreover, the ranges of anomalies represented by the orange, red, and green dumbbells in Figure 4(a)-Figure 4(b) are narrower than the gray regions and blue dumbbells, which are created using the $RPS$. Only the results provided by the $RPS$, $IQR$, $Z_{S_n}$, and $Z_{Q_n}$ in Figure 4(d) may be acceptable, albeit imperfect. These results suggest that the $RPS$ metric is more effective than the $IQR$, $Z_{S_n}$, and $Z_{Q_n}$ in detecting anomalies in the collected reservoir water quality data.

The Ministry of Environment in Taiwan manages reservoir water quality using the CTSI. It is calculated from Secchi-disk transparency, chlorophyll-a concentration, and total phosphorus concentration. These three variables are excluded from the training data used to train a random forest classifier because their inclusion would result in target leakage. Specifically, the CTSI is directly calculated from these measurements; therefore, using them as training data would allow the classifier to reconstruct the class label rather than identify reservoir trophic states from independent water-quality information.

According to the official standard, the class label Y is defined as follows[28,29]:

$$Y=\left\{\begin{array}{c}0\mathrm{CTSI}<40(\mathrm{Oligotrophic}\mathrm{reservoir}\mathrm{water})\hfill \\ 140\le \mathrm{CTSI}\le 50(\mathrm{Mesotrophic}\mathrm{reservoir}\mathrm{water})\hfill \\ 2\mathrm{CTSI}>50(\mathrm{Eutrophic}\mathrm{reservoir}\mathrm{water})\hfill \end{array}\right.$$

(17)

Based on Equation (17), the numbers of instances in classes CTSI < 40, 40 ≤ CTSI ≤ 50, and CTSI > 50 are 740, 3,314, and 2,464, respectively. Unlike the classification of air quality data in Section 4.1, the present classification of reservoir water quality data did not exhibit a substantial class-imbalance problem.

As in Equation (13), the $RPS$ values of conductivity, dissolved oxygen, turbidity, and NH₃−N are standardized using their respective medians and Tukey interquartile ranges as anomaly-informed features. The random forest classifier is then trained using the collected reservoir water-quality data augmented with these standardized $RPS$ values. The data are randomly divided into training and test sets, with 75% of the observations used for model training and the remaining 25% reserved for performance evaluation. As in Section 4.1, the random forest classifier is trained independently under different random seeds. Eighteen seeds are selected from the integer range of 30 to 80.

Across 18 runs, the standard deviation of each macro-averaged and weighted-averaged classification metric is less than 5% of its corresponding mean, indicating that the estimated classification performance is sufficiently stable. As in Section 4.1, the confusion matrix for each run is row-normalized by dividing each entry in a row by the sum of that row. Consequently, each row represents the conditional distribution of predicted classes given the true class.

Equation (18) and Table 6 report the means and standard deviations of the row-normalized confusion matrices and classification metrics across the repeated 18 runs. As in Section 4.1, each classification metric is reported as the mean ± standard deviation. Multi-class AUC values are calculated using a one-versus-rest approach in which each class is treated in turn as the positive class, and all remaining classes are treated as the negative class. The resulting class-specific AUC values are then aggregated to obtain an overall multi-class AUC. Likewise, class-specific precision, recall, and F1-score are calculated under the one-versus-rest framework and subsequently summarized using macro and weighted averages.

$$\begin{array}{cc}\hfill & \begin{array}{cc}& \begin{array}{ccc}\mathrm{CTSI}<40& 40\le \mathrm{CTSI}\le 50& \mathrm{CTSI}>50\end{array}\\ \begin{array}{c}\mathrm{CTSI}<40\\ 40\le \mathrm{CTSI}\le 50\\ \mathrm{CTSI}>50\end{array}& \left[\begin{array}{ccc}0.671\pm 0.0021& 0.327\pm 0.021& 0.002\pm 0.0042\\ 0.018\pm 0.0041& 0.930\pm 0.011& 0.052\pm 0.0100\\ 0& 0.171\pm 0.012& 0.829\pm 0.0120\end{array}\right]\end{array}(\mathrm{Raw}\mathrm{data})\hfill \\ \hfill & \begin{array}{cc}& \begin{array}{ccc}\mathrm{CTSI}<40& 40\le \mathrm{CTSI}\le 50& \mathrm{CTSI}>50\end{array}\\ \begin{array}{c}\mathrm{CTSI}<40\\ 40\le \mathrm{CTSI}\le 50\\ \mathrm{CTSI}>50\end{array}& \left[\begin{array}{ccc}0.709\pm 0.0096& 0.290\pm 0.0093& 0.001\pm 0.0018\\ 0.015\pm 0.0060& 0.939\pm 0.0100& 0.046\pm 0.0080\\ 0& 0.170\pm 0.0970& 0.830\pm 0.0960\end{array}\right]\end{array}\begin{array}{c}(\mathrm{Raw}\mathrm{data}\mathrm{and}\\ RPS\mathrm{values})\end{array}\hfill \end{array}$$

(18)

Table 6 and Equation (18) illustrate a pattern that contrasts with that observed in Section 4.1. In the present classification of reservoir water quality data, the instances with $\mathrm{CTSI}>50$, representing reservoir water that is unsuitable as a source of drinking water, are not the least frequent class. Accordingly, their recall value is not the lowest among the three classes. In contrast, the class with $\mathrm{CTSI}<40$ contained the fewest instances and exhibited the lowest recall value.

Overall, although the improvements are modest, augmenting the reservoir water-quality data with $RPS$-based features improved the macro-averaged classification performance and the recall for each class. The improved recall values may imply that the classification of approximately 60 instances has been corrected after the augmentation. Since the reservoir water quality is monitored monthly [28], those 60($\approx (0.709-0.671)\times 740+(0.939-0.930)\times 3,314+(0.83-0.829)\times 2,464$) instances are the measurements of reservoir water quality over 60 months. If further improvements in classification performance are desired, combining the results in Section 4.1 and Section 4.2 suggests that tuning the maximum tree depth is preferred.

4.3. River Water Quality

On the same website [28] used to collect the reservoir water-quality data for implementing Section 4.2, river water-quality data are also available. Although the dataset includes additional variables, such as dissolved oxygen and NH₃−N, this study retained the air temperature, potential of hydrogen (pH), conductivity, coliform groups, and river pollution index (RPI). A total of 19,699 records were collected for this section. Including other available variables is avoided because they either introduced target leakage or reduced classification accuracy.

The collected river water-quality data consist of monthly measurements from 296 monitoring stations across 54 rivers in Taiwan, spanning 2016 to 2025. Air temperature, conductivity, and coliform groups are measured in °C, $\mu \mathrm{S}/\mathrm{cm}$, and CFU/100 mL, respectively. Figure 5 and Table 7 present the distributions and descriptive statistics of these variables.

The ${\gamma }_{Hill}$ values report in Table 7 suggest that all five variables exhibit heavy-tailed characteristics. Furthermore, the skewness values indicate two distinct distributional patterns. Air temperature and pH have left-skewed distributions, whereas conductivity, coliform groups, and RPI have right-skewed distributions.

Equation (19) presents the Pearson correlation matrix for air temperature, pH, conductivity, and the coliform group. Consistent with the findings in Section 4.2, no pair of variables in this subsection exhibit a strong correlation.

$$\begin{array}{cc}& \begin{array}{cccc}\begin{array}{c}\mathrm{Air}\mathrm{tem}\text{-}\\ \mathrm{perature}\end{array}& \mathrm{pH}& \begin{array}{c}\mathrm{Con}\text{-}\\ \mathrm{ductivity}\end{array}& \begin{array}{c}\mathrm{Coliform}\\ \mathrm{groups}\end{array}\end{array}\\ \begin{array}{c}\mathrm{Air}\mathrm{temperature}\\ \mathrm{pH}\\ \mathrm{Conductivity}\\ \mathrm{Coliform}\mathrm{group}\end{array}& \left[\begin{array}{cccc}1& 0.1330& 0.0113& 0.0084\\ 0.1330& 1& -0.0123& -0.0574\\ 0.0113& -0.0123& 1& 0.0016\\ 0.0084& -0.0574& 0.0016& 1\end{array}\right]\end{array}$$

(19)

Figure 6(a)-Figure 6(d) compare the anomalous ranges identified for successive observations of air temperature, pH, conductivity, and coliform groups using the $RPS$, $IQR$, $Z_{S_n}$, $Z_{Q_n}$, 90th and 5th percentile rules. As in Figure 2(a)-Figure 2(d) and Figure 4(a)-Figure 4(d), the black points represent complementary cumulative distribution curves derived from the river water-quality observations. The blue, orange, red, and green dumbbells represent the anomalous ranges identified using the $RPS$, $IQR$, $Z_{S_n}$, and $Z_{Q_n}$ methods, respectively. The dashed lines indicate the 90th and 5th percentiles. The gray regions represent approximately linear segments identified from the slopes between successive black points. Unlike Section 4.2, the 10th percentiles are not shown in Figure 6(a)-Figure 6(d), because the resulting ranges differ substantially from those identified by the other methods.

Figure 6(a)-Figure 6(d), consistent with Figure 4(a)-Figure 4(d), illustrate a limitation of percentile-based rules. Expert judgment is required to determine whether the 90th percentile rule, the 5th percentile rule, or both rules should be applied for anomaly detection. For example, without a thorough understanding of the river system and its normal pH variability, it is difficult to select appropriate percentile thresholds for identifying anomalous pH observations. In contrast, the anomaly ranges identified by the $RPS$ are generally close to the approximately linear regions shown in gray in Figure 6(a)-Figure 6(d), with only minor discrepancies.

In contrast, the consistency of the $IQR$, $Z_{S_n}$, and $Z_{Q_n}$ with the gray regions is less evident. The discrepancies between the orange, red, and green dumbbells and the gray regions are pronounced in Figure 6(a) and Figure 6(c), whereas the discrepancies in Figure 6(d) appear comparatively acceptable. Therefore, when these methods are applied to detect anomalies in new river water-quality data, their resulting anomaly ranges may require cautious interpretation.

To compare the classification performance of the collected river water-quality data with and without $RPS$-based data augmentation, this study formulates a binary classification problem. Based on the official RPI classification criteria, the class variable Y is defined by

$$Y=\left\{\begin{array}{c}0\mathrm{if}RPI\le 2\hfill \\ 1\mathrm{otherwise}\hfill \end{array}\right.$$

(20)

Based on Equation(20), the RPI ≤ 2 and RPI > 2 classes contain 11,385 and 8,584 observations, respectively, corresponding to 57.0% and 43.0% of the collected river water-quality dataset. Thus, this dataset exhibits only a mild class-imbalance problem.

As in Section 4.1 and Section 4.2, the $RPS$ values for air temperature, pH, conductivity, and the coliform group are standardized using their respective medians and Tukey interquartile ranges to construct anomaly-informed features. A random forest classifier is then trained using the collected river water-quality data augmented with these standardized $RPS$ features. These data is randomly divided into training and test sets, with 75% of the observations used for training and the remaining 25% reserved for performance evaluation. As in Section 4.1, the classifier is trained independently under different random seeds. Sixteen seeds are selected from the integer range of 20 to 50.

Across the 16 runs, the standard deviation of each macro-averaged and weighted-average classification metric is less than 5 of its corresponding mean, indicating that the estimated classification performance is stable across random seeds. As in Section 4.1 and Section 4.2, the confusion matrix for each run is row-normalized by dividing each cell by the sum of the entries in its corresponding row. Consequently, each row represents the conditional distribution of predicted classes given the true class.

Equation (21) and Table 6 summarize the confusion matrices and classification metrics across the 16 repeated runs. The reported values are expressed as the mean ± standard deviation.

$$\begin{array}{l}\begin{array}{cc}& \begin{array}{cc}\mathrm{RPI}\le 2\hfill & \mathrm{RPI}>2\hfill \end{array}\\ \begin{array}{c}\mathrm{RPI}\le 2\\ \mathrm{RPI}>2\end{array}& \left[\begin{array}{cc}0.770\pm 0.0075& 0.230\pm 0.0075\\ 0.003\pm 0.0011& 0.997\pm 0.0011\end{array}\right](\mathrm{Raw}\mathrm{data})\end{array}\hfill \\ \begin{array}{cc}& \begin{array}{cc}\mathrm{RPI}\le 2& \mathrm{RPI}>2\end{array}\\ \begin{array}{c}\mathrm{RPI}\le 2\\ \mathrm{RPI}>2\end{array}& \left[\begin{array}{cc}0.778\pm 0.0036& 0.222\pm 0.0036\\ 0.001\pm 0.0010& 0.999\pm 0.0010\end{array}\right](\mathrm{Raw}\mathrm{data}\mathrm{and}RPS\mathrm{values})\end{array}\hfill \end{array}$$

(21)

Accordingly, Equation (21) and Table 8 show that augmenting the collected river-water-quality data with their $RPS$ values improves both macro-level and class-specific classification performance. The corresponding recall improvements suggest that approximately 91($\approx (0.778-0.770)\times 11,385$) additional instances in the $\mathrm{RPI}\le 2$ class, and approximately 25($\approx (0.999-0.996)\times 8,584$) additional instances in the $\mathrm{RPI}>2$ class were correctly classified after augmentation. Although the total of approximately 116 corrected instances may appear modest, each instance represents a distinct monthly measurement. Thus, the improvement corresponds to correctly classifying water-quality conditions in 116 different months over the monitoring period.

However, we may be interested in the relatively low recall value for the class $\mathrm{RPI}\le 2$. Inspecting the $RPS$ values of the misclassified instances finds that some of them are anomalous, but some of them are not. To enhance the classification accuracy, the current classification criteria [28] for defining the RPI may be insufficient.

4.4. Sensitivity Study

Training a random forest classifier requires the specification of several hyperparameters. Their effects on predictive performance must therefore be examined through a sensitivity analysis. Following the last paragraph of Section 3.2, this study evaluates four key hyperparameters: the maximum tree depth (max_depth), the number of trees in the forest (n_estimators), the minimum number of samples required to split an internal node (min_samples_split), and the minimum number of samples required at a leaf node (min_samples_leaf). Like Section 4.2, 18 random seeds are selected from the integer range of 30 to 80. For each seed, the reservoir water-quality data are randomly divided into training and test sets, with 75% of these data used for training and the remaining 25% reserved for testing. Random forest classifiers are trained with and without the inclusion of standardized $RPS$ features. Figure 7(a)-Figure 7(d) show the AUC values obtained under different hyperparameter settings. Each point represents the mean AUC across the 18 random seeds, whereas the error bars indicate the corresponding standard deviations.

Figure 7(a)-Figure 7(d) indicate that the AUC is most sensitive to the maximum tree depth (max_depth). In contrast, variations in the other hyperparameters have only modest effects on AUC performance. Furthermore, augmenting the reservoir water-quality data with $RPS$-based features consistently improves the AUC across the evaluated hyperparameter settings.

5. Discussion

In light of the contributions stated in Section 1 and the numerical results presented in Section 4.1 and Section 4.2, the findings of this study are discussed as follows:

• The computation of $RPS$ in Equation (2) involves relatively few steps: calculating $Pr\left(t_{obs}\right)$ and identifying ${Pr}_{max}\left(t_{ij}\right);i=1,2,\dots ,N;j=1,2,\dots ,n$. Therefore, its computational procedure is more straightforward than those required for calculating Rousseeuw–Croux estimators $Z_{Q_n}$ and $Z_{S_n}$ in Equations (11) and (12). This computational simplicity enhances the reproducibility and practical accessibility of the proposed measure.
• Augmenting the collected sensor data with $RPS$-based features improves both macro-level and class-specific classification performance of a random forest classifier. Section 4.1, Section 4.2 and Section 4.3 show this improvement across three environmental sensor datasets subjected to different degrees of class imbalance. Although the improved macro-level and class-specific classification performance may be modest for both normal and anomalous sensor readings, it should be noted that these sensor readings were not collected over a short time period. Therefore, the enhanced macro-level and class-specific classification performance may reduce the risk of misinterpreting long-term environmental conditions. Surprisingly, the $RPS$-based features enabled the random forest classifier to identify nearly all anomalous river-pollution events in Section 4.3. These results suggest that the proposed augmentation can provide useful evidence for environmental policymaking and regional governance.
• Calculating the $RPS$ robustly detects anomalous environmental measurements in Section 4.1, Section 4.2 and Section 4.3. Despite the accuracy of detected anomalies, calculating the Tukey interquartile range $IQR$ and Rousseeuw–Croux estimators $Z_{Q_n}$ and $Z_{S_n}$ fails to detect extremely low conductivity, turbidity (in Section 4.2), and anomalous high air temperature (in Section 4.3). If we don’t understand the physical meaning of an environmental sensor measurement, it is difficult to determine which percentile rule should be imposed for locating the range of anomalous environmental sensor readings.
• Section 4.3 highlights the importance of explainability in evaluating classification performance. In this section, the recall for the class RPI ≤ 2 is lower than that for the class RPI > 2, suggesting that misclassified instances are more likely to occur in the former class. We further examine these misclassified cases to improve classification accuracy. However, because the official regulation does not define the river pollution index for the dataset used in this section, it is difficult to explain the misclassifications without augmenting the feature set with $RPS$-based variables. Figure 2(a)–Figure 2(d), Figure 4(a)–Figure 4(d), and Figure 6(a)–Figure 6(d) show that $RPS$ may be a more effective anomaly measure for identifying the anomalous range of a variable than $IQR$, $Z_{Q_n}$, and $Z_{S_n}$, thereby providing a more interpretable basis for explaining misclassified instances.

6. Conclusions

This study applies a random forest classifier to classify air quality, reservoir water quality, and river water quality sensor data augmented with $RPS$-based features. The $RPS$ serves as an anomaly measure to facilitate a random forest classifier in separating normal and anomalous data. To ensure the $RPS$ is an effective anomaly measure, additional tests show that it is a simpler and effective anomaly measure in comparison with the $IQR$, $Z_{S_n}$, and $Z_{Q_n}$. Furthermore, the $RPS$ improves the explainability of classification results. Referring to the corresponding $RPS$ values, we can understand why an instance is misclassified. Since we must avoid the target leakage in defining a class variable, raw data can’t provide similar explainability of misclassification results.

Conclusively, the contributions of this study are improving the weakness of a random forest classifier in classifying heavy-tailed data, demonstrating that the $RPS$ is a simple and effective anomaly measure, and enhancing the explainability of classification results provided by a random forest classifier. Based on these contributions, this study concludes:

• The $RPS$ provides a computational simple anomaly measure for identifying anomalous sensor readings. Its calculation requires fewer steps than the calculations of $Z_{Q_n}$ and $Z_{S_n}$. This computational simplicity enhances the reproducibility and practical accessibility of the proposed measure.
• The $RPS$ serves as an alternative anomaly measure for environmental monitoring data in which conventional mean- and standard-deviation-based measures may be sensitive to extreme values. Numerical results of this study show that it can identify anomalous observations in heavy-tailed environmental sensor datasets.
• Augmenting heavy-tailed environmental datasets with $RPS$-based features improves the performance of a random forest classifier. The improvement was observed in both macro-level and class-specific classification results. After augmentation, the classifier made fewer classification errors for both normal and anomalous sensor readings, implying the reduce of misinterpreting long-time environmental conditions.
• The $(RPS$-based features can enhance the identification of rare but environmentally important events. The collected three datasets for implementing this study are subjected to different degrees of class imbalance. However, augmenting with $RPS$-based features improves the accuracy in classifying all three datasets. In particular, the augmented classifier identified nearly all anomalous river-pollution events.
• The collected sensor data represented long-term environmental monitoring records. The proposed approach. Therefore, it may support environmental assessment, routine data screening, and evidence-based environmental policymaking.

However, the study may have several limitations. First, the collected data for implementing this study were sourced from environmental sensors. The generalizability of $RPS$ must be further examined using data from other types of sensors, monitoring locations, and sampling frequencies. Second, this study uses a random forest classifier to generate classification results. Future studies may evaluate the usefulness of $RPS$-based features in other machine learning, such as the Naive-Bayes classifier. Investigating their use in real-time environmental monitoring systems is also suggested.