I would like to estimate the distribution of a very large population of known size but unknown mean and variance. I cannot assume anything about the underlying distribution. The values of observations are always discrete. I cannot sample the entire population but I would like to estimate the probability density function. I would also like some level of assurance of the correctness of the estimated distribution. What is the best way to go about this?

I can't comment so this will have to be in the form of an answer  my apologies for this. First, I assume that what you want to do is to estimate the distribution of some property of the population, no? One way to go about it would be to construct the empirical distribution function: $\mathbb{F}(x) = \frac{1}{n}\sum _{i=1} ^{n} I(X_i \leq x) $ ($n$ the number of observations). Once you have done this, you can use various test to compare this to any distribution function that you'd like (i.e., one that you suspect produced the data). If you are after a direct estimation of the density fucntion, and not any "known" distribution that it looks similar to, you can use a kernel density estimation: http://en.wikipedia.org/wiki/Kernel_density_estimation There are tools in e.g. Matlab that will make this task fairly straightforward. I would say that what method to use, and there are surely others, depends on that you want to do with the end result. I've only done this sort of thing in easy settings and there are probably more sophisticated ways to solve your problem. However, the above could perhaps set you off to a good start 

