We present a technique to find guaranteed upper and lower bounds on the posterior distribution of probabilistic programs with loops. It relies on probability generating functions to represent distributions. Lower bounds can easily be obtained by Kleene iteration (unrolling the loops a finite number of times). We attack the harder problem of finding upper bounds by searching for an inductive invariant of a certain shape. This can be reduced to finding a solution to a system of polynomial inequalities (a decidable problem in the theory of the reals). Our method is provably sound and, if a solution can be found, it yields bounds on the probability masses, moments, and tail behavior of the program distribution. Our prototype implementation finds guaranteed bounds fully automatically for many examples from the literature that previously required human input to be solved.