<figure><img src="image_1.png" /></figure>

HAL
open science

# A Specialized Probability Density Function for the Input of Mixture of Gaussian Processes

Longbo Zhao, Jinwen Ma

► To cite this version:

Longbo Zhao, Jinwen Ma. A Specialized Probability Density Function for the Input of Mixture of Gaussian Processes. 2nd International Conference on Intelligence Science (ICIS), Nov 2018, Beijing, China. pp.70-80, 10.1007/978-3-030-01313-4_8 . hal-02118817

HAL Id: hal-02118817

https://inria.hal.science/hal-02118817

Submitted on 3 May 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Specialized Probability Density Function for the Input of Mixture of Gaussian Processes

Longbo Zhao and Jinwen Ma*

Department of Information Science
School of Mathematical Sciences & LMAM
Peking University, Beijing, 100871, P.R. China

**Abstract.** Mixture of Gaussian Processes (MGP) is a generative model being powerful and widely used in the fields of machine learning and data mining. However, when we learn this generative model on a given dataset, we should set the probability density function (pdf) of the input in advance. In general, it can be set as a Gaussian distribution. But, for some actual data like time series, this setting or assumption is not reasonable and effective. In this paper, we propose a specialized pdf for the input of MGP model which is a piecewise-defined continuous function with three parts such that the middle part takes the form of a uniform distribution, while the two side parts take the form of Gaussian distribution. This specialized pdf is more consistent with the uniform distribution of the input than the Gaussian pdf. The two tails of the pdf with the form of a Gaussian distribution ensure the effectiveness of the iteration of the hard-cut EM algorithm for MGPs. It demonstrated by the experiments on the simulation and stock datasets that the MGP model with these specialized pdfs can lead to a better result on time series prediction in comparison with the general MGP models as well as the other classical regression methods.

**Keywords:** Gaussian distribution, Mixture of Gaussian processes, Hard-cut EM algorithm, Probability density function, Time series prediction.

## 1 Introduction

Gaussian process (GP) is a powerful model and widely used in machine learning and data mining [1–3]. However, there are two main limitations. Firstly, it cannot fit the multi-modal dataset well because GP model employs a global scale parameter [4]. Secondly, its parameter learning consumes $O(N^3)$ computational time [5, 6], where N is the number of training samples. In order to overcome those difficulties, Tresp [4] proposed mixture of Gaussian processes (MGP) in 2000, which was developed from the mixture of experts. Since then, many kinds of MGP model have been proposed and can be classified into two main forms: the generative model [7–10] and conditional model [4, 6, 11–13]. In comparison with the conditional model, the generative model has two main advantages: (1).

<footnote>* Corresponding author. Email:jwma@math.pku.edu.cn</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** The sketch of the eLoad data. **Fig. 2.** The sketch of the transformed eLoad data.</figcaption></figure>

The missing features can be easily inferred from the outputs; (2). The influence of the inputs on the outputs is more clear [8]. Therefore, many scholars have studied the generative model [14-20].

However, when we learn the generative model on a given dataset, we should set the probability density function (pdf) of the input in advance. In general, it can be set as a Gaussian distribution [14-20]. But, for some actual data like time series, this setting or assumption is not so reasonable and effective. When we learn MGP model on these actual data, we usually need to utilize the ARMA model [14-21] to transform the data, and then use the transformed data on the MGP model. However, this transformation can destroy the correlation of samples, which is very important for MGP model. Fig. 1 shows the eLoad data [14] from which we can see that samples in three different colors (blue, black, and red) represent three temporally sequential samples, respectively. Fig. 2 shows the transformed eLoad data from which we can find that three temporally sequential samples are mixed together and cannot be classified effectively. In this paper, we propose a specialized pdf for the input of the MGP model to solve this problem. As shown in Fig. 3, this pdf consists of three components. The left and right side parts are Gaussian distributions, while the middle is a uniform distribution. For the training of the MGP model, we use the hard-cut EM algorithm [17] as the basic learning framework for parameter estimation. Actually, the hard-cut EM algorithm can get better result than some popular learning algorithms.

The rest of the paper is organized as follows. Section 2 introduces the GP and MGP models. We describe the specialized probability density function in Section 3. We further propose the learning algorithm for the MGP model of the specialized pdfs in Section 4. The experimental results are contained in Section 5. Finally, we make a brief conclusion in Section 6.
<figure><img src="image_4.png" /><figcaption>Fig. 3. The sketch of the specialized input distribution.</figcaption></figure>

# 2 GP and MGP Models

## 2.1 GP Model

We mathematically define the GP model as follows:

$$Y \sim N(m(\mathbf{X}), K(\mathbf{X}, \mathbf{X})) \qquad (1)$$

where $D=\{\mathbf{X}, \mathbf{Y}\}=\{(\mathbf{x}_i, y_i):i=1,2,...,N\}$, $\mathbf{x}_i$ denotes a d-dimensional input vector, and $y_i$ is the corresponding output. $m(\mathbf{X})$ and $K(\mathbf{X}, \mathbf{X})$ denote the mean vector and covariance matrix, respectively. Without loss of generality, we assume $m(\mathbf{X}) = 0$. There are many choices for covariance function, such as linear, Gaussian noise, squared exponential function and so on. Here, we adopt the squared exponential (SE) covariance function [10]:

$$K(\mathbf{x}_i, \mathbf{x}_j; \theta) = \sigma_f^2 \exp(-\frac{\sigma_l^2}{2} \|\mathbf{x}_i - \mathbf{x}_j\|^2) + \sigma_n^2 I_{(i=j)} \qquad (2)$$

where $\theta=\{\sigma_f^2, \sigma_l^2, \sigma_n^2\}$ denote the vector. On the given sample dataset $D$, the log-likelihood function can be expressed as follows:

$$\log p(\mathbf{Y} | \mathbf{X}, \theta) = \log N(\mathbf{Y} | \mathbf{0}, K(\mathbf{X}, \mathbf{X})) \qquad (3)$$

In order to obtain the estimation of parameters $\theta$, we perform the maximum likelihood estimation (MLE) procedure [10], that is, we get

$$\hat{\theta} = \underset{\theta}{\operatorname{argmax}} \log N(\mathbf{Y} | \mathbf{0}, K(\mathbf{X}, \mathbf{X})) \qquad (4)$$

## 2.2 MGP Model

Denote $C$ and $N$ as the number of GP components and training samples in the MGP model, respectively. On the basis of the GP model, we define MGP model by the following steps:
Step 1. Partition samples into each GP components by the Multinomial distribution:

$$p(z_n = c) = \pi_c \quad (5)$$

where $c=1,...,C$ and $n=1,...,N$.

Step 2. Accordingly, each input $x_i$ fulfills the following distribution:

$$p(x_i|z_n = c) \sim p(x|\psi_c) \quad (6)$$

where $\{\psi_c: c=1,...,C\}$ is the parameter set. In general, $p(x|\psi_c)$ is a Gaussian distribution.

Step 3. Denote $I_c=\{n|z_n = c\}$, $X_c=\{x_n|z_n = c\}$, $Y_c=\{y_n|z_n = c\}$ ($c=1,...,C$, $n=1,...,N$) as the sample indexes, inputs and outputs of the training samples in the $c$-th component, respectively. Given $X_c$, the corresponding $c$-th GP component can be mathematically defined as follows:

$$Y_c \sim N(0, K(X_c, X_c)) \quad (7)$$

where $K(X_c, X_c)$ is given by Eq.(2) with the hyper-parameter $\theta_c=\{\sigma_{fc}^2, \sigma_{lc}^2, \sigma_{nc}^2\}$.

Based on Eqs.(5),(6),(7), we mathematically define the MGP model. The log-likelihood function is derived as follows:

$$\log(p(\mathbf{Y}_c|\mathbf{X}_c, \mathbf{\Theta}, \mathbf{\Psi})) = \sum_{c=1}^{C} \left( \sum_{n \in I_c} \left( \log(\pi_c p(x_n|\mu_c, S_c)) \right) \right) + \log(p(\mathbf{Y}_c|\mathbf{X}_c, \mathbf{\theta}_c)) \quad (8)$$

where $\mathbf{\Theta}=\{\theta_c: c=1,...,C\}$ and $\mathbf{\Psi}=\{\psi_c, \pi_c: c=1,...,C\}$ denote the hyper-parameters and parameters of the MGP model, respectively.

# 3 Specialized Input Distribution and Its Learning Algorithm

For many real world datasets, such as UCI machine learning repository, Gaussian distribution is not appropriate for the input. In order to solve this problem, we propose a specialized distribution for this situation.

## 3.1 Specialized PDF

This specialized distribution is a piecewise-defined continuous function, which consists of three parts, the middle part is a uniform distribution density, both sides are Gaussian distribution densities, shown in Fig. 3. We mathematically defined the specialized distribution as follows:

$$P(x; \psi) = \begin{cases} \frac{\lambda_1}{(\sqrt{2\pi}\tau_1)} \exp^{-\frac{(x-a)^2}{2\tau_1^2}} & x < a \\ \lambda & a \le x \le b \\ \frac{\lambda_2}{(\sqrt{2\pi}\tau_2)} \exp^{-\frac{(x-b)^2}{2\tau_2^2}} & x > b \end{cases} \quad (9)$$

where we redefine $\psi=\{\lambda, \lambda_1, \lambda_2, \tau_1, \tau_2, a, b\}$ as the parameter vector.
## 3.2 Learning Algorithm for the Specialized PDF

In order to learn $\psi$, we set that the input interval $(a,b)$ contains the number of the samples with probability $p_0$. Denote $\mathbf{X}$ and $N$ as the training sample set and the number of training sample, respectively. We summarize the algorithm framework as following steps:

Step 1. Learn $a$, $b$, and $\lambda$

$$ \mathbf{a} = X_{\frac{N(1-p_0)}{2}}; \mathbf{b} = X_{\frac{N(1+p_0)}{2}}; \lambda = \frac{p_0}{(\mathbf{b} - \mathbf{a})} \quad (10) $$

where $p(x<X_{\frac{N(1-p_0)}{2}}|x \in \mathbf{X})=\frac{(1-p_0)}{2}$. In order to reduce the effect of the misclassified (or outlier) point on the middle part, we estimate $\mathbf{a}$ and $\mathbf{b}$ as Eq. (10) do.

Step 2. Estimate $\lambda_1$, $\lambda_2$, $\tau_1$ and $\tau_2$.

Denote $p_1$ and $p_2$ as the sample ratio at both left side and right side, respectively. The probability density function is continuously integrable, and the integral of the probability density function is equal to 1. In other word:

$$ \int P(\mathbf{x}; \psi)dx = \begin{cases} p_1 & \mathbf{x} < \mathbf{a} \\ p_0 & \mathbf{a} \le \mathbf{x} \le \mathbf{b} \\ p_2 & \mathbf{x} > \mathbf{b} \end{cases}; \quad p_0 + p_1 + p_2 = 1 \quad (11) $$

According to the continuity of the probability density function, we only need do same simple calculations to get $\{\lambda_1, \lambda_2, \tau_1, \tau_2\}$:

$$ \lambda_1 = 2p_1; \lambda_2 = 2p_2; \tau_1 = \frac{\lambda_1}{\sqrt{2\pi}\lambda}; \tau_2 = \frac{\lambda_2}{\sqrt{2\pi}\lambda} \quad (12) $$

# 4 The MGP Model of the Specialized PDFs and Its Learning Algorithm

We now consider the MGP model with these specialized pdfs. For the parameter learning of the MGP model, there are main three kinds of learning algorithms: MCMC methods [22, 23], variational Bayesian inference [24, 25], and EM algorithm [5, 9, 11]. However, the MCMC methods and variational Bayesian inference methods have their own limitations: the time complexity of the MCMC method is very high, and variational Bayesian inference may lead to a rather deviation from the true objective function. EM algorithm is an important and effective iterative algorithm to do maximum likelihood or maximum a posterior(MAP) estimates of parameters for mixture model. However, for such a complex MGP model, the posteriors of latent variables and Q function are rather complicated. In order to overcome this difficulty, we implement the hard-cut EM algorithm [17] to learn parameter, which makes certain approximations in E-step.

Denote $z_{nc}$ be the latent variables, where $z_{nc}$ is a Kronecker delta function, $z_{nc}=1$, if the sample $(\mathbf{x}_n, y_n)$ belongs to the c-th GP component. Therefore,
we can obtain the log likelihood function of the complete data from Eq.(8) as follows:

$$ \log(p(\mathbf{Y}, \mathbf{Z} | \mathbf{X}, \mathbf{\Theta}, \mathbf{\Psi})) = \sum_{c=1}^{C} \left( \sum_{n=1}^{N} (z_{nc} \log(\pi_c p(\mathbf{x}_n | \mathbf{\psi}_c))) \right. \\ \left. + \log(p(\mathbf{Y}_c | \mathbf{X}_c, \mathbf{\theta}_c)) \right) \quad (13) $$

The main idea of hard-cut EM algorithm can be expressed as the following steps:

E-step. Assign the samples to the corresponding GP component according to the maximum a posterior (MAP) criterion:

$$ \hat{k}_n = \operatorname*{argmax}_{1 \le c \le C} \{ \pi_c p(\mathbf{x}_n | \mathbf{\psi}_c) p(y_n | \mathbf{\theta}_c) \} \quad (14) $$

that is, latent variable $\hat{z}_{\hat{k}_n n}=1$.

M-step. With the known partition, we can estimate the parameters $\mathbf{\Psi}$ and hyper-parameters $\mathbf{\Theta}$ via the MLE procedure:

1) For learning the parameters $\{\psi_c\}_c$, we perform the learning algorithm in the last section.

2) For estimating the hyper-parameter $\mathbf{\Theta}$, we perform the MLE procedure on each c-th component to estimate $\mathbf{\theta}_c$ as shown in Eq.(4).

# 5 Experimental Results

In order to test the accuracy and effectiveness of the specialized pdf for MGP model, we carry out several experiments on the simulation and stock datasets. we employ the root mean squared error(RMSE) to measure the prediction accuracy, which is defined as follows:

$$ RMSE = \sqrt{\frac{\sum_{n=1}^{N} (\mathbf{y}_n - \hat{\mathbf{y}}_n)^2}{N}} \quad (15) $$

where $\hat{\mathbf{y}}_n$ and $\mathbf{y}_n$ denote the predicted value and true value, respectively. We also compare our algorithm with some classical machine learning algorithms: kernel, RBF, SVM and denote 'OURS' as our proposed model with the hard-cut EM algorithm.

## 5.1 Simulation Experiments

In the simulation experiments, we generate three groups of synthetic datasets from MGP model. those three MGP models contain 4, 6, 8 GP components, respectively. The number of samples in each group is 2600, 3900, 5000, respectively. In each group, there are three datasets, which are the same except the degree of overlap. Fig. 4 shows the dataset with the smallest degree of overlap with 4 GP components. On each group dataset, we run each algorithm 100 times,
<figure><img src="image_5.png" /><figcaption>Fig. 4. The dataset with the least degree of Fig. 5. The distribution of the probability overlap from MGP with 4 components. density function of the input on each Gaussian processes component.</figcaption></figure>

and randomly extract training samples and test samples, where 1/3 are training samples and other 2/3 are test samples. The RMSE of each algorithm is listed in Table 1. From the Table 1, We can see that our proposed algorithm obtains the better results. Fig. 5 shows the specialized pdfs on the first group dataset with the smallest overlapping degree. We can obtain that the specialized pdf at both ends of the data in the form of a Gaussian distribution of attenuation, the specialized pdf in the middle of the data is a uniform distribution. This shape of the specialized pdf is more consistent with the uniform distribution than Gaussian distribution. The attenuation of both ends of the specialized pdf in the form of a Gaussian distribution ensures the effectiveness of the iteration of the hard-cut EM algorithm. Then, the class label of the samples can be updated according to the MAP criteria in the iteration of hard-cut EM algorithm. If we apply uniform distribution only, the iterative steps of hard-cut EM algorithm is invalid.

## 5.2 Prediction on Stock Data

In this section, we obtain the closing price data of three stocks from Shanghai Stock Exchange, and the IDs are 300015, 002643, and 601058, respectively.

From Eq.(10), we can know that the specialized pdf is closely related to the interval length of the middle data. In order to check the effect of different input lengths on the prediction accuracy of the algorithm, we do some transformations on the input. Since the range of output changes is too large, we use a linear function to narrow the output down to the same range as the synthetic data. In summary, we transform the datasets as follows:

(i) Transform the input as following equation:

$$X_n = \frac{n}{\delta} \qquad (16)$$

where $i=1,...,N$, $N$ is the sample number, $\delta = \{101, 51, 23, 11, 7, 3, 1\}$.
<table><caption>**Table 1.** The RMSEs of the four algorithms on the three groups.</caption><thead><tr><th>C=4</th><th>Data41</th><th>Data42</th><th>Data43</th></tr></thead><tbody><tr><td>Kernel</td><td>0.21430.0000</td><td>0.68710.0000</td><td>0.65370.0000</td></tr><tr><td>RBF</td><td>0.19650.0000</td><td>0.70650.0000</td><td>0.65940.0000</td></tr><tr><td>SVM</td><td>0.25480.0002</td><td>0.37390.0060</td><td>0.37040.0038</td></tr><tr><td>OURS</td><td>0.06040.0003</td><td>0.29910.0319</td><td>0.33260.0174</td></tr></tbody><thead><tr><th>C=6</th><th>Data61</th><th>Data62</th><th>Data63</th></tr></thead><tbody><tr><td>Kernel</td><td>0.52120.0000</td><td>0.52110.0000</td><td>0.53920.0000</td></tr><tr><td>RBF</td><td>0.54110.0000</td><td>0.54980.0000</td><td>0.55540.0000</td></tr><tr><td>SVM</td><td>0.26780.0013</td><td>0.33780.0033</td><td>0.39420.0080</td></tr><tr><td>OURS</td><td>0.25580.0366</td><td>0.36530.0384</td><td>0.33920.0210</td></tr></tbody><thead><tr><th>C=8</th><th>Data81</th><th>Data82</th><th>Data83</th></tr></thead><tbody><tr><td>Kernel</td><td>0.42470.0000</td><td>0.42450.0000</td><td>0.49930.0000</td></tr><tr><td>RBF</td><td>0.45300.0000</td><td>0.42650.0000</td><td>0.49340.0000</td></tr><tr><td>SVM</td><td>0.36420.003</td><td>0.42380.0023</td><td>0.48310.0017</td></tr><tr><td>OURS</td><td>0.32200.0335</td><td>0.36390.0898</td><td>0.46680.3434</td></tr></tbody></table>

(ii) Transform the output by a linearly compressed, and the compressed interval is [-4.5, 4.5].

$$ \tilde{y} = \frac{9y}{M-m} + \frac{4.5}{M-m} \qquad (17) $$

where *M* and *m* denote the maximum value and minimum value of the stock, respectively.

Through the above transformations, each stock can produce 7 datasets. In each 7 datasets of three stock datasets, we repeat each regression algorithm 100 times, and randomly extracted 1/3 as training samples and the other 2/3 as test samples. The RMSE of each algorithm on those three transformed stock datasets is listed in Table 2. From Table 2, we can obtain that our proposed algorithm can get a better predict accuracy than other classical regression algorithms, and our algorithm obtain the better result with the smaller *δ*, but this is not absolute.

# 6 Conclusion

We have designed a specialized pdf for the input of MGP model which consists of three parts: the right and left side parts still take the form of Gaussian distributions, while the middle part takes the form of a uniform distribution. This specialized pdf has the advantages of both the Gaussian distribution and the uniform distribution. That is, the tail Gaussian distributions in the left and right side parts ensure that the hard-cut EM algorithm can perform more efficiently during each iteration, and the uniform distribution in the middle part is more reasonable for the time series data. The experiments are conducted on three groups of synthetic datasets and stock datasets. It is demonstrated by the experimental results that the hard-cut EM algorithm for the MGPs with the
<table><caption>**Table 2.** The RMSEs of the four algorithms on the three groups of the transformed stock datasets.</caption><tbody><tr><td>300015</td><td>X1</td><td>X2</td><td>X3</td><td>X4</td><td>X5</td><td>X6</td><td>X7</td></tr><tr><td colspan="8">Kernel 0.6347 0.5439 0.4121 0.2753 0.2241 0.1846 0.2149</td></tr><tr><td colspan="8">RBF 0.5218 0.3717 0.2918 0.2671 0.3435 0.6930 0.9934</td></tr><tr><td colspan="8">SVM 0.3634 0.2532 0.1971 0.1707 0.1869 0.1707 0.1842</td></tr><tr><td colspan="8">OURS 0.3531 0.2325 0.1784 0.1676 0.1601 0.1467 0.1583</td></tr><tr><td>002643</td><td>X1</td><td>X2</td><td>X3</td><td>X4</td><td>X5</td><td>X6</td><td>X7</td></tr><tr><td colspan="8">Kernel 0.6487 0.4722 0.3116 0.2323 0.1972 0.1533 0.1542</td></tr><tr><td colspan="8">RBF 0.5357 0.4462 0.2744 0.2280 0.2116 0.2110 0.8882</td></tr><tr><td colspan="8">SVM 0.2815 0.2487 0.2364 0.2175 0.1998 0.1902 0.1792</td></tr><tr><td colspan="8">OURS 0.2681 0.2024 0.1835 0.1799 0.1930 0.1805 0.1544</td></tr><tr><td>601058</td><td>X1</td><td>X2</td><td>X3</td><td>X4</td><td>X5</td><td>X6</td><td>X7</td></tr><tr><td colspan="8">Kernel 0.7200 0.5267 0.3939 0.2738 0.2150 0.1503 0.1256</td></tr><tr><td colspan="8">RBF 0.5296 0.4547 0.3636 0.2644 0.2279 0.2785 0.9705</td></tr><tr><td colspan="8">SVM 0.4079 0.2325 0.1748 0.1502 0.1439 0.1439 0.1458</td></tr><tr><td colspan="8">OURS 0.3271 0.1966 0.1318 0.1507 0.1545 0.1315 0.1472</td></tr></tbody></table>

specialized pdfs can obtain a better prediction accuracy than the other classical regression algorithms. This specialized input pdf is more effective for the time series data.

## Acknowledgment

This work was supported by the National Science Foundation of China under Grant 61171138.

## References

1. Rasmussen C E. Evaluation of Gaussian processes and other methods for non-linear regression[M]. University of Toronto, 1999.

2. Williams C K I, Barber D. Bayesian classification with Gaussian processes[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1998, 20(12): 1342-1351.

3. Rasmussen C E, Kuss M. Gaussian Processes in Reinforcement Learning[C]//NIPS. 2003, 4: 1.

4. Tresp V. Mixtures of Gaussian processes[C]//Advances in neural information processing systems. 2001: 654-660.

5. Yuan C, Neubauer C. Variational mixture of Gaussian process experts[C]//Advances in Neural Information Processing Systems. 2009: 1897-1904.

6. Stachniss C, Plagemann C, Lilienthal A J, et al. Gas Distribution Modeling using Sparse Gaussian Process Mixture Models[C]//Robotics: Science and Systems. 2008, 3.

7. Yang Y, Ma J. An efficient EM approach to parameter learning of the mixture of Gaussian processes[C]//International Symposium on Neural Networks. Springer, Berlin, Heidelberg, 2011: 165-174.
8. Meeds E, Osindero S. An alternative infinite mixture of Gaussian process experts[C]//Advances in Neural Information Processing Systems. 2006: 883-890.

9. Sun S, Xu X. Variational inference for infinite mixtures of Gaussian processes with applications to traffic flow prediction[J]. IEEE Transactions on Intelligent Transportation Systems, 2011, 12(2): 466-475.

10. Williams C K I, Rasmussen C E. Gaussian processes for machine learning[J]. the MIT Press, 2006, 2(3): 4.

11. Nguyen T, Bonilla E. Fast allocation of Gaussian process experts[C]//International Conference on Machine Learning. 2014: 145-153.

12. Lzaro-Gredilla M, Van Vaerenbergh S, Lawrence N D. Overlapping mixtures of Gaussian processes for the data association problem[J]. Pattern Recognition, 2012, 45(4): 1386-1395.

13. Ross J, Dy J. Nonparametric mixture of Gaussian processes with constraints[C]//International Conference on Machine Learning. 2013: 1346-1354.

14. Wu D, Ma J. A Two-Layer Mixture Model of Gaussian Process Functional Regressions and Its MCMC EM Algorithm[J]. IEEE Transactions on Neural Networks and Learning Systems, 2018.

15. Wu D, Chen Z, Ma J. An mcmc based em algorithm for mixtures of gaussian processes[C]//International Symposium on Neural Networks. Springer, Cham, 2015: 327-334.

16. Wu D, Ma J. A daem algorithm for mixtures of gaussian process functional regressions[C]//International Conference on Intelligent Computing. Springer, Cham, 2016: 294-303.

17. Chen Z, Ma J, Zhou Y. A precise hard-cut EM algorithm for mixtures of Gaussian processes[C]//International Conference on Intelligent Computing. Springer, Cham, 2014: 68-75.

18. Chen Z, Ma J. The hard-cut em algorithm for mixture of sparse gaussian processes[C]//International Conference on Intelligent Computing. Springer, Cham, 2015: 13-24.

19. Zhao L, Chen Z, Ma J. An Effective Model Selection Criterion for Mixtures of Gaussian Processes[C]// International Symposium on Neural Networks. Springer International Publishing, 2015:345-354.

20. Zhao L, Ma J. A dynamic model selection algorithm for mixtures of Gaussian processes[C]//Signal Processing (ICSP), 2016 IEEE 13th International Conference on. IEEE, 2016: 1095-1099.

21. Liu S, Ma J. Stock price prediction through the mixture of gaussian processes via the precise Hard-cut EM algorithm[C]//International Conference on Intelligent Computing. Springer, Cham, 2016: 282-293.

22. Shi J Q, Murray Smith R, Titterington D M. Bayesian regression and classification using mixtures of Gaussian processes[J]. International Journal of Adaptive Control and Signal Processing, 2003, 17(2): 149-161.

23. Tayal A, Poupart P, Li Y. Hierarchical Double Dirichlet Process Mixture of Gaussian Processes[C]//AAAI. 2012.

24. Chatzis S P, Demiris Y. Nonparametric mixtures of Gaussian processes with power-law behavior[J]. IEEE transactions on neural networks and learning systems, 2012, 23(12): 1862-1871.

25. Kapoor A, Ahn H, Picard R W. Mixture of gaussian processes for combining multiple modalities[C]//International Workshop on Multiple Classifier Systems. Springer, Berlin, Heidelberg, 2005: 86-96.