For precise application of Bhatnagar–Gross–Krook (BGK) methods, assessing its accuracy in non-equilibrium flows is necessary. Generally, this assessment relies on macroscopic parameters, which are moments of the velocity distribution function (VDF). However, in non-equilibrium flows, the significance of each moment diminishes as the VDF deviates from the Maxwellian VDF. This study investigates the VDF in non-equilibrium flows. Two Prandtl-corrected BGK methods, the ellipsoidal statistical BGK and Shakhov BGK (SBGK), are compared with the direct simulation Monte Carlo method. To observe the VDF while excluding the effects of convection, the homogeneous relaxation of the initial non-equilibrium state is analyzed. The VDF in Couette flow and normal shock waves, where collision and convection coexist, is then examined. When comparing the accuracy of the BGK methods using higher-order moments, inconsistencies are observed. However, when comparing the VDFs, the SBGK method reproduces the non-equilibrium VDF more accurately. The results demonstrate the importance of the VDF in the evaluation of non-equilibrium flows.