  The precision and predictive power of perturbative QCD (pQCD) prediction
depends on both a precise, convergent fixed-order series and a reliable way of
estimating the contributions of unknown higher-order (UHO) terms. It has been
shown that by applying the Principal of Maximum Conformality (PMC), which
applies the renormalization group equation recursively to set the effective
magnitude of $\alpha_s$ of the process, the remaining conformal coefficients
will be well matched with the corresponding $\alpha_s$ at each orders, leading
to a scheme-and-scale invariant and convergent perturbative series. Thus
different from conventional scheme-and-scale dependent fixed-order series, the
PMC series will provide a more reliable platform for estimating UHO
contributions. In this paper, by using the total decay width
$\Gamma(H\to\gamma\gamma)$ which has been calculated up to N$^4$LO QCD
corrections, we derive its PMC series by using the PMC single-scale setting
approach and estimate its unknown N$^5$LO contributions by using the Bayesian
analysis. The Bayesian-based approach estimates the magnitude of the UHO
contributions based on an optimized analysis of probability density
distribution, and the predicted UHO contribution becomes more accurate when
more loop terms have been known to tame the probability density function. Using
the top-quark pole mass $M_t$=172.69 GeV and the Higgs mass $M_H$=125.25 GeV as
inputs, we obtain $\Gamma(H\to\gamma\gamma) =9.56504~{\rm keV}$ and the
estimated N$^5$LO contribution to the total decay width is
$\Delta\Gamma_H=\pm1.65\times10^{-4}~{\rm keV}$ for the smallest credible
interval of $95.5\%$ degree-of-belief.
