Recent space plasma observations have provided us with three-dimensional velocity distributions having multiple peaks. We propose a method for analyzing such velocity distributions via a multivariate Maxwellian mixture model where each component of the model represents each of the multiple peaks. The parameters of the model are determined through an iterative nonlinear optimization technique, specifically the expectation-maximization (EM) algorithm. For the automatic judgment of the preferable number of components in the mixture model, we introduce a method of examining the number of extrema of a resulting mixture model. We show applications of our method to ion observations in the plasma sheet boundary layer (PSBL) and in the central plasma sheet (CPS) of the Earth's magnetotail. From an analysis of the PSBL and CPS, low-energy ions that have properties similar to those of lobe ions were detected also in the PSBL and CPS. In the PSBL, middle-energy ion component which is flowing dawnward and directed the neutral sheet was extracted. We suggested that magnetic field fluctuations in the PSBL can be explained only when the two ion components were properly treated. ¿ 2001 American Geophysical Union