  In this work, we formalize a new technique to investigate joint posterior
density of Cosmic Microwave Background (CMB) signal and its theoretical angular
power spectrum given the observed data, using the global
internal-linear-combination (ILC) method first proposed by Sudevan & Saha
(2017). We implement the method on low resolution CMB maps observed by WMAP and
Planck satellite missions, using Gibbs sampling, assuming that the detector
noise is negligible on large angular scales of the sky. The main products of
our analysis are best fit CMB cleaned map and its theoretical angular power
spectrum along with their error estimates. We validate the methodology by
performing Monte Carlo simulations that includes realistic foreground models
and noise levels consistent with WMAP and Planck observations. Our method has
an unique advantage that the posterior density is obtained without any need to
explicitly model foreground components. Secondly, the power spectrum results
with the error estimates can be directly used for cosmological parameter
estimations.
