  We present a flexible framework for constructing physical models of quasar
evolution that can incorporate a variety of observational constraints, such as
multi-wavelength luminosity functions, estimated masses and accretion rates of
active black holes, space densities of quasar hosts, and the local black hole
mass function. We focus on the accretion rate distribution p(mdot|M,z), the
probability that a black hole of mass M at redshift z accretes at a rate mdot
in Eddington units. Given the radiative efficiency as a function of mdot, the
quasar luminosity function (QLF) is determined by a convolution of p(mdot|M,z)
with the black hole mass function n(M,z). In the absence of mergers p(mdot|M,z)
also determines the full evolution of n(M,z), given a "boundary value" of n(M)
at some redshift. Matching the observed decline of the QLF break luminosity at
z<2 requires either a shift towards lower characteristic accretion rates or an
evolving mass dependence of p(mdot) that preferentially shuts off accretion
onto high mass black holes. These two scenarios make different predictions for
the masses and accretion rates of active black holes. If the first mechanism
dominates, then the QLF changes character between z=2 and z=0, shifting from a
sequence of black hole mass towards a sequence of L/Ledd. We construct and
compare five models that illustrate different assumptions about the quasar
population: short and long lifetime models dominated by unobscured thin-disk
accretion, a model with a high fraction of obscured quasars, and models in
which mergers or ADAF accretion produce substantial black hole growth at low
redshift. We discuss the observational advances that would be most valuable for
distinguishing such models and for pinning down the physics that drives black
hole and quasar evolution. (Abridged)
