Abstract We report on the global structure of the Milky Way (MW) stellar halo up to its outer boundary based on the analysis of blue horizontal-branch stars (BHBs). These halo tracers are extracted from the (g, r, i, z)-band multi-photometry in the internal data release of the ongoing Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP) surveyed over a ∼550 deg2 area. In order to select the most likely BHBs by removing blue straggler stars (BSs) and other contamination in a statistically significant manner, we have developed and applied an extensive Bayesian method, instead of the simple color cuts adopted in our previous work, where each of the template BHBs and non-BHBs obtained from the available catalogs are represented as a mixture of multiple Gaussian distributions in the color–color diagrams. We found from the candidate BHBs in the range of 18.5 < g < 23.5 mag that the radial density distribution over a Galactocentric radius of r = 36–360 kpc can be approximated as a single power-law profile with an index of $α =3.74^{+0.21}_{-0.22}$ or a broken power-law profile with an index of $α _{\rm in}=2.92^{+0.33}_{-0.33}$ at r below a broken radius of $r_{\rm b}=160^{+18}_{-19}\:$kpc and a very steep slope of $α _{\rm out}=15.0^{+3.7}_{-4.5}$ at r > rb. The latter profile with a prolate shape having an axial ratio of $q=1.72^{+0.44}_{-0.28}$ is most likely and this halo may hold a rather sharp boundary at r ≃ 160 kpc. The slopes of the halo density profiles are compared with those from the suite of hydrodynamical simulations for the formation of stellar halos. This comparison suggests that the MW stellar halo may consist of the two overlapping components: the in situ inner halo as probed by RR Lyrae stars showing a relatively steep radial density profile and the ex situ outer halo with a shallow profile probed by BHBs here, which is made by accretion of small stellar systems.