We present a new method based on the N-point probability distribution (pdf) to study non-Gaussianity in cosmic microwave background (CMB) maps. Likelihood and Bayesian estimation are applied to a local non-linear perturbed model up to third order, characterized by a linear term which is described by a Gaussian N-pdf, and a second and third order terms which are proportional to the square and the cube of the linear one. We also explore a set of model selection techniques (the Akaike and the Bayesian Information Criteria, the minimum description length, the Bayesian Evidence and the Generalized Likelihood Ratio Test) and their application to decide whether a given data set is better described by the proposed local non-Gaussian model, rather than by the standard Gaussian temperature distribution. As an application, we consider the analysis of the WMAP 5-year data at a resolution of around 2 degrees. At this angular scale (the Sachs-Wolfe regime), the non-Gaussian description proposed in this work defaults (under certain conditions) to an approximative local form of the weak non-linear coupling inflationary model (e.g. Komatsu & Spergel 2001) previously addressed in the literature. For this particular case, we obtain an estimation for the non-linear coupling parameter of -94 < F_nl < 154 at 95% CL. Equally, model selection criteria also indicate that the Gaussian hypothesis is favored against the particular local non-Gaussian model proposed in this work. This result is in agreement with previous findings obtained for equivalent non-Gaussian models and with different non-Gaussian estimators. However, our estimator based on the N-pdf is more efficient than previous estimators and, therefore, provides tighter constraints on the coupling parameter at degree angular resolution.