Several cosmological measurements have attained significant levels of
maturity and accuracy over the last decade. Continuing this trend, future
observations promise measurements of the statistics of the cosmic mass
distribution at an accuracy level of one percent out to spatial scales with
k~10 h/Mpc and even smaller, entering highly nonlinear regimes of gravitational
instability. In order to interpret these observations and extract useful
cosmological information from them, such as the equation of state of dark
energy, very costly high precision, multi-physics simulations must be
performed. We have recently implemented a new statistical framework with the
aim of obtaining accurate parameter constraints from combining observations
with a limited number of simulations. The key idea is the replacement of the
full simulator by a fast emulator with controlled error bounds. In this paper,
we provide a detailed description of the methodology and extend the framework
to include joint analysis of cosmic microwave background and large scale
structure measurements. Our framework is especially well-suited for upcoming
large scale structure probes of dark energy such as baryon acoustic
oscillations and, especially, weak lensing, where percent level accuracy on
nonlinear scales is needed.