Dynamics and density function for a stochastic anthrax epidemic model

1.   Introduction

With the rapid development of big data and artificial intelligence, data-driven predictions now permeate nearly every aspect of modern science ^[1,2]. Regression analysis, a key component of data-driven science, has wide applications across various domains, including new material design, stock prediction, medical diagnosis, and geological exploration, etc. ^{[3,4,5,6,7,8,9]}. Besides, machine learning offers numerous methods for regression analysis, such as support vector regression (SVR), decision trees (DTs), multi-layer perceptrons (MLPs), deep neural networks, etc. Although those methods can achieve high accuracy, the "black box" nature of neural networks limits their interpretability.

According to the Weierstrass approximation theorem, continuous functions defined on closed intervals can be uniformly approximated by polynomial functions ^[10]. Moreover, a polynomial function has a simpler form compared to other complex functions, making polynomial regression a natural choice for practical applications ^[11,12]. Polynomial regression can be divided into univariate function fitting and multi-variate function regression. While univariate function fitting is highly effective supported by popular data analysis software such as Origin, MATLAB, and SPSS based on the least squares algorithm, there remains a lack of fast, effective, and stable algorithms for multivariate function regression.

The theme of this paper is to design a fast, efficient and accurate algorithm to solve the multivariate polynomial regression problems (Eq (1.1)).

In the regression process, SVR introduces a kernel function to construct the nonlinear relation. The MLP method employs the weighted summation of variables ($\sum_{i}\omega_{ji}x_{i}+\theta_{j}$) to achieve the aggregation of information. As the weighted summation only passes a linear relationship between variables, the nonlinear knowledge, for example the product relation $x_{j}^{\mu_{j}}x_{k}^{\nu_{k}}$, should be approximated by using an activation function, such as the sigmoid function, $f(x) = 1/(1+e^{-x})$, the tanh function, $f(x) = (e^{x}-e^{-x})/(e^{x}+e^{-x})$, etc. However, the applications of kernel functions in SVR and activation functions in MLPs make it difficult to interpret the machine learning model.

Recently, Kolmogorov-Arnold networks (KANs) have gained attention due to their superior accuracy and interpretability compared to multi-layer perceptrons ^[13]. MLPs have fixed activation functions on nodes while KANs have learnable activation functions on edges. The Kolmogorov-Arnold representation theorem established that, if $f$ is a multivariate continuous function on a bounded domain, then $f$ can be written as a finite composition of continuous functions of a single variable and the binary operation of addition. The regression model of KANs is in the form of

where, $L$ represents the number of network layers and $\phi_{k, i_{k+1}, i_{k}}$ is a continuous functions of signal variable. That is to say, KANs train and select suitable signal variable continuous functions in each network layer from the given basic function library including polynomial functions, exponential functions, logarithmic functions, etc. Although, KANs have higher accuracy and interpretability, the regression of product relation in multivariate polynomial functions still need to transform to the coupling of univariate functions, which will undoubtedly increase the complexity of the regression.

The traditional artificial neural network uses nodes whose information transmission is via a sigma unit, i.e., a linear weighted sum of inputs, $a = \sum_{i}\omega_{i}x_{i}$. A pi unit is constructed by replacing the weighted summation in the activation with a weighted product, $a = \prod_{i}x_{i}^{\omega_{i}}$. Units of this type are also designated "higher order" in that they contain polynomial terms of order greater than one (linear). The simplest way of modeling such a product relation in multivariate polynomial functions is to increase the complexity of the node by a pi unit. Moreover, a sigma-pi unit can formulate a sum of product terms,

Inspired by this architecture of a sigma-pi unit, we find the high-order sigma-pi neural network (SPNN) ^[14,15] can efficiently achieve the regression of multivariate polynomial functions. The coefficients and exponents of the multivariate polynomial can be determined through the trained neural network weight parameters ($c_{j}$ and $\omega_{ji}$).

Performing the regression of multivariate polynomial functions by SPNN has the following advantages compared with that of MLPs and KANs.

● Polynomial regression by the SPNN model has better interpretability than that of MLPs, due to the fact that high-order terms are generated by pi units instead of complex activation functions.

● Polynomial regression by the SPNN model has faster learning efficiency compared with that of KANs, owing to the application of a fixed activation function in polynomial form.

● The optimization and control of the parameters in the SPNN model are more convenient than that in MLPs and KANs. Besides, the forms of the SPNN model are more friendly to real applications.

2.   Method and modeling

2.1.   The high-order sigma-pi neural network

In an MLP neural network, information from different input signals is transmitted through weighted summation. For example, in a two-layer forward neural network (see Figure 1(a)), the hidden layer

and the output signal

where $\sigma_1$ is the activation function of the hidden layer, $\sigma_2$ is the activation function of the output layer, $\alpha$ is the weight of the first layer network, $\beta$ is the weight of the second layer network, and $b$ and $\theta$ are the biases. The information transmission from $x_{j}$ to $h_{i}$, and from $h_{i}$ to $y$, are all by weighted summation. So, this kind of structure unit in MLP neural networks can be called sigma-sigma units.

To embed high-order coupled information into a neural network, high-order neural networks are designed. For instance, the sigma-pi neural network (SPNN) enables the embedding of product information into the network. In SPNN, the hidden layer

and the output layer

The information transmission from $x_{j}$ to $g_{i}$ is by weighted product, and from $g_{i}$ to $y$ is by weighted summation, see Figure 1(b).

2.2.   Multivariate polynomial regression

As previously mentioned, the basic unit of SPNN is the power function of input signals. We aim to achieve multivariate polynomial regression using the explainable SPNN. In this model, if the activation functions, $\sigma_{1}$ and $\sigma_{2}$, are linear, the output signal

Cleary, Eq (2.5) is a multivariate polynomial function, where $\beta_{i}$ is the polynomial coefficient, $\alpha_{ij}$ is the polynomial exponent, and $\sum_{i = 1}^{m}\beta_{i}b_{i}+\theta$ is the constant term. Thus, this explainable SPNN is essentially a type of Kolmogorov-Arnold network.

To illustrate the regression process of a multivariate polynomial function more clearly, we present a simple regression problem as an example:

Here, $\{x_{1}, x_{2}, \cdots, x_{n}\}$ represents the state variables, and $y$ is the target variable. Since Eq (2.6) lacks a constant term, bias parameters are not used in the designed network. A two-layer SPNN is employed to solve this regression problem. The number of neurons in the input layer equals $n$, while the hidden layer contains $n+n(n-1) = n^2$ neurons. The neural connections from the input layer to the hidden layer primarily capture the regression of power functions, such as $x_{i}^{b_{i}}$ and $x_{i}^{p_{i}}x_{j}^{q_{j}}$, corresponding to each term in the polynomial function. The network then performs a weighted summation of all terms from the hidden layer to the output layer.

Additionally, to accelerate the parameter optimization, we convert the above sigma-pi neural network to the sigma-sigma neural network, and use the back-propagation algorithm to update the parameters. The input signals are first transformed using logarithmic functions,

Next, a sigma unit operates on the $u_{i}$,

Here $\alpha_{ji}$ represents the weight from the $i^{th}$ neuron in the input layer to the $j^{th}$ neuron in the hidden layer. A signed exponential function is then applied to $v_{j}$,

Therefore, the conversion from the pi unit to the sigma unit can be achieved using logarithmic and exponential operators. The output signal, $y$, is then obtained by the sigma unit from the hidden layer to the output layer.

Here, $\beta_{j}$ is the weight from the $j^{th}$ neuron in the hidden layer to the output signal. The model error is defined as

where $y$ is the model output and $\hat{y}$ is the real value. For parameter updates, the gradient descent method is employed,

Here, $\eta$ is the learn rate,

and

The optimized parameters $\beta$ and $\alpha$ are the corresponding multivariate polynomial regression coefficients and exponents.

2.3.   Algorithm pseudo code

The pseudo code of multivariate polynomial regression by our designed neural network is presented in Algorithm 1 below.

3.   Examples and application

In this section, we present two examples to demonstrate the feasibility and effectiveness of the proposed model. Moreover, we illustrate the application of our method in analyzing multivariate correlated regression, specifically concerning the maximum stress on concrete pipes under traffic loads.

3.1.   Example 1: Integer order multivariate polynomial regression

In this example, the model is used to determine the regression coefficients and exponents of Eq (3.1) based on generated data. The variable $x_{1}$ is linearly sampled from the interval $[1, 3]$, $x_{2}$ is randomly sampled from the interval $[1, 4]$, and $y$ is calculated according to Eq (3.1). These generated data are briefly presented in Figure 2, where Figure 2(a) gives the time series of $x_{1}$, $x_{2}$, and $y$; Figure 2(b), (c) display scatter plots of $y$ versus $x_{1}$, and $y$ versus $x_{2}$, respectively.

The number of neurons in the hidden layer equals 3, i.e., the number of terms of polynomial function Eq (3.1). Let $u = \log (x_{1}), v = \log (x_{2})$, and the signals $u$ and $v$ are fed into the sigma unit as inputs. It should be noted that the variables $x_{1}$ and $x_{2}$ in this example are restricted to positive values. We use the linear activation function in each architecture, thus, the neurons in the hidden layer satisfy

So, by the exponential operation,

The output signal is the weighted sum of $e^{h_{1}}, e^{h_{2}}$, and $e^{h_{3}}$ according to the sigma unit from the hidden layer to the output layer.

Here, $\beta_{i}$ represents the corresponding polynomial regression coefficient, and $\alpha_{ij}$ denotes the related polynomial regression exponent. Since we are addressing an integer-order polynomial regression problem, the parameters $\alpha_{ij}$ are constrained to integer values during the parameter iteration process. The training is conducted for 500 epochs, with a maximum learning rate of 0.01. Compared with other optimizers, such as BGD, SGD, and MBGD, the combination of the Adam optimizer with a one-cycle scheduler has faster convergence speed and better modeling performance, so we chose this optimization algorithm in the training process. The regression performance is evaluated using the mean squared error (MSE) and the coefficient of determination ($R^2$),

Here, $y_{i}$ is the real value, $\hat{y_{i}}$ is the regression value, and $\bar{y}$ is the mean value of $y$. The modeling results are shown in Figure 3, the train loss and valid loss are depicted by the MSE, and the coefficient of determination $R^{2} = 1.0$.

The weight parameters of the proposed neural network are listed in Table 1. The parameter $\alpha_{ij}$ is related to the polynomial exponent, and parameter $\beta_{i}$ corresponds to the polynomial coefficient.

From the results above, it is evident that our designed SPNN effectively solves the integer-order multivariate polynomial regression problem. The coefficients and exponents of the integer-order multivariate polynomial function align with the weight parameters of the designed SPNN. In the next example, we use the SPNN to solve a fractional-order multivariate polynomial regression problem.

3.2.   Example 2: Fractional-order multivariate polynomial regression

In this example, we use the designed SPNN to solve the fractional-order polynomial regression problem. The values of $x_{1}$ and $x_{2}$ are generated in the same manner as in Example 1, and the output signal is computed according to Eq (3.3). The dataset contains 3000 samples, and the training is conducted for 10,000 epochs. The Adam optimizer is employed combined with a one-cycle scheduler. The maximum learning rate is set to 0.009, while the optimizer's learning rate is initialized at 0.004. The corresponding time series are plotted in Figure 4(a), and the regression result is shown in Figure 4(b).

The values of mean squared error (MSE), coefficient of determination ($R^2$), and corresponding weight parameters are listed in Table 2.

From the regression results, we observe small discrepancies between the actual polynomial coefficients (or exponents) and the weight parameters. However, the parameters closely align with the polynomial coefficients and exponents when considering their approximate values. Thus, we conclude that the designed SPNN can effectively solve fractional-order multivariate polynomial regression problems.

3.3.   Application: Multi-variable regression of the maximum stress

In this part, we apply the designed SPNN to analyze the maximum stress of concrete sewage pipelines under the combined influence of traffic load, earth pressure, and groundwater level, with a consideration of 12 physical parameter variables. These variables encompass corrosion depth ($Cd$), corrosion width ($Cw$), corrosion length ($Cl$), void width ($Vw$), void length ($Vl$), burial depth ($H$), traffic load ($P$), pipe diameter ($D$), wall thickness ($t$), bedding modulus ($Eb$), backfill soil modulus ($Es$), and groundwater level ($hw$). The detailed parameter ranges for each variable are provided in Table 3.

These physical variables are randomly generated within the specified ranges, resulting in a dataset of 250 samples. The maximum stress signals are formulated by finite element simulation based on these datasets. Figure 5 illustrates the dependence of the maximum stress on each physical variable. As the maximum stress of concrete sewage pipelines is determined by the 12 physical variables, it is challenging to analyze the relationship between individual variables and stress through simple data fitting. As shown in Figure 5, the relationship between maximum stress and each variable is nonlinear and complex, indicating that this is a multivariate nonlinear regression problem. Existing software and algorithms are unable to efficiently and accurately address this challenge; however, our designed SPNN method is capable of handling it effectively.

According to the univariate empirical formula ^[7,9], we define the multivariate regression equation of the maximum stress as

The data set is first normalized to eliminate the influence of dimensionality, and then randomly split into training, validation, and test sets, with respective ratios of $64\%$, $20\%$, and $16\%$. It makes almost no difference to divide the data set at the other ratios, for example, $70\%$, $15\%$, and $15\%$, because of the excellent model performance. The regression model is developed using the explainable SPNN. During the training process, a total of 1000 epochs are set, and the Adam optimizer with with a learning rate equal to 0.006 is employed together with a one-cycle scheduler with the maximum learn rate equal to 0.01. Figure 6(a) shows the training loss and valid loss depicted by MSE, and Figure 6(b) presents a comparison between the real and predicted maximum stress values.

The regression error MSE $= 3.36\times 10^{-5}$, and the coefficient of determination $R^2 = 0.9998$. We save model parameters to two decimal places, and the regression equation is

It should be noted that the values of coefficients and exponents in Eq (3.4) may exhibit minor fluctuations due to the broad feasible region of the solutions. However, the contributions of each physical variable to the maximum stress derived from Eq (3.4) can be well interpreted. It is evident that the corrosion depth, corrosion width, corrosion length, void width, void length, burial depth, traffic load, pipe diameter, wall thickness, and bedding modulus have positive effects on the maximum stress, while the backfill soil modulus and groundwater level exhibit negative effects. These findings align with existing empirical knowledge.

Additionally, the SHAP approach ^[16,17,18] is utilized to analyze the impact of input features on the model's output. The SHAP value for each sample ($x$) corresponding to feature ($f_{i}$) is computed using Eq (3.5).

where $f_{i}$ represents the $i^{th}$ feature (i.e., the $i^{th}$ physical variable), $F$ is the total number of features, and $f$ is a feature subset that includes $f_{i}$. $|f|$ denotes the number of elements in subset $f$, and $\binom{F}{|f|}$ is the number of combinations of $F$ items taken $|f|$ at a time. $P_{f}(x)$ is the predicted value when all the physical variables in $f$ are selected, and $P_{f\setminus f_{i}}(x)$ is the predicted value with all the variables in $f$ except $f_{i}$. The SHAP values for the samples related to each physical variable are shown in Figure 7.

In Figure 7, if the higher feature value corresponds to larger SHAP values, it indicates a positive effect of that feature on the model output. Conversely, if the higher feature value results in smaller SHAP values, it suggests a negative effect. Thus, we conclude that features $P, D, t, Cd, Vl, Eb, Cw, Cl, H$, and $Vw$ have a positive effect on the maximum stress, while $Es$ and $hw/H$ have a negative effect. These conclusions align with the regression analysis results from Eq (3.4). Additionally, the importance ranking of the input features is $P > D > t > Cd > Vl > Eb > Cw > Cl > H > Vw > Es > hw/H$.

According to the regression results, we can further analyze the evolution of the maximum stress of concrete sewage pipelines with respect to the univariate input. By selecting a sample data point and employing the control variable method, the maximum stress signal can be calculated using Eq (3.4), with only one variable being varied. Figure 8 shows the evolution of maximum stress with respect to each physical variable. These univariate evolutionary trends provide theoretical guidance for strategies aimed at repairing and improving the properties of concrete sewage pipelines.

4.   Conclusions

In this manuscript, we propose an explainable sigma-pi neural network to address the multi-variable nonlinear polynomial regression problem. The coefficient and exponent parameters of polynomials are effectively represented by the corresponding weight parameters in the SPNN. To accelerate the regression process, the back-propagation algorithm is employed for parameter optimization. The examples reveal that our designed SPNN can efficiently and accurately solve both integer-order and fractional-order multivariate polynomial regression problems. In practical application, the SPNN provides high-precision fitting of the maximum stress in concrete sewage pipelines under the combined influence of 12 physical parameter variables. Furthermore, feature importance ranking and additional analyses on the relationship between maximum stress and these variables can be conducted based on this explainable machine learning model.

Building upon the framework for solving multi-variable polynomial regression, future algorithms for addressing more complex nonlinear regression problems could be developed through higher-order neural networks, such as sigma-pi or sigma-pi-sigma networks. This study lays the theoretical foundation for developing generalized tools to solve multi-variable nonlinear regression problems.

Use of AI tools declaration

The authors declare that this manuscript is the authors' original work, and they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by the China Postdoctoral Science Foundation Funded Project (Grant No. 2022M712902), the Natural Science Foundation of Henan (Grant No. 232300421345), and the Zhongyuan Youth Top Talent Plan (Zhongyuan Youth Postdoctoral Innovative Talents).

Conflict of interest

The authors declare no conflict of interest with respect to the research, authorship, and/or publication of this article.