Accurate phase diagram calculation from molecular dynamics requires systematic treatment and convergence of statistical averages. In this work we propose a Gaussian process regression based framework for reconstructing the free-energy functions using data of various origins. Our framework allows for propagating statistical uncertainty from finite molecular dynamics trajectories to the phase diagram and automatically performing convergence with respect to simulation parameters. Furthermore, our approach provides a way for automatic optimal sampling in the simulation parameter space based on a Bayesian optimization approach. We validate our methodology by constructing phase diagrams of two model systems, the Lennard-Jones and soft-core potential, and compare the results with the existing studies and our coexistence simulations. Finally, we construct the phase diagram of lithium at temperatures above 300 K and pressures below 30 GPa from a machine-learning potential trained on ab initio data. Our approach performs well when compared to coexistence simulations and experimental results.