{'fraction_non_alphanumeric': 0.049484038380302936, 'fraction_numerical': 0.039376018345302036, 'mean_word_length': 4.652908067542214, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Evaluating the optical classification of Fermi BCUs using machine learning

In the third catalog of active galactic nuclei detected by the Fermi-LAT (3LAC) Clean Sample, there are 402 blazars candidates of uncertain type (BCU). Due to the limitations of astronomical observation or intrinsic properties, it is difficult to classify blazars using optical spectroscopy. The potential classification of BCUs using machine learning algorithms is essential. Based on the 3LAC Clean Sample, we collect 1420 Fermi blazars with 8 parameters of {\gamma}-ray photon spectral index, radio flux, flux density, curve significance, the integral photon flux in 100 to 300 MeV, 0.3 to 1 GeV, 10 to 100 GeV and variability index. Here, we apply 4 different supervised machine learning (SML) algorithms (\emph{Decision trees, Random forests, support vector machines and Mclust Gaussian finite mixture models}) to evaluate the classification of BCUs based on the direct observational properties. All the 4 methods can perform exceedingly well with a more accuracy and can effective forecast the classification of Fermi BCUs. The evaluating results show the results of these methods (SML) are valid and robust, where, about 1/4 sources are FSRQs and 3/4 are BL Lacs in 400 BCUs, which are consistent with some other recent results. Although a number of factors influence the accuracy of SML, the results are stable at a fixed ratio 1:3 between FSRQs and BL Lacs, which suggests that the SML can provides an effective method to evaluate the potential classification of BCUs. Among the 4 methods, Mclust Gaussian Mixture Modelling has the highest accuracy for our training sample (4/5, seed=123).


INTRODUCTION
Blazars are a peculiar sub-class of radio-loud active galactic nuclei (AGNs), whose broadband emission is mainly dominated by non-thermal components produced in a relativistic jet pointed at a small viewing angle to the line of sight (Urry & Padovani 1995). According to the features of optical emission-line in blazars, they are traditionally sub-divided into two groups: flat spectrum radio quasars (FSRQs) and BL Lacertae objects (BL Lacs), where the BL Lacs have weak or no emission lines (e.g., equivalent width, EW, of the emission line in rest frame is less than 5Å) while FSRQs show stronger emission lines (EW ≥ 5Å) (Urry & Padovani 1995;Stocke et al. 1991;Stickel et al. 1991) in their optical spectra. The multi-wavelength spectral energy distributions (SEDs) from radio to γ-ray bands of blazars dominantly comes from the non-thermal emission, where the SED normally exhibits a two-hump structure in the logν − logνF ν space. The lower energy hump (peaked at between millimeter and soft X-ray waveband) is nor-mally attributed to the synchrotron emission produced by the non-thermal electrons in the jet, while the second hump (peaked at the MeV-GeV range) mainly comes from inverse Compton (IC) scattering. The location of the peak for the lower energy bump in the SED, ν S p , is used to classify the sources as low (LSP, e.g., ν S p < 10 14 Hz), intermediate (ISP,e.g.,10 14 Hz < ν S p < 10 15 Hz) and high-synchrotron-peaked (HSP, e.g., ν S p > 10 15 Hz) blazars (Abdo et al. 2010a).
Classified FSRQs and BL Lacs are sources with their optical classifications can be well identified from the literature and/or optical spectrum in the 3FGL catalog Acero et al. 2015). BCUs are the sources with their counterparts have been established. However, their optical classifications have not been identified as a FSRQ or a BL Lac from the weaker or lacking an optical spectrum, and/or their synchrotron peak frequencies of SED, and/or their broadband emission shows blazar-type characteristics with a flat radio spectrum (see Ackermann et al. 2015;Acero et al. 2015 for the details and references therein). Such a large sample of blazars provides a good chance to explore the nature of γ-ray emission of blazars (e.g., Singal et al. 2012;Xiong & Zhang 2014;Singal 2015;Xiong et al. 2015a,b;Chen et al. 2016;Fan et al. 2016a,b;Ghisellini 2016;Lin & Fan 2016;Lin et al. 2017;Chen 2018;Kang et al. 2018;Lin & Fan 2018). In the 3LAC Clean Sample, there are about 30% of blazars (BCUs) that have no optical classification. Evaluating potential classification of the BCUs is a meaningful topics, which have been extensively explored based on the Fermi source catalogs (e.g., see Hassan et al. 2013;Doert & Errando 2014;Chiaro et al. 2016;Einecke 2016;Saz Parkinson et al. 2016;Lefaucheur & Pita 2017;Salvetti et al. 2017;Yi et al. 2017 for the reviews and references therein).
At present, machine learning and data mining techniques are developing rapidly, which has been widely used in the study of astronomy and astrophysics (e.g., see the review in Ball & Brunner 2010;Feigelson &Babu 2012 andWay et al. 2012; also see Ackermann et al. 2012;Mirabal et al. 2012;Hassan et al. 2013;Doert & Errando 2014;Chiaro et al. 2016;Einecke 2016;Saz Parkinson et al. 2016;Lefaucheur & Pita 2017;Salvetti et al. 2017;Yi et al. 2017 ;Bai et al. 2018;Ma et al. 2018). Supervised machine learning (SML) is the most common technique, which aim is to build a classifier (or a decision rule) from the observations of known classification, to classify others (an observation with an unknown class membership to one of K known classes). In the 3LAC Clean Sample, there are about 70% sources with known optical classification (414 FSRQs ∼ 30%, 604 BL Lacs ∼ 40%), however, there are about 30% of blazars (BCUs) that have no optical classification. Evaluating potential classification of the BCUs using supervised machine learning is an interesting work.
In this work, we employ 4 supervised machine learning algorithms (Decision trees (DT), Random forests (RF), support vector machines (SVMs) and M clust Gaussian finite mixture models (Mclust)) to evaluate the potential classification of BCUs only based on (only focus on) the direct observational properties of the 3LAC Clean Sample. We give some description on the sample selection in Section 2, and the supervised machine learning techniques are introducted in Section 3. Section 4 reports the results of supervised machine learning. The discussion and conclusion are presented in Section 5.

SAMPLE
From the 3FGL catalog (Acero et al. 2015) and 3LAC Clean catalog , we select 1420 Fermi Clean blazars (including 414 FSRQs, 604 BL Lacs and 402 BCUs) with 37 variables. In order to select suitable parameters for supervised machine-learning, and to built an available supervised classifier, the independence of these 37 parameters distributions between two subsamples (414 FSRQs and 604 BL Lacs) are calculated using two sample test (KS test, t-test and Wilcox-test) (e.g., Acuner & Ryde 2018). Based on the two sample test results (see Table 1), excluding the same, similar, and related parameters, or some parameters that are directly related to classification (e.g., redshift), 8 parameters (the γ-ray photon spectral index (Γ ph ), radio flux (logF R ), flux density (logF D ), curve significance (C S ), the integral photon flux in 100 to 300 MeV (logF 1 ), 0.3 to 1 GeV (logF 2 ), 10 to 100 GeV (logF 3 ) and variability index (VI)) with the better test results (e.g, D > 0.2 in KS test, or p 3 < 1.00E − 18) are selected in this work. Here, some of these 8 parameters (e.g., spectral Note-Column 1 shows the parameters selected in sample; Column 2 and Column 3 give the value of the test statistic (D) and p-value (p1) for the two-sample Kolmogorov−Smirnov test; The value of the t-statistic (t), the degrees of freedom for the t-statistic (df) and the p-value (p2) for the Welch Two Sample t-test are listed in Column 4, Column 5 and Column 6 respectively; Column 7 and Column 8 report the value of the test statistic (W ) and p-value (p3) for the Wilcoxon rank sum test with continuity correction; Column 9 lists the Gini coefficient (Gini) that is a natural measure of variable importance in Random forests.
index and variability index) are also used in other recent works (e.g., Doert & Errando 2014, Chiaro et al. 2016and Lefaucheur & Pita 2017. However, their research focus (e.g., aims and/or selected parameters and/or methods) are different from that of our work. For instance, Doert & Errando (2014) focused and identified "AGN" or "non-AGN" from 576 unassociated sources of the 2FGL catalogue using a neural network and a random forest SML algorithms; Chiaro et al. (2016)  The fields of unsupervised and supervised machine learning provide many classification methods for predicting categorical outcomes, including logistic regression, decision trees, random forests, support vector machines, neural networks, Bayesian networks, Gaussian finite mixture models, and many others (e.g., see Feigelson & Babu 2012;Kabacoff 2015 for the reviews).
In supervised learning (e.g., see Feigelson & Babu 2012; Kabacoff 2015 for more detail), a dataset containing values for both the predictor variables and the outcome is divided into a training sample and a validation sample. Then one uses the training sample to develop a predictive model, while uses the validation sample to verify the accuracy. This dividing of data is essential for creating an effective model, since one needs a separate validation sample to make a realistic estimation of the effectiveness of the classification schemes developed on a training sample. Once an effective predictive model is created, one can use it to predict outcomes when only the predictor variables are known (e.g., see Feigelson & Babu 2012;Kabacoff 2015 for more detail).
In this section, a brief introduction to DT, RF, SVMs and Mclust is provided. "Decision trees" aim to build a tree that can be used to classify new observations into one of two groups, by creating a set of binary splits on the predictor variables. They are popular in datamining techniques (Utgoff 1989;Duda & Stork 2001), but very often, they tend to produce a large tree and suffer from overfitting (e.g., see Breiman et al. 1984;Duda & Stork 2001).
"Random forests" involve a large number of decision trees from a single training sample. The strategy is to enhance the classification by conducting votes among those many trees. The method is presented by Breiman  . Classification plots for the γ-ray photon spectral index (Γ ph ), radio flux (logFR), flux density (logFD), curve significance (CS), the integral photon flux in 100 to 300 MeV (logF1), 0.3 to 1 GeV (logF2), 10 to 100 GeV (logF3) and variability index (VI) according to the estimated classification from M clustDA for fermi BCUs dataset (predicted dataset), where red empty squares represent FSRQs and blue points represent BL Lacs identified using the function "predict.M clustDA()".
(2001) and is applied to an astronomical dataset by Breiman et al. (2003). It is a highly effective method, producing classification with better accuracy compared with other classification methods (e.g., Decision trees) (e.g., see Fernández-Delgado et al. 2014). Additionally, it can handle problems with many observations and variables, and it can handle the cases where there are large amounts of missing data in the training set and where the number of variables is much greater than the number of observations. Another advantage of RF is that it produces OOB (out-of-bag) error rates and measures of variable importance. On the other hand, due to the large number of trees (default 500 trees), it is difficult to understand the classification rules and make communications. Support vector machines are a group of SML models that can be used for classification and regression (Vapnik 1995(Vapnik , 2000. The mathematical theory behind it is to find the optimal hyperplane or a set of hyperplanes for separating classes in a high-dimensional space. This method produces accurate predicting models and is popular at present. M clust (Scrucca et al. 2016) is a powerful R package for model-based clustering, classification, density estimation to discriminant analysis. It is based on finite Gaussian mixture modelling and provides several tools for finite mixture models, including functions for parameter estimation using the EM algorithm (e.g., Fraley &Raftery 2002 andScrucca et al. 2016).
In order to estimate which approach is most accurate, and therefore to choose the best predictive solution, we define the quantity accuracy in binary classifications context. A function perf ormance() for calculating these statistics (Kabacoff 2015) is provided. The perf ormance() function takes a table containing the true outcome (rows) and predicted outcome (columns) and returns the five accuracy measures. First, the number of true positives (a true positive is an outcome where the model correctly predicts the positive class), true negatives (a true negative is an outcome where the model correctly predicts the negative class), false positives (a false positive is an observation with a positive classification is correctly identified as positive) and false negatives (a false negative is an observation with a negative classification is correctly identified as negative) are extracted. Next, these counts are used to calculate the sensitivity, specificity, positive and negative predictive values, and the accuracy (Kabacoff 2015).
In this work, we will use the rpart packages in R to create decision trees; the randomF orest package to fit random forests; and the e1071 package to build support vector machines. the M clust packages to fit Gaussian Mixture Modelling with the M clustDA() function in the base R 3 installation.

RESULTS
The selected sample (1416 sources) includes 1016 sources (413 FSRQs and 603 BL Lacs) with known optical classification and 400 BCUs with unidentified the optical classification in F ermi 3FGL catalog. In supervised learning, we randomly (random seed =123) assign approximately 4/5 of the observations (known classification: 603 BL Lacs and 413 FSRQs) to the training dataset, and the remaining ones to the validation dataset (test set) in the 8-dimensional parameter space as described in Section 2. All the 400 BCUs are viewed as a sample for prediction (as forecast dataset). The training set has 813 blazars (476 BL Lacs, 337 FSRQs), and the validation set has 203 blazars (127 BL Lacs, 76 FSRQs). The training dataset is used to create classification schemes using a decision tree, a random forest, a support vector machine and a Gaussian Mixture Modelling (M clust). Where, in order to simplify calculating, all the default settings for each of 4 classification function (e.g., M clustDA(), randomF orest(), rpart() and svm() function) are used in this work. The validation dataset is used to evaluate the effectiveness of these schemes. Using an effective predictive model that is developed using the data in the training set to forecast dataset, one can predict outcomes (a BCU belongs to BL Lacs or FSRQs) in situations where only the predictor variables are known. Here, the main R steps can be obtained from a public website 4 .
In M clust discriminat analysis, we use the function M clustDA() with MclustDA model (modelType = "MclustDA", where each known classification is modeled by a finite mixture of Gaussian distributions with a number of components and covariance matrix structures being different between classes, named as MclustDA , see, e.g., Fraley &Raftery 2002 andScrucca et al. 2016) to the training dataset. The largest BIC (Bayesian Information Criterion) value of -6370.124 was obtained using the VEV model (assuming clusters having ellipsoidal distributions described by variable volumes, equal shapes and variable orientations) with a 4-component mixture distribution for 476 BL Lacs; and EVE model (assuming clusters having ellipsoidal distributions described by equal volumes, variable shapes and equal orientations) with a 4-component mixture distribution for 337 FSRQs; based on the training sample. The train- Note-Column 1 shows the different combinations of parameters (8, 4 or 3 parameters, see Section 5), different quality (e.g., seed=123, 321 or 1234) and different quantity (e.g., 4/5 or 2/3) of the training samples. Column 2 gives the classifiers; Column 3 and Column 4 give the number of BL Lacs and FSRQs predicted by a supervised classifier (using supervised machine-learning techniques) to the BCUs (predicted dataset). A classifier's sensitivity, specificity, positive predictive power, and negative predictive power are shown in Column 5, Column 6, Column 7 and Column 8 respectively; The accuracy of a classifier is reported in Column 9.
ing error rate ≃ 0.116 is also obtained based on the the function M clustDA().
In order to test the result of M clust supervised learning (M clustDA() discriminate analysis), the predict.M clustDA() function are used to the test set, so the test error rate ≃ 0.084 is reported. Here, we also compute the classification error using cross-validation. A cross-validation error ≃ 0.138 (which is approximately consistent with the training error rate ≃ 0.116) can be computed using the cvMclustDA() function, which by default use nfold = 10 for a 10-fold cross-validation. The classification for the training dataset from M clustDA are shown in Figure 1. In order to evaluate the utility of a classification scheme, the perf ormance() function is performed and returns Sensitivity = 0.895, Specificity = 0.929, Positive Predictive Value = 0.883, Negative Predictive Value = 0.937, Accuracy = 0.916 (see Table 2 and Figure 2). Using the function "predict.M clustDA()" for classifying predicted dataset (BCUs), we obtain 295 BL Lacs and 105 FSRQs (see Figure 3, Table 2 and machine-readable supplementary material in Table 4) from the 400 BCUs (2 sources have no radio data are excluded).
In Randomf orests discriminat analysis, using the randomF orest() function (the default number of trees is 500) in the random-Forest R package (Liaw & Wiener 2002) to the training dataset, OOB (out-of-bag) estimate of error rate = 0.124 was obtained. Where ran-dom forests also provide a natural measure of variable importance: the Gini coefficient (see Table 1) of Spectral Index = 90.89, flux density = 74.98, Radio flux = 55.55, variability index = 45.65, flux 10−100 GeV = 47.94, flux 100−300 MeV = 32.89, flux 0.3−1 GeV = 24.04, curve significance = 22.20, which suggests that Spectral Index is the most important variable and curve significance is the least important among the 8 selected parameters. Applying the predictive model obtained from the random forest to the validation sample, the validation sample is classified and the predictive Sensitivity = 0.842, Specificity = 0.890, Positive Predictive Value = 0.821, Negative Predictive Value = 0.904 and Accuracy = 0.872, are calculated (see Table 2 and Figure 2). Applying the random forest predictive model to the forecast dataset, we obtain 308 BL Lacs and 92 FS-RQs (see Table 2 and machine-readable supplementary material in Table 4) from the 400 BCUs.
For the training sample, a decision tree is grown using the rpart() function in R package (Therneau & Atkinson 2018). However, unfortunately, the tree sometimes becomes too large and suffers from overfitting (e.g., Breiman et al. 1984;Duda & Stork 2001). To make up for the deficiency, a prune() function is used to prune back the tree in the rpart package. And then a tree with the desired size can be obtained. Using it to the validation sample, the Sensitivity = 0.829, Specificity = 0.858, Positive Predictive Value = 0.778, Negative Predictive Value = 0.893, Accuracy = 0.847, are shown (see Table 2 and Figure 2). Then using it to the forecast dataset, 279 BL Lacs and 121 FSRQs are obtained (see Table 2 and machine-readable supplementary material in Table 4).
Finally, Support vector machines is also applied to the training sample. The svm() function in the e1071 R package (Meyer et al. 2018) is used. Using the optimal predictive model obtained from SVMs to the validation sample, the Sensitivity = 0.829, Specificity = 0.921, Positive Predictive Value = 0.863, Negative Predictive Value = 0.900, Accuracy = 0.887, are printed (see Table 2 and Figure 2). Also using the optimal predictive SV M s model to the forecast dataset, 301 BL Lacs and 99 FSRQs are obtained (see Table 2 and machinereadable supplementary material in Table 4).

DISCUSSIONS AND CONCLUSIONS
In this work, one try to evaluate the potential classification of F ermi BCUs using the supervised machine learning (discriminant analysis). We use 4 methods (DT, RF, SVMs and Mclust ) to perform the discriminant analysis for 8 parameters (Γ ph , logF R , logF D , C S , logF 1 , logF 2 , logF 3 and log VI). All the 4 classifiers per-form exceedingly well and produce accurate and effective forecast of the classification of F ermi BCUs. Compared with the results of these methods, M clust Gaussian Mixture Modelling is the most promising (see Table  2) for our training sample (4/5, seed=123).
FSRQs have stronger emission lines (EW ≥ 5Å), while the BL Lacs have weak (EW < 5Å) or no emission lines (e.g., Urry & Padovani 1995); FSRQs show higher luminosity than that of BL Lacs (e.g., see Fossati et al. 1998;Ghisellini et al. 2011 andGhisellini et al. 2017); Based on the 3LAC catalogue, Ackermann et al. (2015) argued that FSRQs tend to have softer spectra, stronger variability and lower peak frequencies in both synchrotron and IC components than BL Lacs; And many others.
These distinctions suggest different physical origin between in FSRQs and in BL Lacs (e.g., Bhattacharya et al. 2016;Fan et al. 2016b;Yang et al. 2018;Boula et al. 2019). The synchrotron radiation peak frequency of FSRQs is significantly lower than that of BL Lacs (e.g., Fossati et al. 1998;Abdo et al. 2009Abdo et al. , 2010bAckermann et al. 2011Ackermann et al. , 2015Ghisellini et al. 2017), this imply that more electron populations lose their energy through synchrotron cooling in FSRQs. In this scenario, we could expect stronger radio emission in FSRQs. For gamma-ray band, it is commonly believed that the gamma-ray radiation in BL Lacs originate from a pure synchrotron self-Compton (SSC) process, (e.g., Mastichiadis & Kirk 1997;Krawczynski et al. 2004;Zheng & Zhang 2011;Zhang et al. 2014;Zheng et al. 2014;Chen 2017;Zheng et al. 2018), while that in FSRQs come from SSC+EC (external Compton) processes (e.g., Sambruna et al. 1999;Böttcher & Chiang 2002;Chen & Bai 2011;Kang et al. 2014;Kang et al. 2016;Zheng et al. 2017). This indicates, for FSRQs, a complex physical process can be expected in Fermi energy bands. The Fermi energy spectrum in FSRQs could be resulted from the spectrum that is superposed other spectra components (e.g., Zheng & Kang 2013;Kang 2017). The Fermi band of FSRQs locating at intersection of both synchrotron self-Compton component and external Compton component could result to a more complex observational features (e.g., Abdo et al. 2009Abdo et al. , 2010bAckermann et al. 2011Ackermann et al. , 2015. Other physical origins (e.g., mass accretion rate on to the central black hole) are also discussed in resent works (e.g., Boula et al. 2019). The more fundamental physical origins between in FSRQs and in BL Lacs require further discussion in the future.
One checks the results for different combinations of parameters. Based on the Gini coefficient (a natural measure of variable importance) in random forests su- Note-Column 1 shows the different classification (− represents the number of mismatch by cross comparison, "bll" ,"fsrq" and "unc" indicate BL Lac, FSRQ and uncertain type respectively, where a and b indicate the number of the disagreement (misjudged as FSRQs or BL Lacs)), and the mismatch rate (e.g., rate = 4/(4+47)% ∼ 7.8% ). pervised learning (see Column 9 in Table 1), one selects part parameters with a higher Gini coefficient (4 parameters: Γ ph , logF R , logVI and logF D ; or 3 parameters: Γ ph , logF R , logF D ) to discriminant analysis also using the 4 methods. We find that the predictive accuracy will be smaller than that of 8 parameters. However, M clust Gaussian Mixture Modelling also tends to be more accurate compared with other classification methods for the different testing variables in combination (see Table  2, e.g., 8 parameters, 4 parameters and 3 parameters) for our training sample. In general, which implies more parameters and more accuracy, but it is unstable for the different classification methods (see Table 2 for the details). However, we should note that the predictive accuracy and results may be affected by the training dataset and validation dataset. When one randomly (seed=123) assigns approximately 2/3 of the known classification blazars (603 BL Lacs and 413 FSRQs) to the training dataset (677 blazars: 399 BL Lacs and 278 FS-RQs) and the remaining ones to the validation dataset (339 blazars: 204 BL Lacs and 135 FSRQs) in the 8-dimensional parameter space as in Section 4. The predictive accuracy (Accuracy = 0.897) and results (304 BL Lacs and 96 FSRQs predicted from the 400 BCUs) are slightly different with that of the (4/5, seed=123) training and validation samples (see Table 2) in M clust discriminant analysis. And other methods also show similar results. Also, for the randomly samples (e.g., randomly seed =123, =321, or =1234), the results and accuracy are also different (see Table 2), where the most accuracy are obtained in the support vector machines (seed=321) or in the Random forests (seed=1234) respectively. These suggest the results of discriminant analysis (supervised learning) are significantly affected by the quality (e.g., seed=123, 321 or 1234) and quantity (e.g., changed from 813 blazars (4/5, see Section 4) to 677 blazars (2/3)) of the training samples. Where, sometimes, the support vector machines or the Random forests yield a higher accuracy, which is consistent with other works (e.g., Fernández-Delgado et al. 2014).
In addition, we also should note that all the default settings for each of the 4 classification function (e.g., M clustDA(), randomF orest(), rpart() and svm() func-tion) are used in Section 4. For each different classification method, choosing of calculation model and setting of each parameter in fitting function (e.g., the "model-Type = MclustDA or EDDA (Eigenvalue Decomposition Discriminant Analysis, e.g., see Scrucca et al. 2016))" in M clustDA() function , the "tree = 500 or 10000" in randomF orest() function, if a prune() function is used in rpart(), "gamma=0.1 or 0.01" and "cost=1 or 1000" in svm() function) can also affect predictive models, accuracy and results (see the labels "EDDA", "Forest", "no pruned" and "cost" in Table 2). About how to select the appropriate parameter settings are need to further address in future, which also is beyond the scope of this work.
We compare the results of M clust SML algorithms with that of other three SML algorithms. We find that, for BL Lacs, 253, 277 and 274 (mean about 91%) BL Lac candidates in DT, RF and SVMs match with M clust (295) BL Lacs candidate sample (see Table 3); unfortunately, 42, 18 and 21 (mean about 9%) sources classed as FSRQ do not match the M clust (295) BL Lacs candidates respectively. For FSRQs, 79, 74 and 78 (mean about 73%) FSRQ candidates match the results (105 FSRQ candidates) of M clust method; but 26, 31 and 27 (mean about 27%) sources do not match the subset of M clust 105 FSRQ candidates respectively (see Table  3).
We also compare the results of M clust algorithms with other resent similar results (e.g., Chiaro et al. 2016;Massaro et al. 2016;Lefaucheur & Pita 2017;Yi et al. 2017). After cross comparison with the results of Lefaucheur & Pita (2017) using multivariate classifications, in the subset of 295 BL Lac candidates (see M clust method), we find that 3 sources do not match sources and 28 sources did not provide a clear classification in Lefaucheur & Pita (2017). M clust prediction is in accordance with Lefaucheur & Pita (2017) for 228 objects (about 77%) and is inconsistent for 36 (about 12%). For the subset of 105 FSRQ candidates, 83 objects are in agreement with and 14 objects are in disagreement with M clust prediction, and 8 objects do not provide a clear classification (Lefaucheur & Pita 2017). Also, in the subset of 295 BL Lac candidates, M clust prediction is in accordance with Chiaro et al. (2016) using artificial neural networks (ANN) machine-learning techniques for 244 sources and is inconsistent for 22, 29 objects do not provide a clear classification (Chiaro et al. 2016). For the subset of 105 FSRQ candidates, only 54 objects are in agreement with and 29 objects are in disagreement with M clust prediction, and 22 objects do not provide a clear classification (Chiaro et al. 2016). For comparison with Yi et al. (2017) performed a statistical analysis of the broadband spectral properties (e.g., spectral indices in the gamma-ray, X-ray, optical, and radio bands) of blazars, the similar results are also shown. 218 BL Lac candidates and 97 FSRQ candidates are in agreement with M clust prediction (295 BL Lac and 105 FSRQ candidates); 4 BL Lac or 8 FSRQ candidates objects are in disagreement with M clust prediction (see Table  3). Most of the results of M clust SML are consistent with that (47 BL Lacs ∼92.2% and 7 FSRQs ∼77.8%) of Massaro et al. (2016) using optical spectroscopic observations. However, a fraction of sources (4 BL Lacs and 2 FSRQs) are misjudged using M clust SML. These results suggest that it is a good overall agreement in these SMLs and other resent results. However, SML algorithms probably lead to some misjudgments in evaluating the potential (optical) classification of blazars. Only the optical spectroscopic observations is still most efficient and accurate way to determine the real nature of these sources.
Although the discriminant analysis can return the probabilities P Bi and P F i that a BCU i belongs to the BL Lacs (B) or FSRQs (F) classifications, respectively (e.g., in M clust method, see Table 4 and a machinereadable supplementary material). However, it should be noted that the error of the supervised machine learning is still very large (accuracy is still not high enough) in the work, where, the accuracy is less than 92%. It also probably leads to some misjudgments that some FSRQs are falsely classified as BL Lacs, and vice versa (see discussion above). These result may be biased, or only be an apparent phenomenon, or hide the essential difference between BL Lacs and FSRQs (e.g., see Blandford et al. 2018 for the reviews and Boula et al. 2019), which needs further consideration. Although we did not conclusively evaluate their potential classifications (FSRQs and BL Lacs) using SML, it may be helpful for source selections in the spectroscopic observation campaigns in the future performing a spectroscopic and photometric campaign, further diagnosing their optical classification of BCUs (e.g., see Yi et al. 2017;Massaro et al. 2013 for some discussion), or provide some clues for future studying of spectroscopic and photometric.
Finally, it must be highlighted that, in this work, our results are obtained only by supervised machine learning the data obtained from Fermi catalogue, adding no external data (obtained from other archives), and have not done any of the fittings ourselves. The limited sample (made of mostly bright, firmly classified sources, excluded fainter sources) is used to diagnose the optical classifications of the BCUs in 3LAC Clean Sample. Selection effects of the direct observational data because of detection thresholds and energy bands in instruments may affect the source distributions and affect the results of the analysis in this work.
However, it should point out that each of these classifiers (DT, RF, SVMs and Mclust ) performed exceedingly well on each of the accuracy measures. The results of SML are in agreement with each classifiers and other resent results (Chiaro et al. 2016;Massaro et al. 2016;Lefaucheur & Pita 2017;Yi et al. 2017), which suggests that the SML can provide an effective and easy method to evaluate the potential classification of BCUs. The evaluating results show the approach (SML) is valid and robust (see Table 2). It is about 1:3 ratio between FS-RQs and BL Lacs predicted from the 400 BCUs for any SML algorithms. Here, we also should note that 1:4 ratio (64 FSRQs and 246 BL Lacs) was obtained by combinating the results of these 4 methods (see Table  3). Whether the true ratio is 1:3 or 1:4 or others needs further verification. However, M clust Gaussian Mixture Modelling tend to be more accurate compared with other classification methods for the different testing the variables in combination (see Table 2, e.g., 8 parameters, 4 parameters and 3 parameters) for our training sample. Although there are a number of factors influencing the accuracy of SML. However, this work provides some sim-ple methods to distinguish the BL Lacs or FSRQs with the probabilities P Bi and P F i (see Table 4) from BCUs based on the direct observational data. A more preferable statistical approach, that uses the a large and more complete sample (e.g., the upcoming 4LAC) are needed to further test and address the issue.