We present XookSuut, a Python implementation of the DiskFit algorithm, optimized to perform robust Bayesian inference on parameters describing models of circular and noncircular rotation in galaxies. XookSuut surges as a Bayesian alternative for kinematic modeling of 2D velocity maps; it implements effcient sampling methods, specifically Markov Chain Monte Carlo (MCMC) and Nested Sampling (NS), to obtain the posteriors and marginalized distributions of kinematic models including circular motions, axisymmetric radial flows, bisymmetric flows, and harmonic decomposition of the LoS velocity. In this way, kinematic models are obtained by pure sampling methods, rather than standard minimization techniques based on the Χ2. All together, XookSuut represents a sophisticated tool for deriving rotational curves and to explore the error distribution and covariance between parameters.