  Turbulence is a dominant feature operating in gaseous flows across nearly all
scales in astrophysical environments. Accordingly, accurately estimating the
statistical properties of such flows is necessary for developing a
comprehensive understanding of turbulence. We develop and employ a hierarchical
Bayesian fitting method to estimate the parameters describing the scaling
relationships of the velocity power spectra of supersonic turbulence. We
demonstrate the accuracy and other advantages of this technique compared with
ordinary linear regression methods. Using synthetic power spectra, we show that
the Bayesian method provides accurate parameter and error estimates. Commonly
used normal linear regression methods can provide estimates that fail to
recover the underlying slopes, up to 70% of the instances, even when
considering the 2std.dev. uncertainties. Additionally, we apply the Bayesian
methods to analyse the statistical properties of compressible turbulence in 3D
numerical simulations. We model driven, isothermal, turbulence with rms Mach
numbers in the highly supersonic regime M~15. We study the influence of purely
solenoidal (divergence-free) and purely compressive (curl-free) forcing on the
scaling exponent of the power spectrum. In simulations with solenoidal forcing
and 1024^3 resolution, our results indicate that there is no extended inertial
range with a constant scaling exponent. The bottleneck effect results in a
curved power spectrum at all wave numbers and is more pronounced in the
transversal modes compared with the longitudinal modes. Therefore, this effect
is stronger in stationary turbulent flows driven by solenoidal forcing compared
to the compressive one. The longitudinal spectrum driven with compressive
forcing is the only spectrum with constant scaling exponent z=-1.94 +- 0.01,
corresponding to slightly shallower slopes than the Burger's prediction.
