A method for analytic continuation of imaginary-time correlation functions (here obtained in quantum Monte Carlo simulations) to real-frequency spectral functions is proposed. Stochastically sampling a spectrum parametrized by a large number of delta-functions, treated as a statisticalmechanics problem, it avoids distortions caused by (as demonstrated here) configurational entropy in previous sampling methods. The key development is the suppression of entropy by constraining the spectral weight to within identifiable optimal bounds and imposing a set number of peaks. As a test case, the dynamic structure factor of the S = 1/2 Heisenberg chain is computed. Very good agreement is found with Bethe Ansatz results in the ground state (including a sharp edge) and with exact diagonalization of small systems at elevated temperatures. Obtaining real-frequency dynamic response functions from imaginary-time correlations remains one of the outstanding challenges for quantum Monte Carlo (QMC) and related simulation methods (e.g., lattice QCD). The general form of the problem is to invert the relationshipwhere a QMC estimateG(τ ) of the correlation function G(τ ) is available, A(ω) is the spectral function sought, and the kernel K(τ, ω) depends on the type of spectral function. Similar to an inverse Laplace transform, there is no closed form for A(ω). Only broad features of A(ω) can be resolved in numerical analytic continuation, because information on fine structure is only present at a level of precision of G(τ ) which is not attainable in practice. Nevertheless, one can extract important dynamical features and the key question is how to do that with the maximum fidelity, givenG(τ ) and its statistical errors. Significant progress will be presented here. The Maximum Entropy (ME) method [1] was adapted to the particulars of QMC some time ago [2]. Overcoming problems of previous approaches [3,4], it quickly became a standard tool [5]. The ME method has an appealing footing in probability theory, but in many cases the entropic prior regularizes the spectrum too heavily, leading to excessive broadening and distortions. To avoid this, an alternative line of methods has been developed [6][7][8][9][10] (and applied to diverse systems [11-14]) which do not impose the entropic prior, instead using stochastic sampling of A(ω) with the probability distributionwhere χ 2 is the standard measure of the goodness of the fit of G(τ ) obtained from A(ω) according to Eq. (1) to the QMC-computedG(τ ) with its full covariance matrix [5,8] for a set {τ i }. The spectrum is typically parametrized as a sum of a large number of δ-functions, though other forms have also been used [10]. The sampling temperature Θ in Eq. (2) acts as a regularizing parameter.An important insight was gained by Beach [7], showing that a mean-field treatment of the sampling approach gives the ME method, with Θ corresponding to the entropic weight. Subsequently, Syljuåsen argued for fixing Θ = 1 [8] (as had also been done by White in earlier work Here a previously overlook...