We investigate the semi-leptonic decays of B→D^(∗)ℓν in terms of the Heavy-Quark-Effective-Theory (HQET) parameterization for the form factors, which is described with the heavy quark expansion up to O(1/m_c^2) beyond the simple approximation considered in the original CLN parameterization. An analysis with this setup was first given in the literature, and then we extend it to the comprehensive analyses including (i) simultaneous fit of ∣Vcb∣ and the HQET parameters to available experimental full distribution data and theory constraints, and (ii) New Physics (NP) contributions of the V_2 and T types to the decay distributions and rates. For this purpose, we perform Bayesian fit analyses by using Stan program, a state-of-the-art public platform for statistical computation. Then, we show that our ∣Vcb∣ fit results for the SM scenarios are close to the PDG combined average from the exclusive mode, and indicate significance of the angular distribution data. In turn, for the SM+NP scenarios, our fit analyses find that non-zero NP contribution is favored at the best fit point for both SM+V_2 and SM+T depending on the HQET parameterization model. A key feature is then realized in the B→D^(∗)τν observables. Our fit result of the HQET parameters in the SM produces a consistent value for RD while smaller for RD^∗, compared with the previous SM prediction in the HFLAV report. On the other hand, SM+V_2 points to smaller and larger values for RD and RD^∗ than the SM predictions. In particular, the RD^* deviation from the experimental measurement becomes smaller, which could be interesting for future improvement on measurements at the Belle II experiment.