This paper explores the use of the substitution method to derive the most probable distribution in the context of teaching statistical physics, serving as a complementary approach to the traditional Lagrange multiplier method. By explicitly incorporating constraints into the thermodynamic probability and extremizing it, the equilibrium distributions for classical, Bose, and Fermi systems are derived. Although computationally more involved, this method offers a clearer geometric interpretation of constraints and reveals more physical details, helping students better understand the mathematical nature of extremization problems and their connection to statistical distributions. The second-order condition confirming entropy maximization is also verified, underscoring the method’s foundational and pedagogical value.