The Sloan Digital Sky Survey (SDSS) is the first dense redshift survey encompassing a volume large enough to find the best analytic probability density function that fits the galaxy Counts-in-Cells distribution $f_V(N)$, the frequency distribution of galaxy counts in a volume $V$. Different analytic functions have been previously proposed that can account for some of the observed features of the observed frequency counts, but fail to provide an overall good fit to this important statistical descriptor of the galaxy large-scale distribution. Our goal is to find the probability density function that better fits the observed Counts-in-Cells distribution $f_V(N)$. We have made a systematic study of this function applied to several samples drawn from the SDSS. We show the effective ways to deal with incompleteness of the sample (masked data) in the calculation of $f_V(N)$. We use LasDamas simulations to estimate the errors in the calculation. We test four different distribution functions to find the best fit: the Gravitational Quasi-Equilibrium distribution, the Negative Binomial Distribution, the Log Normal distribution and the Log Normal Distribution including a bias parameter. In the two latter cases, we apply a shot-noise correction to the distributions assuming the local Poisson model. We show that the best fit for the Counts-in-Cells distribution function is provided by the Negative Binomial distribution. In addition, at large scales the Log Normal distribution modified with the inclusion of the bias term also performs a satisfactory fit of the empirical values of $f_V(N)$. Our results demonstrate that the inclusion of a bias term in the Log Normal distribution is necessary to fit the observed galaxy Count-in-Cells distribution function.