  In Bayesian probabilistic programming, a central problem is to estimate the
normalised posterior distribution (NPD) of a probabilistic program with
conditioning via score (a.k.a. observe) statements. Most previous approaches
address this problem by Markov Chain Monte Carlo and variational inference, and
therefore could not generate guaranteed outcomes within a finite time limit.
Moreover, existing methods for exact inference either impose syntactic
restrictions or cannot guarantee successful inference in general.
  In this work, we propose a novel automated approach to derive guaranteed
bounds for NPD via polynomial solving. We first establish a fixed-point theorem
for the wide class of score-at-end Bayesian probabilistic programs that
terminate almost-surely and have a single bounded score statement at program
termination. Then, we propose a multiplicative variant of Optional Stopping
Theorem (OST) to address score-recursive Bayesian programs where score
statements with weights greater than one could appear inside a loop. Finally,
we use polynomial solving to implement our fixed-point theorem and OST variant.
To improve the accuracy of the polynomial solving, we further propose a
truncation operation and the synthesis of multiple bounds over various program
inputs. Our approach can handle Bayesian probabilistic programs with unbounded
while loops and continuous distributions with infinite supports. Experiments
over a wide range of benchmarks show that compared with the most relevant
approach (Beutner et al., PLDI 2022) for guaranteed NPD analysis via recursion
unrolling, our approach is more time efficient and derives comparable or even
tighter NPD bounds. Furthermore, our approach can handle score-recursive
programs which previous approaches could not.
