In this study a family of estimators is developed for local maxima, or modes, of a multivariate probability density function. The mode estimators are computationally feasible iterative optimization procedures utilizing nonparametric techniques of probability density estimation which generalize easily to sample spaces of arbitrary dimension. The estimators are proven to be strongly consistent for any distribution possessing mild continuity properties. Three specific mode estimators are evaluated by extensive Monte Carlo testing upon samples from both classical unimodal and nonstandard unimodal and biomodal distributions. Detection of the presence of multiple modes is a matter of special concern in many investigations. Thus a global strategy is developed and tested to demonstrate the potential of the estimators for complete characterization of sample modality.