We investigate the semi-leptonic decays of $bar B to D^{(*)} ellbar u$ in terms of the Heavy-Quark-Effective-Theory (HQET) parameterization for the form factors, which is described with the heavy quark expansion up to $mathcal 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 $|V_{cb}|$ 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 $|V_{cb}|$ 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 $text{SM} + text{NP}$ scenarios, our fit analyses find that non-zero NP contribution is favored at the best fit point for both $text{SM} + V_2$ and $text{SM} + T$ depending on the HQET parameterization model. A key feature is then realized in the $bar B to D^{(*)} taubar u$ observables. Our fit result of the HQET parameters in the $text{SM} (+T)$ produces a consistent value for $R_D$ while smaller for $R_{D^*}$, compared with the previous SM prediction in the HFLAV report. On the other hand, $text{SM}+V_2$ points to smaller and larger values for $R_D$ and $R_{D^*}$ than the SM predictions. In particular, the $R_{D^*}$ deviation from the experimental measurement becomes smaller, which could be interesting for future improvement on measurements at the Belle II experiment.