  The Maximum Entropy Spectral Analysis (MESA) method, developed by Burg,
offers a powerful tool for spectral estimation of a time-series. It relies on
Jaynes' maximum entropy principle, allowing the spectrum of a stochastic
process to be inferred using the coefficients of an autoregressive process
AR($p$) of order $p$. A closed-form recursive solution provides estimates for
both the autoregressive coefficients and the order $p$ of the process. We
provide a ready-to-use implementation of this algorithm in a Python package
called \texttt{memspectrum}, characterized through power spectral density (PSD)
analysis on synthetic data with known PSD and comparisons of different criteria
for stopping the recursion. Additionally, we compare the performance of our
implementation with the ubiquitous Welch algorithm, using synthetic data
generated from the GW150914 strain spectrum released by the LIGO-Virgo-Kagra
collaboration. Our findings indicate that Burg's method provides PSD estimates
with systematically lower variance and bias. This is particularly manifest in
the case of a small (O($5000$)) number of data points, making Burg's method
most suitable to work in this regime. Since this is close to the typical length
of analysed gravitational waves data, improving the estimate of the PSD in this
regime leads to more reliable posterior profiles for the system under study. We
conclude our investigation by utilising MESA, and its particularly easy
parametrisation where the only free parameter is the order $p$ of the AR
process, to marginalise over the interferometers noise PSD in conjunction with
inferring the parameters of GW150914.
