In this paper, we study the problem of learning dynamical properties of ensemble systems from their collective behaviors using statistical approaches in reproducing kernel Hilbert space (RKHS). Specifically, we provide a framework to identify and cluster multiple ensemble systems through computing the maximum mean discrepancy (MMD) between their aggregated measurements in an RKHS, without any prior knowledge of the system dynamics of ensembles. Then, leveraging on a gradient flow of the newly proposed notion of aggregated Markov parameters, we present a systematic framework to recognize and identify an ensemble systems using their linear approximations. Finally, we demonstrate that the proposed approaches can be extended to cluster multiple unknown ensembles in RKHS using their aggregated measurements. Numerical experiments show that our approach is reliable and robust to ensembles with different types of system dynamics.