This paper presents a method for detecting the ensemble means, spreads, and occurrence probabilities for each of the multiple Kuroshio states. This is accomplished by classifying the forecasts of the ensemble members with a Gaussian mixture distribution model, a machine learning method. Ensemble simulations with 80 members are conducted to reproduce possible occurrences of the multiple Kuroshio states, targeting the large meander event in 2017. To test its performance, first, the method is applied for the southernmost latitude, a conventional index that represents meander intensity. The results show that the Kuroshio initially taking the nearshore non-large meander state bifurcates into the large meander and offshore non-large meander states, which occur with the similar probabilities. Both developments are accompanied by positive potential energy extraction rates, consistent with baroclinic instability. As a more objective approach, the method is then applied for the dominant modes derived from empirical orthogonal function (EOF) analysis of the sea surface height field in the entire Kuroshio region. Importantly, almost identical results can be achieved. In particular, the bimodality between the large meander and non-large meander is shown to appear on the axis of the first EOF mode. From a mathematical perspective, this mode can be interpreted as the singular vector which grows most rapidly following the time-evolution operator. Finally, the multimodality of the Kuroshio is reinterpreted as a phase transition phenomenon where the nearshore non-large meander constitutes the basic state.