We give a method to compute guaranteed estimates of Bayesian a pos- teriori distributions in a model where the relation between the observation y and the parameters θ is a function, possibly involving additive noise parameters ψ, say y = f(θ) + h(ψ). This model covers the case of (noisy) ode parameters estimation and the case when f is computed by a neural network. Applying a combination of methods based on uncertain probability (P-boxes), Interval Arith- metic (IA) and Monte Carlo (MC) simulation, we design an efficient randomized algorithm that returns guaranteed estimates of the posterior CDF of the parameters θ, and moments thereof, given that the observation y lies in a (small) rectangle. Guarantees come in the form of confidence intervals for the CDF values and its moments. Comparison with state-of-the-art approaches on odes benchmarks shows significant improvement in terms of efficiency and accuracy.