Context. The interstellar medium (ISM) is now widely acknowledged to display features ascribable to magnetized turbulence. With the public release of Planck data and the current balloon-borne and ground-based experiments, the growing amount of data tracing the polarized thermal emission from Galactic dust in the submillimetre provides choice diagnostics to constrain the properties of this magnetized turbulence. Aims: We aim to constrain these properties in a statistical way, focussing in particular on the power spectral index βB of the turbulent component of the interstellar magnetic field in a diffuse molecular cloud, the Polaris Flare. Methods: We present an analysis framework based on simulating polarized thermal dust emission maps using model dust density (proportional to gas density nH) and magnetic field cubes, integrated along the line of sight (LOS), and comparing these statistically to actual data. The model fields are derived from fractional Brownian motion (fBm) processes, which allows a precise control of their one- and two-point statistics. The parameters controlling the model are (1)-(2) the spectral indices of the density and magnetic field cubes, (3)-(4) the RMS-to-mean ratios for both fields, (5) the mean gas density, (6) the orientation of the mean magnetic field in the plane of the sky (POS), (7) the dust temperature, (8) the dust polarization fraction, and (9) the depth of the simulated cubes. We explore the nine-dimensional parameter space through a Markov chain Monte Carlo analysis, which yields best-fitting parameters and associated uncertainties. Results: We find that the power spectrum of the turbulent component of the magnetic field in the Polaris Flare molecular cloud scales with wavenumber as k-βB with a spectral index βB = 2.8 ± 0.2. It complements a uniform field whose norm in the POS is approximately twice the norm of the fluctuations of the turbulent component, and whose position angle with respect to the north-south direction is χ0 ≈-69°. The density field nH is well represented by a log-normally distributed field with a mean gas density ≈40 cm-3, a fluctuation ratio σnH/≈1.6, and a power spectrum with an index βn=1.7-0.3+0.4. We also constrain the depth of the cloud to be d ≈ 13 pc, and the polarization fraction p0 ≈ 0.12. The agreement between the Planck data and the simulated maps for these best-fitting parameters is quantified by a χ2 value that is only slightly larger than unity. Conclusions: We conclude that our fBm-based model is a reasonable description of the diffuse, turbulent, magnetized ISM in the Polaris Flare molecular cloud, and that our analysis framework is able to yield quantitative estimates of the statistical properties of the dust density and magnetic field in this cloud.