The renormalization group is an important method to understand the critical behaviors of a statistical system. In this thesis, we develop a stochastic method to perform the renormalization group calculations non-perturbatively with Monte Carlo simulations. The method is variational in nature, and involves minimizing a convex functional of the renormalized Hamiltonians. The variational scheme overcomes critical slowing down, by means of a bias potential that renders the coarse-grained variables uncorrelated. When quenched disorder is present in the statistical system, the method gives access to the flow of the renormalized Hamiltonian distribution, from which one can compute the critical exponents if the correlations of the renormalized couplings retain finite range. The bias potential again reduces dramatically the Monte Carlo relaxation time in large disordered systems. With this method, we also demonstrate how to extract the higher-order geometrical information of the critical manifold of a system, such as its tangent space and curvature. The success of such computations attests to the existence and robustness of the renormalization group fixed-point Hamiltonians. In the end, we extend the method to continuous-time quantum Monte Carlo simulations, which allows an accurate determination of the sound velocity of the quantum system at criticality. In addition, a lattice energy-stress tensor emerges naturally, where the continuous imaginary-time direction serves as a ruler of the length scale of the system.