  Measurement of macroscopic properties of neutron stars, whether in binary or
in an isolated system, provides us a key opportunity to place a stringent
constraint on its equation of state. In this paper, we perform Bayesian
model-selection on a wide variety of neutron star equation of state using
multi-messenger observations. In particular, (i) we use the mass and tidal
deformability measurement from two binary neutron star merger event, GW170817
and GW190425; (ii) simultaneous mass-radius measurement of PSR J0030+0451 and
PSR J0740+6620 by NICER collaboration, while the latter has been analyzed by
joint NICER/radio/XMM-Newton collaboration. Among the 31 equations of state
considered in this analysis, we are able to rule out different variants of MS1
family, SKI5, H4, and WFF1 EoSs decisively, which are either extremely stiff or
soft equations of state. The most preferred equation of state model turns out
to be AP3 (or MPA1), which predicts the radius and dimensionless tidal
deformability of a $1.4 M_{\odot}$ neutron star to be 12.10 (12.50) km and 393
(513) respectively.
