Intrinsic stochasticity plays an essential role in gene regulation because of the small number of involved molecules of DNA, mRNA and protein of a given species. To better understand this phenomenon, small gene regulatory systems are mathematically modeled as systems of coupled chemical reactions, but the existing exact description utilizing a Chapman-Kolmogorov equation or simulation algorithms is limited and inefficient. The present work introduces a much more efficient yet accurate modeling approach, which allows analyzing stochasticity in the system in terms of the underlying distribution function. The novel modeling approach is motivated by the analysis of a single gene regulatory module with three sources of stochasticity: intermittent gene activity, mRNA transcription/decay and protein translation/decay noise. Although the corresponding Chapman Kolmogorov equation cannot be solved when a large number of molecules are considered, it is used to analytically derive the first two moments of the underlying distribution function. The mRNA and protein variance is found decomposable into additive terms resulting from the respective sources of stochasticity, which allow quantifying their significance in the process. The variance decomposition is asserted by constructing two approximations that establish a novel modeling approach: First, the continuous approximation, which considers only the stochasticity due to the intermittent gene activity. Second, the mixed approximation, which in addition attributes stochasticity to the mRNA transcription/decay process. Introduced approximations yield systems of first order partial differential equations for the underlying distribution function, which can be efficiently solved using developed numerical methods. Single cell simulations and numerical two-dimensional mRNA-protein stationary distribution functions are presented to confirm accuracy of introduced models. Further simplifications in the model allow considering regulation of the two- (possibly three-) gene systems for which two-dimensional protein-protein distributions are calculated. Finally, the assumption that gene activity is due to the binding and dissociation of a single regulatory molecule is relaxed. Based on the gene expression data, the models developed are applied to hypothesize the existence of a sequential activation mechanism of NF-kappaB dependent genes important in cell survival and inflammation. Future applications include analysis of small genetic networks, which are being currently engineered based on the prokaryotic and eukaryotic components.