CUORE (Cryogenic Underground Observatory for Rare Events) is a tonne-scale cryogenic detector located at LNGS exploiting bolometric technique to search for neutrinoless double beta decay ($0 \nu \beta \beta$) of ${}^{130}$Te. The experimental signature is a sharp peak at the Q value in the summed energy spectrum of the emitted electrons. In this contribution, we present a thorough description of the software for the $0 \nu \beta \beta$ fit that was developed for the last data release. It is based on BAT (Bayesian Analysis Toolkit) and follows a Bayesian statistical approach. We show how the model is implemented, the treatment of systematics and the interpretation of results.