We develop a multiscale approach to estimate high-dimensional probability
distributions from a dataset of physical fields or configurations observed in
experiments or simulations. In this way we can estimate energy functions (or
Hamiltonians) and efficiently generate new samples of many-body systems in
various domains, from statistical physics to cosmology. Our method -- the
Wavelet Conditional Renormalization Group (WC-RG) -- proceeds scale by scale,
estimating models for the conditional probabilities of "fast degrees of
freedom" conditioned by coarse-grained fields. These probability distributions
are modeled by energy functions associated with scale interactions, and are
represented in an orthogonal wavelet basis. WC-RG decomposes the microscopic
energy function as a sum of interaction energies at all scales and can
efficiently generate new samples by going from coarse to fine scales. Near
phase transitions, it avoids the "critical slowing down" of direct estimation
and sampling algorithms. This is explained theoretically by combining results
from RG and wavelet theories, and verified numerically for the Gaussian and
${\phi}^{4}$ field theories. We show that multiscale WC-RG energy-based models
are more general than local potential models and can capture the physics of
complex many-body interacting systems at all length scales. This is
demonstrated for weak-gravitational-lensing fields reflecting dark matter
distributions in cosmology, which include long-range interactions with
long-tail probability distributions. WC-RG has a large number of potential
applications in non-equilibrium systems, where the underlying distribution is
not known {\it a priori}. Finally, we discuss the connection between WC-RG and
deep network architectures.