A Bayesian Approach with Type-2 Student-t Membership Function for T-S Model Identification

TSK fuzzy model provides a rule-based structure for modeling a complex non-linear system. If $G(x,y)$ is the system to be identified, where $\mathbf{x}\in R^{n}$ be the input vector and $y\in R$ be the output. Then, the $i^{th}$ rule is written as

Rule i : IF $x_{i}$ is $A^{i}_{1}$ and $\cdots$ and $x_{m}$ is $A^{i}_{m}$ THEN

where, $i=1,\cdots,c$ is the number of fuzzy rule and $y^{i}$ is the $i^{th}$ output. Using these rules, we can infer the final model output as follows:

where, $w^{i}$ denotes the overall firing strength of the $i^{th}$ rule.

In the interval type-2 FCM, two different objective function that differ in their degrees of fuzziness are optimized simultaneously using Lagrange multipliers to obtain the upper and lower membership function [3]. Let $m_{1}$ and $m_{2}$ be the two degree of fuzziness, then the two objective function are described as

The main motivation of inter type-2 fuzzy c-regression algorithm is to partition the set of $n$ data points $(\mathbf{x}_{k},y_{k})$ $(k=1\cdots n)$ into c clusters. The data points in every cluster $i$ can be described by a regression model as

where, $\mathbf{x}_{k}=[x_{k1},\cdots,x_{km}]$ be the $k^{th}$ input vector, $j=1,\cdots,m$ be the number of features, $i=1,\cdots,c$ be the number of clusters and $\mathbf{\zeta}_{i}=[b^{i}_{1},\cdots,b^{i}_{m},b^{i}_{0}]$ be the coefficient vector of the $i^{th}$ cluster. In [4], the coefficient vectors are optimized by weighted least square method, whereas, in our approach we used SGD. The primary objective for using SGD is to make the algorithm robust even for the cases where $[\mathbf{x}^{T}\mathbf{P}_{i}\mathbf{x}]$ become singular [6].

Firstly, we formulate the task of fuzzy-space partitioning as a Maximum A-Posterior (MAP) over a squared error function [8]. Exploiting Bayes rule the MAP estimator is defined as

where, $p(x/y)$ be the posterior, $p(y/x)$ be the likelihood and $p(x)$ be the prior distribution. Using the above equation the MAP estimator is expressed in term of regression problem as

where, ${E_{ik}}(\zeta_{i})$ be the MAP estimator, $(y_{k}-g^{i}(x_{k},{\zeta_{i}}))^{2}$ be the likelihood, $\sum_{p=0}^{m}(b^{i}_{p})^{2}$ be the prior or called as a regularizer, which is equivalent to the Bayesian notion of having a prior on the regression weights $b^{i}$ for each cluster $i$ and $\lambda$ be the regularizer control parameter. The regularizer reduces the overfitting in cluster assignment by constraining the small regression weights.

In the proposed approach, we first define two degrees of fuzziness $m_{1}$ and $m_{2}$, initializes the number of clusters $c$ and a termination threshold $\epsilon$. We also initialize the parameters $\overline{\zeta_{i}}$ and $\underline{\zeta_{i}}$, which are the upper and lower regression coefficient vectors of $i^{th}$ cluster. Then, the equation (7) is written in the term upper and lower error function MAP estimator as follows:

To reduce the complexity of the system, a weighted average type reduction technique is applied to obtain ${E_{ik}}(\zeta_{i})$ as

Through MAP estimate on the posterior of defuzzified error function ${E_{ik}}(\zeta_{i})$, the upper and lower membership function for every data points in each cluster are obtained as similar to [3] and they are given as follows:

The above equation can be interpreted as that for a MAP problem formulated in (8). To estimate the parameters $\overline{\zeta_{i}}$ and $\underline{\zeta_{i}}$, we formulate the problem as a locally weighted linear regression with an objective function:

Here, $u_{ik}$ denotes the membership value of $k^{th}$ data point in the $i^{th}$ cluster. The parameter $\overline{\zeta_{i}}$ and $\underline{\zeta_{i}}$ are estimated by SGD using the above objective function by appropriately finding $\overline{u}_{ik}$ and $\underline{u}_{ik}$. Then the regression coefficient ($\zeta_{i}$) are obtained by a type reduction technique as follow:

The steps in this subsection are run and the parameters are updated until the convergence of $||\zeta_{i}^{current}-\zeta_{i}^{previous}||$ $\geq$ $\epsilon$ to obtain the optimal value of the regression coefficient as briefly described in Algorithm 1.

The MF developed in [5] is hyperplane shaped, which cannot successfully incorporate the relevant information of data distributions within different clusters. To overcome this issue, we proposed a modified Gaussian based MF combined with a student-t density function. The student-t distribution is widely used as a prior in Bayesian inference for sparse data modelling [7]. Here, we weigh the Gaussian and the student-t part by a hyper-parameter $\alpha$. If the data we are modelling is very sparse then, $\alpha$ should be set very low so as to give more weight to student-t density membership value.

In the above, $d_{ik}$ is the distance between $k^{th}$ input vector and $i^{th}$ cluster hyperplane.

where, $r=\max\{d_{ik}(\overline{\zeta}_{i}),\;\;i=1\cdots c\}$ is the maximum distance of $k^{th}$ input vector from the $i^{th}$ cluster, $v_{i}$ and $\sigma_{i}$ denotes the average distance and variance of each data points from the cluster hyperplane respectively.

The lower MF ($\underline{\mu}_{A_{i}}(\mathbf{x}_{k})$) and upper MF ($\overline{\mu}_{A_{i}}(\mathbf{x}_{k})$) are called as weights of the TSK fuzzy model corresponding to $k^{th}$ input belonging to the $i^{th}$ cluster.

In the most of the literature the defuzzification of weights is computed before determining the model output $\hat{y}_{k}$. The problem with these approaches are that they do not consider effect of model output which will affect the over all performance of the model. To overcome this problem we evaluated the $\underline{y}_{k}$ and $\overline{y}_{k}$ corresponding to the $\underline{\mu}_{A_{i}}(\mathbf{x}_{k})$ and $\overline{\mu}_{A_{i}}(\mathbf{x}_{k})$ using the KM algorithm [2]. The values of $\underline{y}_{k}$ and $\overline{y}_{k}$ are optimized parallelly until the convergence. The another advantage of this approach is that it become more robust in handling noise and provide a confidence interval for every output data points. The model output $\underline{y}_{k}$ and $\overline{y}_{k}$ corresponding to the weights $\underline{\mu}_{A_{i}}(\mathbf{x}_{k})$ and $\overline{\mu}_{A_{i}}(\mathbf{x}_{k})$ are calculate using (1) and (2) as follows:

where, $p$ and $q$ are switching points and computed by KM algorithm. We run above mentioned steps until the convergence of $\overline{y}_{k}$ and $\underline{y}_{k}$. Finally , the model output is determined by applying a type reduction technique as

The house prices dataset (https://www.kaggle.com/lespin/house-prices-dataset) is used to predict the sale price of a particular property. Through experimentation, we have demonstrated the robustness of proposed method on this sparse data. The dataset is divided in training (70$\%$) and testing (30$\%$) sets and five-fold cross-validation is used while training. The hyper-parameters of the model are initially set as: $c=3$, $m_{1}=1.6$, $m_{2}=4.7$, $\lambda=0.3$, $\alpha=0.15$ and $\eta=3.7$. It should be noted that the value of $\alpha=0.15$ is small because the dataset is sparse. So, in MF, the contribution of student-$t$ function should be high, which is ensured by a smaller value of $\alpha$ i.e., larger value of $1-\alpha$ as defined by (12) and (13). The mean square error (MSE) is $0.008$ on the test data, which is lower than state-of-the-art methods as shown in Table I. The absolute value of error as shown in Fig. 2 is also small in compare to the absolute house prices as shown in Fig. 1. We postulate that this is due to the student-$t$ MF used in our model, which helps in robustly quantifying the effects of sparse data. Also, the higher test accuracy is due to greater generalization owing to $L_{2}$ regularizer used in our model. The coefficient of variation which is the ratio of explained variance to total variance is very high ($0.85$). This suggests that our model captures variations in the data robustly and is not susceptible to faulty performance in the presence of outliers.

The second-order non-linear difference equation as given in (20) is used in order to draw comparison with other benchmark models as given in Table II.

where, $v(k)=\sin(2k/25)$ is the input for validation of model, $z(k)$ is the model output whereas, $z(k-1)$, $z(k-2)$ and $u(k)$ are the model inputs respectively. The hyper-parameters are tuned by grid search and finally set as: $c=4$, $m_{1}=1.5$, $m_{2}=7$ and $\eta=3.14$. The obtained MSE of the model on $500$ test data points is $7.2\times 10^{-5}$ using only four rules which is much smaller compared to other models. Through simulations, we have shown that proposed model outperforms with other state-of-the-art model. The Fig. 3 shows that our model output closely tracks the actual output at every time-step. As observed in Fig. 4, the error fluctuates with data point, but the absolute error is consistently less than 0.1 with no rapid surge at stationary points of time series data. This is a crucial requirement for a stable system. Therefore, we conclude that our algorithm yields a dynamically stable model.

In this subsection a non-linear sinc function is used to present the effectiveness of the proposed model;

where, $x\in[-40,0)\bigcup(0,40]$. We have sampled 121 data points uniformly for this one dimensional function. As similar to previous case study, the number of rules is taken as four. The hyper-parameters are tuned through grid-search and finally fixed as: $m_{1}=1.5$, $m_{2}=7$ and $\eta=3.14$. The MSE of the proposed model is $2.4\times 10^{-3}$, which is lower in compare to modified inter type-2 FRCM (MIT2-FCRM) [6], which is $7.7\times 10^{-3}$ on the test data of 121 samples. The Table III provides a detailed comparison of performance with state-of-the-art methods.