The rate of earthquakes in a region is a fundamental input to Probabilistic Seismic Hazard Analysis. We present a Monte Carlo method for computing that rate from seismicity catalogs while including a range of data and analysis uncertainties. This method is applied to regions for which the b value is assumed to be spatially invariant. Each region is broken down into epochs for which each epoch is estimated to have a uniform magnitude of completeness (⁠⁠M_c). The distribution of earthquake rates for M ≥ M_c is determined for each epoch by considering the Poisson likelihood of rates given the number of observed earthquakes with M ≥ M_c⁠. We use a Monte Carlo process to include the uncertainty in b, ⁠M_c, and individual event magnitudes. The result for each epoch is the joint distribution of the Poisson rate of earthquakes with magnitudes larger than the minimum value used to calculate hazard (⁠⁠M₁) and the Gutenberg–Richter b values, which control the extrapolation to other magnitudes. The rate for each region is either the duration‐weighted average over the epochs or, to better capture temporal variations, we also consider mixture models. The mixture models also provide an avenue to allow temporal variations in b values. To implement this joint distribution in a logic tree, we use the mean and 95% confidence branches, each of which is parameterized with an M ≥ M₁ rate and b value. We explore different ways of defining those branches, as well as non‐Gutenberg–Richter branches, and their impact on hazard estimates. The mean hazard, but not the fractiles, is robust with respect to these choices. To illustrate these new methods, we use synthetic data and catalogs from recent U.S. Geological Survey National Seismic Hazard Models for the Central and Eastern United States and for Puerto Rico and the U.S. Virgin Islands.