Marginal likelihood based model comparison in Fuzzy Bayesian Learning

TL;DR

This paper extends Fuzzy Bayesian Learning by using marginal likelihood for model comparison, providing a Bayesian alternative to error minimization, validated through synthetic and real-world financial data.

Contribution

It introduces a Bayesian model evidence approach for selecting the best fuzzy rule base within Fuzzy Bayesian Learning, enhancing model comparison capabilities.

Findings

Marginal likelihood effectively distinguishes between competing fuzzy models.

Bayesian model selection outperforms mean squared error minimization in case studies.

Validated approach on synthetic and real financial data.

Abstract

In a recent paper [1] we introduced the Fuzzy Bayesian Learning (FBL) paradigm where expert opinions can be encoded in the form of fuzzy rule bases and the hyper-parameters of the fuzzy sets can be learned from data using a Bayesian approach. The present paper extends this work for selecting the most appropriate rule base among a set of competing alternatives, which best explains the data, by calculating the model evidence or marginal likelihood. We explain why this is an attractive alternative over simply minimizing a mean squared error metric of prediction and show the validity of the proposition using synthetic examples and a real world case study in the financial services sector.

Tables

Table 1. TABLE I: Comparison of model evidences for all the synthetic dataset case

	Single Ellipsoidal		Multi Ellipsoidal
Model	$log\left(Z\right)$	Fcalls	$log\left(Z\right)$	Fcalls
$H_{true}$	-72.087 $\pm$ 0.665	27387	-71.880 $\pm$ 0.657	13257
$H_1$	-738.822 $\pm$ 0.714	9127	-738.699 $\pm$ 0.710	7495
$H_2$	-460.352 $\pm$ 0.672	14862	-460.362 $\pm$ 0.676	8586
$H_3$	-71.867 $\pm$ 0.658	11631	-72.878 $\pm$ 0.677	11631
$H_{GLM1}$	-3564.282 $\pm$ 0.598	1904	-3564.282 $\pm$ 0.598	1904
$H_{GLM2}$	-1330.995 $\pm$ 0.623	2510	-1330.995 $\pm$ 0.623	2510
$H_{GLM3}$	-1547.996 $\pm$ 0.779	6978	-1548.437 $\pm$ 0.782	5220
$H_{GLM4}$	-1269.334 $\pm$ 0.909	7529	-1268.904 $\pm$ 0.903	10612
$H_{GLM5}$	-694.100 $\pm$ 1.130	26739	-693.865 $\pm$ 1.128	63324

Table 2. TABLE II: Rule bases for the different Fuzzy Bayesian models fitted on the real world case

Model	Antecedents			Consequent
	O&M	Loss history	DnTnF	System uptime
$H_{rw1}$	Good	-	-	High
	Good	-	Average	Medium
	Average	-	-	High
	Bad	Good	-	High
	Bad	Bad	Bad	Low
$H_{rw2}$	-	Good	-	High
	-	Good	Average	Medium
	-	Average	-	High
	Good	Bad	-	Low
	Bad	Bad	Bad	Low
$H_{rw3}$	-	-	Good	High
	Average	-	Good	Medium
	-	-	Average	Medium
	-	Good	Bad	High
	Bad	Bad	Bad	Low

Table 3. TABLE III: Comparison of model evidences for the real world example: σ 𝜎 \sigma estimated from parameters

		Single Ellipsoidal		Multi Ellipsoidal
Model	Prior	$log\left(Z\right)$	Fcalls	$log\left(Z\right)$	Fcalls
$H_{rw1}$	[0,10]	-59.314 $\pm$ 0.350	1464	-59.314 $\pm$ 0.350	1464
$H_{rw2}$	[0,10]	-62.412 $\pm$ 0.281	969	-62.412 $\pm$ 0.281	969
$H_{rw3}$	[0,10]	-66.469 $\pm$ 0.223	504	-66.469 $\pm$ 0.223	504
$H_{GLM6}$	[-100,100]	-74.591 $\pm$ 0.642	3533	-73.679 $\pm$ 0.630	4584
	[-10,10]	-65.368 $\pm$ 0.489	1898	-65.339 $\pm$ 0.486	1849
$H_{GLM7}$	[-100,100]	-95.391 $\pm$ 0.922	10420	-94.500 $\pm$ 0.907	18289
	[-10,10]	-79.761 $\pm$ 0.731	4500	-78.925 $\pm$ 0.719	5876
$H_{GLM8}$	[-100,100]	-109.288 $\pm$ 1.048	53316	-109.722 $\pm$ 1.053	39006
	[-10,10]	-86.613 $\pm$ 0.804	19700	-87.005 $\pm$ 0.807	18399

Table 4. TABLE IV: Comparison of model evidences for the real world example: σ = 0.25 𝜎 0.25 \sigma=0.25

		Single Ellipsoidal		Multi Ellipsoidal
Model	Prior	$log\left(Z\right)$	Fcalls	$log\left(Z\right)$	Fcalls
$H_{rw1}$	[0,10]	-596.783 $\pm$ 0.673	12454	-598.993 $\pm$ 0.691	6392
$H_{rw2}$	[0,10]	-780.824 $\pm$ 0.765	28171	-780.058 $\pm$ 0.745	16852
$H_{rw3}$	[0,10]	-1903.106 $\pm$ 0.642	10605	-1904.153 $\pm$ 0.661	5648
$H_{GLM6}$	[-100,100]	-795.069 $\pm$ 0.735	5564	-795.367 $\pm$ 0.741	4115
	[-10,10]	-787.166 $\pm$ 0.620	2222	-786.612 $\pm$ 0.612	2358
$H_{GLM7}$	[-100,100]	-643.082 $\pm$ 1.077	18943	-643.200 $\pm$ 1.073	20175
	[-10,10]	-652.957 $\pm$ 0.918	8948	-652.871 $\pm$ 0.919	9109
$H_{GLM8}$	[-100,100]	-589.990 $\pm$ 1.237	546259	-589.281 $\pm$ 1.237	1066203
	[-10,10]	-598.993 $\pm$ 1.060	21860	-597.475 $\pm$ 1.050	65359

Table 5. TABLE V: Comparison of model evidences for the real world example: σ = 1.0 𝜎 1.0 \sigma=1.0

		Single Ellipsoidal		Multi Ellipsoidal
Model	Prior	$log\left(Z\right)$	Fcalls	$log\left(Z\right)$	Fcalls
$H_{rw1}$	[0,10]	-73.250 $\pm$ 0.449	5451	-73.575 $\pm$ 0.463	3295
$H_{rw2}$	[0,10]	-88.356 $\pm$ 0.500	13774	-88.637 $\pm$ 0.527	8447
$H_{rw3}$	[0,10]	-150.194 $\pm$ 0.330	1443	-150.636 $\pm$ 0.323	929
$H_{GLM6}$	[-100,100]	-94.155 $\pm$ 0.659	3888	-94.727 $\pm$ 0.666	2966
	[-10,10]	-86.649 $\pm$ 0.532	1697	-85.996 $\pm$ 0.518	1596
$H_{GLM7}$	[-100,100]	-110.170 $\pm$ 0.972	9579	-110.018 $\pm$ 0.970	11597
	[-10,10]	-95.198 $\pm$ 0.775	6355	-94.820 $\pm$ 0.771	5959
$H_{GLM8}$	[-100,100]	-123.227 $\pm$ 1.130	50185	-123.128 $\pm$ 1.118	57867
	[-10,10]	-102.192 $\pm$ 0.895	13508	-100.428 $\pm$ 0.876	32887

Table 6. TABLE VI: Posterior parameter estimates for the synthetic simulation cases, mean (standard deviation)

Model	Method	theta0	theta1	theta2	theta3	theta4	theta5	theta6	theta7	theta8
$H_{true}$	Single Ellip	4.98 (0.12)	4.99 (0.21)	5.02 (0.07)	5.02 (0.21)	4.97 (0.21)	4.95 (0.26)	49.22 (4.46)	50.06 (0.98)	49.85 (2.11)
	Multi Ellip	4.98 (0.12)	4.98 (0.22)	5.03 (0.07)	5.04 (0.23)	4.95 (0.22)	4.97 (0.25)	49.14 (4.64)	49.96 (1.00)	49.74 (2.34)
$H_1$	Single Ellip	0.48 (0.28)	0.14 (0.13)	3.66 (0.09)	9.95 (0.04)	9.13 (0.12)	3.68 (2.07)	3.16 (0.29)	18.64 (1.16)	72.28 (1.30)
	Multi Ellip	0.49 (0.29)	0.13 (0.12)	3.68 (0.09)	9.95 (0.05)	9.16 (0.15)	3.79 (2.02)	3.14 (0.29)	18.71 (1.19)	72.32 (1.42)
$H_2$	Single Ellip	5.08 (2.81)	3.92 (0.30)	4.63 (0.11)	7.50 (0.34)	4.81 (0.18)	4.98 (2.89)	3.90 (0.27)	43.26 (1.90)	52.13 (1.96)
	Multi Ellip	4.95 (2.86)	3.91 (0.27)	4.63 (0.10)	7.46 (0.37)	4.80 (0.17)	5.04 (2.78)	3.91 (0.27)	43.42 (2.01)	51.94 (1.96)
$H_3$	Single Ellip	4.98 (0.11)	4.97 (0.23)	5.03 (0.07)	5.05 (0.22)	4.95 (0.20)	4.98 (0.26)	49.21 (4.53)	50.07 (1.05)	49.78 (2.10)
	Multi Ellip	4.99 (0.12)	4.98 (0.21)	5.03 (0.07)	5.04 (0.21)	4.94 (0.19)	4.96 (0.25)	49.31 (4.42)	50.08 (1.01)	50.03 (2.10)
$H_{GLM1}$	Single Ellip	49.99 (0.01)	0.12 (0.00)	-	-	-	-	-	-	-
	Multi Ellip	49.99 (0.01)	0.12 (0.00)	-	-	-	-	-	-	-
$H_{GLM2}$	Single Ellip	49.94 (0.06)	2.83 (0.04)	-1.91 (0.04)	-	-	-	-	-	-
	Multi Ellip	49.94 (0.06)	2.83 (0.04)	-1.91 (0.04)	-	-	-	-	-	-
$H_{GLM3}$	Single Ellip	49.98 (0.02)	0.24 (0.01)	-0.18 (0.01)	0.11 (0.02)	-	-	-	-	-
	Multi Ellip	49.98 (0.02)	0.24 (0.01)	-0.18 (0.01)	0.11 (0.02)	-	-	-	-	-
$H_{GLM4}$	Single Ellip	49.97 (0.03)	3.69 (0.16)	-3.62 (0.17)	-0.11 (0.02)	0.12 (0.02)	0.13 (0.02)	-	-	-
	Multi Ellip	49.97 (0.03)	3.67 (0.16)	-3.61 (0.16)	-0.11 (0.02)	0.12 (0.02)	0.13 (0.02)	-	-	-
$H_{GLM5}$	Single Ellip	49.92 (0.08)	16.70 (0.38)	-10.86 (0.44)	-3.17 (0.10)	1.73 (0.11)	-0.50 (0.06)	0.18 (0.01)	-0.10 (0.01)	0.09 (0.01)
	Multi Ellip	49.92 (0.08)	16.64 (0.38)	-10.82 (0.40)	-3.16 (0.10)	1.71 (0.10)	-0.49 (0.06)	0.18 (0.01)	-0.09 (0.01)	0.09 (0.01)

Equations

$$p\left(H_i\mid \mathbf{D}\right)=\frac{p\left(\mathbf{D}\mid H_i\right)p\left(H_i\right)}{p\left(\mathbf{D}\right)}$$

$$p\left(\theta \mid \mathbf{D},H_i\right)=\frac{p\left(\mathbf{D}\mid \theta ,H_i\right)p\left(\theta \mid H_i\right)}{p\left(\mathbf{D}\mid H_i\right)}$$

$${\mathcal{Z}}_i=p\left(\mathbf{D}\mid H_i\right)=\int p\left(\mathbf{D}\mid \theta ,H_i\right)p\left(\theta \mid H_i\right)d\theta$$

$$\frac{p\left(H_2\mid \mathbf{D}\right)}{p\left(H_1\mid \mathbf{D}\right)}=\frac{p\left(\mathbf{D}\mid H_2\right)p\left(H_2\right)}{p\left(\mathbf{D}\mid H_1\right)p\left(H_1\right)}=\frac{{\mathcal{Z}}_{2}p\left(H_2\right)}{{\mathcal{Z}}_{1}p\left(H_1\right)}$$

$$p\left(\theta \mid H\right)d\theta =\pi \left(\theta \right)d\theta =dX$$

$$p\left(\mathbf{D}\mid \theta ,H\right)=L\left(\theta \right)=L\left(X\right)$$

$$\mathcal{Z}=\int L\left(\theta \right)\pi \left(\theta \right)d\theta ={\int }_0^{1}L\left(X\right)dX$$

$$X\left(\lambda \right)={\int }_{L\left(\theta \right)\ge \lambda }\pi \left(\theta \right)d\theta ,X\left(0\right)=1,X\left(\infty \right)=0$$

$$\begin{array}[]{cc}\theta_{n}\sim p\left(\theta\mid\lambda\right)=\frac{\Theta[L(\theta)-\lambda]}{X(\lambda)}\pi(\theta)\\ \text{where }\Theta(x)=\left\{\begin{array}[]{ll}0;x<0\\ 1;x\geq 0\\ \end{array}\right.\end{array}$$

$$H(\lambda \to {\lambda }^{\prime })=\int p(\theta \mid {\lambda }^{\prime })ln\left[\frac{p(\theta \mid {\lambda }^{\prime })}{p(\theta \mid \lambda )}\right]d\theta =ln\left[\frac{X(\lambda )}{X({\lambda }^{\prime })}\right]$$

$$⟨H(\lambda \to {\lambda }^{\prime })⟩=-⟨lnt{⟩}_{t\sim Beta(N,1)}=1/N$$

$${\beta }_k{R}_k:\text{if }{A}_1^k\oplus A_2^k\oplus \dots \oplus A_{T_k}^k\text{ then }{C}_k$$

$${\mathbf{Y}}_N\sim \mathcal{M}\mathcal{V}\mathcal{N}\left(\mathbf{\mu },\Sigma \right)\text{,}$$

$$p\left({\mathbf{y}}_N\mid g(\mathbf{x};\theta ),{\mathbf{X}}_N\right)$$

$$=\frac{1}{\sqrt{(2\pi {)}^N{\sigma }^{2N}}}\times$$

$$\displaystyle\phantom{=}\exp\left(-\frac{1}{2\sigma^{2N}}\big{(}\boldsymbol{y}_{N}-g(\boldsymbol{X}_{N};\theta)\big{)}^{\mathrm{T}}\big{(}\boldsymbol{y}_{N}-g(\boldsymbol{X}_{N};\theta)\big{)}\right)\text{.}$$

$$f_{i}\left(u,\phi_{i}\right)=\left\{\begin{array}[]{cl}\frac{\left(u-\phi_{i_{1}}\right)}{\left(\phi_{i_{2}}-\phi_{i_{1}}\right)}&\text{if }\phi_{i_{1}}\leq u\leq\phi_{i_{2}}\\ \frac{\left(u-\phi_{i_{3}}\right)}{\left(\phi_{i_{2}}-\phi_{i_{3}}\right)}&\text{if }\phi_{i_{2}}\leq u\leq\phi_{i_{3}}\\ 0&\text{otherwise}\end{array}\right.$$

$$\begin{split}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\pgfsys@color@cmyk@stroke{1}{1}{0}{0}\pgfsys@color@cmyk@fill{1}{1}{0}{0}R_{1}}:if&\>loc\_risk==HI\;\vee\;maintenance==POOR\\ &then\;downtime==HI\\ {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\pgfsys@color@cmyk@stroke{1}{1}{0}{0}\pgfsys@color@cmyk@fill{1}{1}{0}{0}R_{2}}:if&\>loc\_risk==MED\;\vee\;maintenance==AVG\\ &then\;downtime==MED\\ {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\pgfsys@color@cmyk@stroke{1}{1}{0}{0}\pgfsys@color@cmyk@fill{1}{1}{0}{0}R_{3}}:if&\>loc\_risk==LOW\;\wedge\;maintenance==GOOD\\ &then\;downtime==LOW\\ {\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color@cmyk@fill{0}{0.89}{0.94}{0.28}R_{4}}:if&\>loc\_risk==LO\;\\ &then\;downtime==HI\\ {\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color@cmyk@fill{0}{0.89}{0.94}{0.28}R_{5}}:if&\>maintenance==POOR\\ &then\;downtime==LOW\\ {\color[rgb]{0,0.88,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.88,0}\pgfsys@color@cmyk@stroke{0.91}{0}{0.88}{0.12}\pgfsys@color@cmyk@fill{0.91}{0}{0.88}{0.12}R_{6}}:if&\>dummy\_covar==POOR\\ &then\;downtime==LOW\\ \end{split}$$

$$\begin{split}H_{true}:\{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\pgfsys@color@cmyk@stroke{1}{1}{0}{0}\pgfsys@color@cmyk@fill{1}{1}{0}{0}R_{1}},{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\pgfsys@color@cmyk@stroke{1}{1}{0}{0}\pgfsys@color@cmyk@fill{1}{1}{0}{0}R_{2}},{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\pgfsys@color@cmyk@stroke{1}{1}{0}{0}\pgfsys@color@cmyk@fill{1}{1}{0}{0}R_{3}}\}\\ H_{1}:\{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\pgfsys@color@cmyk@stroke{1}{1}{0}{0}\pgfsys@color@cmyk@fill{1}{1}{0}{0}R_{1}},{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\pgfsys@color@cmyk@stroke{1}{1}{0}{0}\pgfsys@color@cmyk@fill{1}{1}{0}{0}R_{2}},{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\pgfsys@color@cmyk@stroke{1}{1}{0}{0}\pgfsys@color@cmyk@fill{1}{1}{0}{0}R_{3}}\,{\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color@cmyk@fill{0}{0.89}{0.94}{0.28}R_{4}},{\color[rgb]{0.72,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.72,0,0}\pgfsys@color@cmyk@stroke{0}{0.89}{0.94}{0.28}\pgfsys@color@cmyk@fill{0}{0.89}{0.94}{0.28}R_{5}}\}\\ H_{2}:\{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\pgfsys@color@cmyk@stroke{1}{1}{0}{0}\pgfsys@color@cmyk@fill{1}{1}{0}{0}R_{1}},{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\pgfsys@color@cmyk@stroke{1}{1}{0}{0}\pgfsys@color@cmyk@fill{1}{1}{0}{0}R_{2}},{\color[rgb]{0,0.88,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.88,0}\pgfsys@color@cmyk@stroke{0.91}{0}{0.88}{0.12}\pgfsys@color@cmyk@fill{0.91}{0}{0.88}{0.12}R_{6}}\}\\ H_{3}:\{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\pgfsys@color@cmyk@stroke{1}{1}{0}{0}\pgfsys@color@cmyk@fill{1}{1}{0}{0}R_{1}},{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\pgfsys@color@cmyk@stroke{1}{1}{0}{0}\pgfsys@color@cmyk@fill{1}{1}{0}{0}R_{2}},{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\pgfsys@color@cmyk@stroke{1}{1}{0}{0}\pgfsys@color@cmyk@fill{1}{1}{0}{0}R_{3}},{\color[rgb]{0,0.88,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.88,0}\pgfsys@color@cmyk@stroke{0.91}{0}{0.88}{0.12}\pgfsys@color@cmyk@fill{0.91}{0}{0.88}{0.12}R_{6}}\}\\ \end{split}$$

$$GLM1:\mu ={\alpha }_0+$$

$$GLM2:\mu ={\alpha }_0+$$

$$GLM3:\mu ={\alpha }_0+$$

$$GLM4:\mu ={\alpha }_0+$$

$$GLM5:\mu ={\alpha }_0+$$

$$+$$

$$\begin{array}{rl}GLM6:\mathbf{\mu }={\alpha }_0+& \sum _{i=1}^3{\alpha }_i{x}_i\\ GLM7:\mathbf{\mu }={\alpha }_0+& \sum _{i=1}^3{\alpha }_i{x}_i+\sum _{\begin{array}{c}i,j=1\\ i\ne j\end{array}}^3{\alpha }_{i+3}{x}_i{x}_j\\ +& \sum _{\begin{array}{c}i,j,k=1\\ i\ne j\ne k\end{array}}^3{\alpha }_{i+6}{x}_i{x}_j{x}_k\\ GLM8:\mathbf{\mu }={\alpha }_0+& \sum _{i=1}^3{\alpha }_i{x}_i+\sum _{i,j=1}^3{\alpha }_{i+3}{x}_i{x}_j\\ +& \sum _{\begin{array}{c}i,j,k=1\\ i\ne j\ne k\end{array}}^3{\alpha }_{i+9}{x}_i{x}_j{x}_k\end{array}$$

Full text

I Introduction

Traditional rule based fuzzy systems encode expert opinion in the form of IF-THEN rules and optimise some performance metric (typically mean squared error in predicting a data-set) to obtain the hyper-parameters of the model (like the numeric values defining the shape of the membership functions etc.) [2, 3, 4]. The rule base is one of the core elements driving the model and slight change in the rule base would significantly affect the performance of the fuzzy inference system. Many automated methods have been proposed where the rule base structure and parameters can be automatically determined [5, 6, 7]. However interpretability of such models is an issue and various methods have been proposed to alleviate the issue [8]. In the present paper however, we are interested in the actual metric for comparison between different models having different rule bases derived from expert opinion. The comparison metric can nevertheless be embedded within an automated framework to evolve the best model if so required.

The most straight forward way of comparing the fuzzy rule bases is to optimise the model parameters based on the prediction error (e.g. mean squared error, mean absolute error etc.) on a data-set and to select the rule base with the lowest error as the better model, for the given data-set. Improvisations on this method would be to divide the data into a training and testing set and use the average prediction from k-fold cross validation on the data-set. However, as demonstrated in [9], the optimisation of an utility estimate like cross validation, with scarce data (which is the case for us), is likely to find overfitted models due to high variance in the utility estimates. The authors also show that Bayesian model averaging over all the candidate models takes the model uncertainty into account and is less prone to overfitting than cross-validation. A more generic overview of the different methods of Bayesian model comparison is presented in [10].

Another popular method is to use information criteria (IC) like Akaike’s, Deviance, Bayesian information criteria etc. for model comparison, which penalises both the prediction misfit and the complexity of the model. Now, the complexity of the model is generally penalised based on the number of effective parameters ($N_{\text{eff}}$) [11]; the rationale being that models with more free parameters would have better flexibility at fitting the data and hence more complexity. But this is not necessarily true for all class of models. For example, the fuzzy models considered in this paper have the same number of parameters to be estimated (as the number of membership functions are the same), but differ in the number and definition of rule bases. The rule bases essentially act like a complicated non-linear kernel to transform the input covariates into the output. This would be difficult to capture in $N_{\text{eff}}$ and would need customised application specific methods.

All the ICs are essentially approximations of the marginal likelihood or the Bayesian model evidence since the corresponding integral is very hard to estimate [12]. A rigorous comparison of different ways to evaluate Bayesian model evidence is shown in [13] and it is demonstrated that Bayesian evidence approximations from ICs are often heavily biased and accuracy of the model ranking is significantly affected by the choice of the approximation method. The authors also suggest using bias free numerical methods (e.g. nested sampling [14]), over ICs if it is computationally feasible. The marginal likelihood computation integrates out model parameters rather than maximising them, which automatically prevents overfitting [12]. It also maintains a trade-off between model complexity and fit to the dataset. Also the marginal likelihood can be factored into a product of predictive likelihoods which reflect the model’s ability to make predictions on new data. These reasons make the marginal likelihood an attractive metric for model selection.

In the next sections we briefly summarise the mathematical foundations of Bayesian model selection, the basic nested sampling algorithm with its improvements, integrate them with the Fuzzy Bayesian Learning framework and show demonstrative examples on synthetic and real world data-sets.

II Bayesian evidence for model selection in FBL and nested sampling

II-A Bayesian evidence and model selection using nested sampling

Bayes theorem allows for the calculation of the posterior probability of a model (or hypothesis), from a set of plausible ones, given the data and can be expressed as

[TABLE]

where $p\left(H_{i}\right)$ is the prior probability of the $i^{th}$ hypothesis ($H$) and $p\left(\boldsymbol{D}\right)$ is the probability of the data. The posterior probability of the parameters ($\theta$) of $H_{i}$ given $\boldsymbol{D}$ is

[TABLE]

where $p\left(\boldsymbol{D}\mid\theta,H_{i}\right)$ is the likelihood of the data given the hypothesis and its parameters, $p\left(\theta\mid H_{i}\right)$ is the prior probability of the parameters given the hypothesis. These two terms in the numerator are relevant for the parameter estimation problem. For the case of model selection, the important term is the denominator $p\left(\boldsymbol{D}\mid H_{i}\right)$, also known as the evidence or the marginal likelihood of the model $H_{i}$ [15]. The evidence ($\mathcal{Z}_{i}$), for the $i^{th}$ model is also the normalizing constant and can be written as

[TABLE]

The evidence automatically implements Occam’s razor [16, 17, 18], which implies that a simpler model with less number of parameters would have a higher value of evidence than a more complex model, unless the latter can explain the data much better. The model selection problem between two competing models (hypothesis) $H_{1}$, $H_{2}$ can be expressed as the ratio of how probable the respective models are in light of the given data-set $\boldsymbol{D}$ (i.e. the ratio of their posterior distributions), given by the following equation

[TABLE]

which essentially reduces to the ratio of the model evidences, if the prior probabilities of both the models are the same (i.e. $p\left(H_{2}\right)=p\left(H_{1}\right)$).

The evidence in Eq. (3) is a multi-dimensional integral and has been traditionally difficult to handle with numerical approximation techniques. The canonical method is to apply thermodynamic integration which employs a sophisticated version of MCMC sampling. However it requires large number of function evaluations (at least an order of magnitude more than parameter estimation) [19] and cannot handle phase changes properly [14]. A better way to evaluate the evidence is through the use of nested sampling [14] which overcomes the afore-mentioned problems. It employs two key concepts – the first is that the multi-dimensional evidence integral in Eqn. (3) can be transformed into a one dimensional integral. Following the lines of [14], the total prior mass is denoted by $X$ and each unit in the prior mass as $dX$ where

[TABLE]

The likelihood can then be expressed as

[TABLE]

which transforms the evidence into a one dimensional integral

[TABLE]

summing over the prior mass

[TABLE]

where $\lambda$ is the likelihood contour. The other key concept is that if $N$ particles $\theta_{1},\dots,\theta_{N}$ explore the prior above the likelihood contour $\lambda$

[TABLE]

the corresponding prior masses $X_{n}\equiv X(L_{n})$ enclosed by the contours $L_{n}\equiv L(\theta_{n})$ follow an uniform distribution. $p\left(\theta\mid\lambda\right)$ in Eqn. (9) is the prior with a likelihood constraint $\lambda$. The prior mass can be ordered as $X\left(\lambda\right)<X\left(\lambda^{\prime}\right)$ if $\lambda>\lambda^{\prime}$ and the prior mass enclosed by the worst state/sample can be predicted.

In each iteration, the nested sampling algorithm moves from the current likelihood contour ($\lambda$) to the next $\lambda^{\prime}$, where $\lambda^{\prime}>\lambda$ and the relative entropy ($H$) in Eqn. (10) can be used to measure the overlap between successive truncated priors $p\left(\theta\mid\lambda\right)$ and $p\left(\theta\mid\lambda^{\prime}\right)$ [14].

[TABLE]

The average relative entropy is constant and is dependent on the number of particles $N$, fed in by the user.

[TABLE]

The basic nested sampling algorithm pseudo-code is shown in Algorithm 1.

The key bottleneck in implementing Algorithm 1, is sampling from the truncated prior $p(\theta\mid\lambda)$. Various versions of nested sampling have evolved to solve this issue. The single ellipsoidal nested sampling [20] approximates the iso-likelihood contour by a hyper ellipsoid and works well for uni-modal posteriors, but not multi-modal ones. Multi-modal nested sampling [19] improves on this using multiple ellipsoids and a clustering algorithm to identify isolated modes in the posterior distribution. This works on low to moderate dimensional parameter spaces and other techniques relying on MCMC, like polychord [21] and diffusive nested sampling [22] have been proposed for high dimensional parameter spaces.

II-B Application to Fuzzy Bayesian Learning

Following on from the lines of [1] the inference problem is that of estimating $\theta$ for the underlying non-linear model function $g(\boldsymbol{x};\theta)$ given the data $\boldsymbol{D}$ which comprises of the set of input vectors $\boldsymbol{X}_{N}=\boldsymbol{x}_{i}\forall i\in\left\{1,2,...,N\right\}$ and the set of outputs $\boldsymbol{Y}_{N}=\boldsymbol{y}_{i}\forall i\in\left\{1,2,...,N\right\}$. As introduced in [1], the function $g(\boldsymbol{x};\theta)$ is a fuzzy inference system given by the rule base Eq. (12)

[TABLE]

where $\beta_{k}$ is a dichotomous variable indicating the inclusion of the $k^{th}$ rule in the fuzzy inference system, $A^{k}_{i}\;\forall i\in\left\{1,2,...,T_{k}\right\}$ is a referential value of the $i^{th}$ antecedent attribute in the $k^{th}$ rule, $T_{k}$ is the number of antecedent attributes used in the $k^{th}$ rule, $C_{k}$ is the consequent in the $k^{th}$ rule, $\oplus\in\left\{\vee,\wedge\right\}$ represents the set of connectives (OR, AND operations) in the rules. A referential value is a member of a referential set with the universe of discourse $U_{i}$, which are meaningful representations for specifying an attribute using subjective linguistic terms which are easily understood by an expert. Each of these referential values are characterised by membership functions $f_{i}\left(A_{i},U_{i},\phi_{i}\right)orf_{i}\left(C_{i},U_{i},\phi_{i}\right)$ and for the inference of the function $g(\boldsymbol{x};\theta)$ in Eq. (1), the set of parameters to be estimated is $\theta$, which contains all $\phi_{i}$ and $\beta_{k}$

We assume a Gaussian likelihood function given by

[TABLE]

where $\mathcal{MVN}(\boldsymbol{\mu},\Sigma)$ is the Multivariate-Normal distribution with mean vector $\boldsymbol{\mu}$ and covariance matrix $\Sigma$. The mean vector for our model comes from $g\left(\boldsymbol{X}_{N};\theta\right)$ and the fuzzy inference system from Eq. (12). We use $I_{N}$, the $N\times N$ identity matrix, along with a parameter $\sigma$ to construct the covariance matrix. This leads to the likelihood

[TABLE]

where $\boldsymbol{y}_{N}$ is a realisation of the random variable $\boldsymbol{Y}_{N}$.

Now for comparing different fuzzy models (i.e. different $H_{i}$), the components of the fuzzy inference system (e.g. rule bases, membership functions etc.) would be different, giving rise to a different $g_{i}(\cdot)$ for each case, along with the corresponding set of hyper-parameters $\theta$. For all the demonstrative simulations $\beta_{k}=1,\forall k$.

III Demonstrative examples

III-A Examples with a synthetic dataset

Consider a synthetic example similar to [1] where the downtime of an engineering operation needs to be predicted from two covariates - location risk and maintenance level. For the referential sets $loc\_risk,maintenance,downtime$ the referential values are $\left\{LO,MED,HI\right\}$ , $\left\{POOR,AVG,GOOD\right\}$, $\left\{LO,MED,HI\right\}$ respectively. Also consider another spurious covariate $dummy\_covar$ with referential values $\left\{POOR,AVG,GOOD\right\}$ which is also to be tested for inclusion in model.

The corresponding membership functions are considered to be triangular and can be mathematically expressed as

[TABLE]

where, $u\in U_{i}$ is any value within the universe of discourse and $\phi_{i_{j}}\in\phi_{i}\forall j$ are the parameters of the membership functions to be estimated by the nested sampling algorithm.

Consider the following different rules which would be combined to make up all the different models for evidence calculation.

[TABLE]

The different models or hypotheses ($H$) considered for this example consist of the following sets of individual rules

[TABLE]

The fuzzy model with rule base $H_{true}$ consists of the rules $R_{1}-R_{3}$ and is the true model from which the data is generated. $H_{1}$ consists of two additional rules $R_{4},R_{5}$ which are antagonistic to the original set of rules. In $H_{2}$ one of the original rules $R_{3}$ is replaced with $R_{6}$ which includes the dummy covariate. $H_{3}$ has all the rules as $H_{true}$ with the additional dummy rule $R_{6}$.

The 50 data points for comparing the evidence among the different models is generated by randomly sampling from $H_{true}$ with parameters $\phi_{i}=[5,5,5,5,5,5,50,50,50]$. The universe of discourse for the referential sets $loc\_risk,maintenance,dummy\_covar$ is $[0,10]$ and for $downtime$ is $[0,100]$. The $dummy\_covar$ is also randomly sampled within $[0,10]$ and appended to the input dataset.

We also compare the traditional GLMs with linear predictors to show their evidence on the same data-set. We use the Gaussian likelihood function as given in Eqn. (13) for this. Each element of the mean vector for these GLMs are linear predictors of the form

[TABLE]

where $\alpha_{i}\forall i$ are the parameters of the models. The range of all $\alpha_{i}$ for the GLMs is between -50 to 50. Since the evidence calculation is also dependent on the support of the prior, changing this range would influence the final evidence values. We have shown such simulations later for the real world case.

The evidence calculations from these multiple models for the synthetic dataset is shown in Table I using two different implementations of the nested sampling algorithm - single and multi-ellipsoidal methods. The corresponding mean and standard deviation of the posterior distributions of the model parameters are shown in Table VI.

Both models $H_{2}$ and $H_{3}$ have three additional parameters which estimates the hyper-parameters of the membership function of the dummy variable. The posterior values of these are not shown in the appendix in Table VI as they cannot be compared with respect to a true value. All the simulations in the paper are done using the python nested sampling library nestle [23]. For both the synthetic and the real world case the number of particles or live points are set to 50 and the algorithm terminates when the difference between the estimated total log evidence and the current log evidence is less than the user defined threshold of 0.5.

Table I shows that the model $H_{true}$ has the highest value of evidence $(log(Z))$, with $H_{3}$ following close behind, using the multi-ellipsoidal method which is expected. However with the single ellipsoidal method, $H_{3}$ has a marginally lower $(log(Z))$ value than $H_{true}$ implying that $H_{3}$ is preferred. Now the error estimates in the $log(Z)$ indicates that this difference is within the error margin and indeed looking at Eqn. (17) we see that $H_{3}$ consists of all the rules of $H_{true}$ along with a rule composed of the dummy variable which would have negligible impact on the model output. Either the number of particles in the algorithm needs to be increased to discern this (this increases accuracy with a consequent increase in computational run time) or 50 data points might not be enough and more data is needed. All the other fuzzy models ($H_{1},H_{3}$) and GLMs have a much lower evidence value. Nevertheless, the evidence metric helps in quantitatively ranking the different alternative models and deciding how worse or better they are in comparison to each other. As also shown in Table I, the number of function calls (Fcalls) for the multi-ellipsoidal algorithm is generally lower or the same as that of single ellipsoidal algorithm except in some cases like $H_{GLM5}$. This implies that in cases where the likelihood function is not multi-modal and both the algorithms would give similar answers, the multi-ellipsoidal algorithm would be useful in terms of run time reduction, especially for expensive likelihood functions (e.g. fuzzy rule based models).

III-B A real world example

We use the same dataset and problem setting of modelling insurance risks for power plants as described in [1]. Please refer to [1] for detailed background of the problem and the rationale for applying the Fuzzy Bayesian methodology. The key idea is to encode the underwriters opinions as a fuzzy rule base expressing dependence between the three covariates – operation and maintenance (O&M), loss history, design & technology & fire risk (DnTnF) which are used in predicting the system uptime and consequently helping in making insurance premium calculations and underwriting decisions where the dataset is very small.

Table II shows three fuzzy models with different rule bases that are compared for the real world data-set using the evidence metric. $H_{rw1}$ is the same rule base as used in [1] indicating that the underwriter perceives the O&M to be more indicative of risk than the other covariates. $H_{rw2},H_{rw3}$ indicate that the loss history and the DnTnF are the key covariates respectively.

In addition we also compare three GLM models, as given by Eqn. (19), on the real world data-set. Three sets of simulations are conducted for the real world case, where $\Sigma=\sigma^{2}I_{N}$ in Eqn. (13) is estimated from the data itself and $\sigma$ is set arbitrarily to constant values of 0.25 and 1.0. The corresponding simulation results are reported in Tables III, IV and V respectively. The tables also contain an additional column for the prior ranges. The fuzzy model parameters have been defined to be within the range of $[0,10]$, however the GLM parameters can be any real number and hence two different prior ranges are tested as choice of prior range influences the evidence values. As can be observed from Table III, an uniform prior over a narrower range $[-10,10]$ has a bit higher value of evidence than over a wider range of $[-100,100]$. This might be due to the fact that we are distributing the same probability mass over a narrower range (having more confidence in them as we are ruling out the other values) as opposed to a wider one. However, this is not necessarily the case when $\sigma$ is artificially set to a lower value in Table IV (e.g. $H_{GLM7},H_{GLM8}$).

Table III shows that $H_{rw1}$ has the highest model evidence among the competing alternatives. Arbitrarily setting the value of $\sigma=1.0$ as done in Table V still has $H_{rw1}$ as the best model but the ranking of the other models change. This effect is more predominant in Table IV where $\sigma$ is artificially set to an even lower value of 0.25 and $H_{GLM8}$ has the lowest value of evidence among all the competing alternatives.

IV Discussions and challenges

The example cases showed the model comparison between different rule bases arrived at by different thought processes by multiple experts who do not concur. The method however is much more generic and can also be used to choose other hyper-parameters in the fuzzy model as well. For example, if there is reason to believe that trapezoidal membership functions might be a better representation of the fuzzy sets than the triangular ones used here, then a model with the same rule base but with different shapes of membership functions can be used and the evidence of the model vis-à-vis the other one can decide which type of membership function is best supported by the data.

One advantage of the method was that it presents a unified way for model comparison among widely different modelling philosophies (like the FBL models and the GLMs). This would not have been possible with other approximate methods like AIC, DIC or BIC as the number of ‘effective’ parameters is difficult to find out and have different implications for the different models. Also the use of nested sampling directly produces the evidence along with the posterior parameter estimates making the parameter estimation and model comparison in a single step as opposed to a two stage process.

To lend more credibility to the rankings of the different models based on their evidence, it would be better to perform some sort of sensitivity analysis by checking how the rankings change depending on modelling assumptions like choice of prior, likelihood etc. If the results are not very sensitive to these variations, then the results are more reliable.

V Conclusion

This paper looked at comparing different rule bases in Fuzzy Bayesian Learning using marginal likelihood. Synthetic and real world case studies show how the methodology is able to quantitatively rank the alternative rule bases and also the pitfalls of setting model parameters (like variance) arbitrarily instead of learning from the small dataset. The methodology is generic and can compare other types of models (e.g. GLMs) as demonstrated in the case studies. The marginal likelihood automatically trades-off between model complexity and predictive accuracy and does not need special parameters for tweaking unlike other information theoretic criteria and hence is easy to apply to a large set of models. Future work can look at automatically generating better and interpretable rule bases with higher evidence values using machine intelligence.

Acknowledgment

The authors would like to thank Sciemus Ltd. for sponsoring this work.

Appendix A

GLMs for the real-world dataset.

[TABLE]