We propose an approach called self-learning hybrid Monte Carlo (SLHMC), which is a general method to make use of machine learning potentials to accelerate the statistical sampling of first-principles density-functional-theory (DFT) simulations. The trajectories are generated on an approximate machine learning (ML) potential energy surface. The trajectories are then accepted or rejected by the Metropolis algorithm based on DFT energies. In this way, the statistical ensemble is sampled exactly at the DFT level for a given thermodynamic condition. Meanwhile, the ML potential is improved on the fly by training to enhance the sampling, whereby the training data set, which is sampled from the exact ensemble, is created automatically. Using the examples of alpha-quartz crystal SiO2 and phonon-mediated unconventional superconductor YNi2B2C systems, we show that SLHMC with artificial neural networks (ANN) is capable of very efficient sampling, while at the same time enabling the optimization of the ANN potential to within meV/atom accuracy. The ANN potential thus obtained is transferable to ANN molecular dynamics simulations to explore dynamics as well as thermodynamics. This makes the SLHMC approach widely applicable for studies on materials in physics and chemistry.