The tremendous impact of Cosmic Microwave Background (CMB) radiation experiments on our understanding of the history and evolution of the universe is based on a tight connection between the observed fluctuations and the physical processes taking place in the very early universe. According to the prevalent paradigm, the anisotropies were generated during the era of inflation. The simplest inflationary models predict almost perfectly Gaussian primordial perturbations, but competitive theories can naturally be constructed, allowing for a wide range in primordial non-Gaussianity. For this reason, the test for non-Gaussianity becomes a fundamental means to probe the physical processes of inflation. The aim of the project is to develop a Bayesian formalism to infer the level of non-Gaussianity of local type. Bayesian statistics attaches great importance to a consistent formulation of the problem and properly calculates the error bounds of the measurements on the basis of the actual data. As a first step, we develop an exact algorithm to generate simulated temperature and polarization CMB maps containing arbitrary levels of local non-Gaussianity. We derive an optimization scheme that allows us to predict and actively control the simulation accuracy. Implementing this strategy, the code outperforms existing algorithms in computational efficiency by an order of magnitude. Then, we develop the formalism to extend the Bayesian approach to the calculation of the amplitude of non-Gaussianity. We implement an exact Hamiltonian Monte Carlo sampling algorithm to generate samples from the target probability distribution. These samples allow to construct the full posterior distribution of the level of non-Gaussianity given the data. The applicability of the scheme is demonstrated by means of a simplified data model. Finally, we fully implement the necessary equations considering a realistic CMB experiment dealing with partial sky coverage and anisotropic noise. A direct comparison between the traditional frequentist estimator and the exact Bayesian approach shows the advantage of the newly developed method. For a significant detection of non-Gaussianity, the former suffers from excess variance whereas the Bayesian scheme always provides optimal error bounds.