Given a sample of independent and identically distributed random variables, a novel nonparametric maximum entropy method is presented to estimate the underlying continuous univariate probability density function (pdf). Estimates are found by maximizing a log-likelihood function based on single order statistics after transforming through a sequence of trial cumulative distribution functions that iteratively improve using a Monte Carlo random search method. Improvement is quantified by assessing the random variables against the statistical properties of sampled uniform random data. Quality is determined using an empirically derived scoring function that is scaled to be sample size invariant. The scoring function identifies atypical fluctuations, for which threshold values are set to define objective criteria that prevent under-fitting as trial iterations continue to improve the model pdf, and, stopping the iteration cycle before over-fitting occurs. No prior knowledge about the data is required. An ensemble of pdf models is used to reflect uncertainties due to statistical fluctuations in random samples, and the quality of the estimates is visualized using scaled residual quantile plots that show deviations from size-invariant statistics. These considerations result in a tractable method that holistically employs key principles of random variables and their statistical properties combined with employing orthogonal basis functions and data-driven adaptive algorithms. Benchmark tests show that the pdf estimates readily converge to the true pdf as sample size increases. Robust results are demonstrated on several test probability densities that include cases with discontinuities, multi-resolution scales, heavy tails and singularities in the pdf, suggesting a generally applicable approach for statistical inference.