  We have investigated the systematic differences introduced when performing a
Bayesian-inference analysis of the equation of state of neutron stars employing
either variable- or constant-likelihood functions. The former have the
advantage that it retains the full information on the distributions of the
measurements, making an exhaustive usage of the data. The latter, on the other
hand, have the advantage of a much simpler implementation and reduced
computational costs. In both approaches, the EOSs have identical priors and
have been built using the sound-speed parameterization method so as to satisfy
the constraints from X-ray and gravitational-waves observations, as well as
those from Chiral Effective Theory and perturbative QCD. In all cases, the two
approaches lead to very similar results and the $90\%$-confidence levels are
essentially overlapping. Some differences do appear, but in regions where the
probability density is extremely small and are mostly due to the sharp cutoff
set on the binary tidal deformability $\tilde \Lambda \leq 720$ employed in the
constant-likelihood analysis. Our analysis has also produced two additional
results. First, a clear inverse correlation between the normalized central
number density of a maximally massive star, $n_{\rm c, TOV}/n_s$, and the
radius of a maximally massive star, $R_{\rm TOV}$. Second, and most
importantly, it has confirmed the relation between the chirp mass
$\mathcal{M}_{\rm chirp}$ and the binary tidal deformability $\tilde{\Lambda}$.
The importance of this result is that it relates a quantity that is measured
very accurately, $\mathcal{M}_{\rm chirp}$, with a quantity that contains
important information on the micro-physics, $\tilde{\Lambda}$. Hence, once
$\mathcal{M}_{\rm chirp}$ is measured in future detections, our relation has
the potential of setting tight constraints on $\tilde{\Lambda}$.
