We present new Keck images at 0.9 μm and OVRO 1.3 mm continuum images of five Class I protostars in the Taurus star-forming region. We analyze these data in conjunction with broadband spectral energy distributions (SEDs) and 8-13 μm spectra from the literature using a Monte Carlo radiative transfer code. By fitting models for the circumstellar dust distributions simultaneously to the scattered light images, millimeter continuum data, and the SEDs, we attempt to distinguish between flared disks, infalling envelopes with outflow cavities, and combinations of disks and envelopes. For each of these circumstellar density distributions, we generate grids of models for varying geometries, dust masses, and accretion rates and determine the best fits by minimizing the residuals between model and data. Comparison of the residuals for best-fit disk, envelope, and disk+envelope models demonstrates that, in general, models incorporating both massive envelopes and massive embedded disks fit the imaging+SED data best. The implied envelope infall rates for these disk+envelope models are generally consistent with infall rates derived by previous investigators, although they are approximately an order of magnitude larger than inner disk accretion rates inferred from recent spectroscopic measurements. In addition, the disk masses inferred from our models are close to or larger than the limit for gravitationally stable disks, indicating that Class I disks may undergo periodic episodes of enhanced accretion, perhaps as a result of gravitational instabilities. An important caveat to these results is that in some cases, no single model can fit all of the imaging and SED data well, suggesting that further refinements to models of the circumstellar dust distributions around Class I sources are necessary. We discuss several potential improvements to the models, as well as new constraints that will become available with upcoming millimeter and infrared facilities.