Using fast lattice Monte Carlo (FLMC) simulations both in a canonical ensemble and with Wang-Landau–transition-matrix sampling, we have studied a model system of laterally homogeneous homopolymer brushes in an explicit solvent. Direct comparisons of the simulation results with those from the corresponding lattice self-consistent field (LSCF) theory, both of which are based on the same Hamiltonian (thus without any parameter-fitting between them), unambiguously and quantitatively reveal the fluctuations and correlations in the system. We have examined in detail how the Flory–Huggins interaction parameter χ between polymer segments and solvent molecules and the number of grafted chains n affect both the brush structures and thermodynamics. For our model system, the LSCF theory is exact in the limit of χ → −∞, except that it neglects the correlations among solvent molecules caused by the incompressibility constraint (thus overestimating the solvent entropy). At finite n and χ, the segmental density profile in the direction perpendicular to the grafting substrate obtained from FLMC simulations is flatter than the LSCF prediction, and the free-end density from FLMC simulations is also lower than the LSCF prediction close to the substrate. At finite n and χ > 0, the LSCF theory overestimates the internal energy. In addition, at small χ ≥ 0 it underestimates the difference in free energy but overestimates that in entropy from the reference state in the limit of χ → −∞, and at larger χ the opposite occurs. At large n, FLMC results approach LSCF predictions at a rate of 1/n in most cases.