The chemical master equation (CME) is a mathematical tool utilized to model the stochasticity of the complex biochemical reaction networks. As the direct solution of the CME is notoriously expensive, moment-based approximations are computationally attractive in terms of time and memory to compute the statistics of the CME. The Bayesian method with delayed rejection adaptive Metropolis (DRAM) sampler provides an accurate probabilistic framework for parameter inference, reducing the need for costly computations of likelihood functions. In this study, we derive a system of ordinary differential equations (ODEs) using zero-moment approximations and employ the DRAM sampler algorithm to find the approximate posterior distributions of the model parameters. The quality of the method is evaluated through two systems biology case studies, with the aim of enabling efficient data-driven inferences for more complex biochemical models.