In many applications, the success probability 𝜋 of a Bernoulli distributed variable 𝑌 is influenced by another variable 𝑋. For example for loans granted, it is necessary to rate debtors in different rating classes, where the probability of default (PD) 𝜋 of 𝑌 is assumed to be homogeneous within and heterogeneous between the rating classes. The PD of a debtor is largely influenced by macroeconomic and individual variables (𝑋). In this work, we study a Bernoulli mixture model for 𝑌, where the success probability of 𝑌 changes systematically at the breakpoints. We focus on cross-sectional data and our main objective is to estimate all 𝑘 breakpoints with 𝑘 either known or unknown and their corresponding success probabilities between each pair of neighbouring breakpoints. To the best of our knowledge, an estimator for estimating multiple breakpoints has not yet been developed in this context. Thus, we develop an approach with a view to closing this research gap. We show that our estimator works for independent and identically distributed (i.i.d.) 𝑋 as well as for a linear one-factor model for 𝑋. A theoretical foundation for this estimator is also presented. In practice, the number of breakpoints 𝑘 is often unknown a priori. As the multiple estimator is based on an iterative procedure, we propose stopping criteria for estimating 𝑘 correctly. We conduct a simulation study in the context of credit rating to demonstrate the performance of the developed estimator. Furthermore, we apply the new estimator on credit risk data from the Sächsische Aufbaubank, the Development Bank of Saxony. To simplify the use of the new estimator, we also develop an R package called MultipleBreakpoints.