  We extend to quenched disordered systems the variational scheme for real
space renormalization group calculations that we recently introduced for
homogeneous spin Hamiltonians. When disorder is present our approach 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. Key to the variational approach is the bias
potential found by minimizing a convex functional in statistical mechanics.
This potential reduces dramatically the Monte Carlo relaxation time in large
disordered systems. We demonstrate the method with applications to the
two-dimensional dilute Ising model, the random transverse field quantum Ising
chain, and the random field Ising in two and three dimensional lattices.
