Many recent analyses have indicated that large scale WMAP data display anomalies that appear inconsistent with the standard cosmological paradigm. However, the effects of foreground contamination, which require elimination of some fraction of the data, have not been carefully investigated due to the complexity in the analysis. Here we develop a general formalism of how to incorporate these effects in any analysis of this type. Our approach is to compute the full multi-dimensional probability distribution function of all possible sky realizations that are consistent with the data and with the allowed level of contamination. As an example we apply this method to compute the joint probability distribution function for the possible realizations of quadrupole and octopole using the WMAP data. This 12 dimensional distribution function is explored using the Markov Chain Monte Carlo technique. The resulting chains are used to asses the statistical significance of the low quadrupole using frequentist methods, the quadrupole-octopole alignment using several methods (angular momentum dispersion, multipole vectors and a new method based of feature matching) and the quadrupole-octopole aligment with ecliptic.