Abstract
Large-scale sky surveys reveal noticeable spatial variation in stellar surface density, reflecting the structure of the Milky Way and the clustering of stellar populations. Star counts in sky regions often show substantial overdispersion due to the Galactic structure, unresolved substructure, and observational heterogeneity. Using publicly available Gaia DR3 data, we model star counts in 5∘×5∘ sky bins using the two-parameter Negative Binomial (NB2) distribution. This model handles the extra variation while keeping a clear relationship between the mean and variance. We develop profile-likelihood confidence intervals for the mean stellar density, both globally and as a function of Galactic latitude, which provide a statistically rigorous inference for the stellar surface density near and away from the Galactic plane. We compare our results with Wald confidence intervals, parametric bootstrap confidence intervals, and Poisson-based methods to show why it is important to account for extra variation in modern sky surveys. Our approach offers a reliable and practical way to model overdispersed star counts and can be used in many large-scale astronomy surveys.