  The discovery of the coalescence of binary neutron star GW170817 was a
watershed moment in the field of gravitational wave astronomy. Among the rich
variety of information that we were able to uncover from this discovery was the
first non-electromagnetic measurement of the neutron star radius, and the cold
nuclear equation of state. It also led to a large equation of state
model-selection study from gravitational-wave data. In those studies Bayesian
nested sampling runs were conducted for each candidate equation of state model
to compute their evidence in the gravitational-wave data. Such studies, though
invaluable, are computationally expensive and require repeated, redundant,
computation for any new models. We present a novel technique to conduct
model-selection of equation of state in an extremely rapid fashion (~minutes)
on any arbitrary model. We test this technique against the results of a
nested-sampling model-selection technique published earlier by the LIGO/Virgo
collaboration, and show that the results are in good agreement with a median
fractional error in Bayes factor of about 10%, where we assume that the true
Bayes factor is calculated in the aforementioned nested sampling runs. We found
that the highest fractional error occurs for equation of state models that have
very little support in the posterior distribution, thus resulting in large
statistical uncertainty. We then used this method to combine multiple binary
neutron star mergers to compute a joint-Bayes factor between equation of state
models. This is achieved by stacking the evidence of the individual events and
computing the Bayes factor from these stacked evidences for each pairs of
equation of state.
