Background <p>Continuous traits evolution of a group of taxa that are correlated through a phylogenetic tree is commonly modelled using parametric stochastic differential equations to represent deterministic change of traits through time, while incorporating noises that represent different unobservable evolutionary pressures. A heterogeneous Gaussian process that consists of multiple parametric sub-processes is often used when the observed data come from a very diverse set of taxa. In the maximum likelihood setting, challenges arise when exploring the involved likelihood surface and when interpreting the uncertainty around the parameters.</p> Results <p>We extend the methods to tackle inference problems for mixed Gaussian phylogenetic models (MGPMs) by implementing a Bayesian scheme that can take into account biologically relevant priors. The posterior inference method is based on the Population Monte Carlo (PMC) algorithm that is easily parallelized, and uses an efficient algorithm to calculate the likelihood of phylogenetically correlated observations. A model evaluation method that is based on the proximity of the posterior predictive distribution to the observed data is also implemented. Simulation study is done to test the inference and evaluation capability of the method. Finally, we test our method on a real-world dataset.</p> Conclusion <p>We implement the method in the R package bgphy, available at <a href="https://github.com/bayubeta/bgphy">https://github.com/bayubeta/bgphy</a>. Simulation study demonstrates that the method is capable to infer parameters to evaluate different models, while its implementation on the real-world dataset indicates that a carefully selected model of evolution based on naturally occurring classifications results in a better fit to the observed data.</p>