[en] Searching for new physics in rare B meson decays governed by b → s transitions, we perform a model-independent global fit of the short-distance couplings C₇, C₉, and C₁₀ of the ΔB=1 effective field theory. We assume the standard-model set of b → sγ and b → sl⁺l^- operators with real-valued C_i. A total of 59 measurements by the experiments BaBar, Belle, CDF, CLEO, and LHCb of observables in B→K^*γ, B→K^(*)l⁺l^-, and B_s→μ⁺μ^- decays are used in the fit. Our analysis is the first of its kind to harness the full power of the Bayesian approach to probability theory. All main sources of theory uncertainty explicitly enter the fit in the form of nuisance parameters. We make optimal use of the experimental information to simultaneously constrain theWilson coefficients as well as hadronic form factors - the dominant theory uncertainty. Generating samples from the posterior probability distribution to compute marginal distributions and predict observables by uncertainty propagation is a formidable numerical challenge for two reasons. First, the posterior has multiple well separated maxima and degeneracies. Second, the computation of the theory predictions is very time consuming. A single posterior evaluation requires O(1s), and a few million evaluations are needed. Population Monte Carlo (PMC) provides a solution to both issues; a mixture density is iteratively adapted to the posterior, and samples are drawn in a massively parallel way using importance sampling. The major shortcoming of PMC is the need for cogent knowledge of the posterior at the initial stage. In an effort towards a general black-box Monte Carlo sampling algorithm, we present a new method to extract the necessary information in a reliable and automatic manner from Markov chains with the help of hierarchical clustering. Exploiting the latest 2012 measurements, the fit reveals a flipped-sign solution in addition to a standard-model-like solution for the couplings C_i. The two solutions are related by approximate invariance of the observables under C_i → -C_i. Both solutions contain about half of the posterior probability and provide a good fit to the data. The standard-model prediction is close to the global best-fit point. The Bayes factor strongly favors the simplicity of the standard model over the more complex new-physics scenario. For future searches, we compute two sets of predictions of observables in the angular distributions of B→K^*(→Kπ) l⁺l^-. This includes currently unmeasured optimized observables with reduced form-factor dependence and sensitivity to nonstandard interactions. In the first set, predictions within the standard model are calculated with a Bayesian approach and found to agree well with existing results. In the second set, we make use of the improved posterior knowledge of nuisance parameters to compute predictions based on the new-physics fit output. In the latter case, we observe significantly reduced theory uncertainty for all observables with form-factor dependence. Deviations from the predictions in future measurements of the predicted observables would clearly indicate new physics beyond the considered scenario in the standard-model operator basis.