Antony Lewis wrote:This paper investigates how many samples of full-sky signal [tex]C_l[/tex] estimators drawn from the posterior given the data you need to accurately estimate parameters. I find the result - that you need a (possibly) exponentially growing number of samples as you increase l_max - rather surprising (and disappointing for what seems like an otherwise very neat method).

Yes, the exponential function was rather surprising the first time we saw it, especially since, as you point out, for full sky coverage and no noise, you only need one single sample. The problem is that once you start to get coupling between different sigma_l's (see Fig. 6b of Eriksen et al. 2004, ApJS, 155, 227,

astro-ph/0407028) because of the sky cut, the total volume of the space you need to cover becomes huge with increasing l_max. Now, *we* know that the functional form of the likelihood is really very well behaved, and that perhaps only two- or three-point correlations are important, but the *sample-based estimator* doesn't know that, unless we explicitly tell it so (effectively set most of the likelihood space to zero).

In the very first run we made, we applied the BR estimator to the full range between l=2 and 400, and found that one single sigma-sample completely dominated the estimator. And that single sample was the one that by chance happened to have the smallest RMS power spectrum scatter around some imagined model. Now, it is true that this particular sample would have a very good chi^2, but on the other hand, it should also be assigned a very small volume in the likelihood space. But because we have so few samples, it isn't -- the chi^2 is still very good, but the volume (1/N_samples) is quite large. Therefore this single sample drives the fit.

So, the question is how to solve this problem. One approach is of course to generate an infinite number of samples, but that is a little difficult in practice [:-)]. Another is to impose further constraints on the estimator. One option is to split the likelihood space into bands, and apply the estimator to each of the bands, and multiply the partial likelihoods together. Unfortunately, this loses the inter-bin correlations. Another option is to pull out the univariate distributions from the BR estimator, and apply analytical approximate corrections to handle the correlations. There are lots of options for how to do this, and we are still working on this issue.

I still think the method is quite neat, although we currently do have problems with treating the full thing perfectly exactly. Being able to treat bands of delta_l=30 exactly with only 1000 samples isn't that bad, I think. And this is only the very first application of the method -- we will definitely refine it a lot in the future.

Antony Lewis wrote:
I'm also not clear whether *samples* in this paper correpsonds to *independent* samples, or *Gibbs* samples. If the latter, is the exponential bahaviour just the fact that more iterations are required to obtain independent samples?

For S/N >> 1, consequtive samples are independent, so the difference doesn't matter. You can have a look at Figure 8 of Eriksen et al. 2004, ApJS, 155, 227 to see the correlation length as a function of S/N. We also verified this in the current paper.

Antony Lewis wrote:
If one sampled directly from the parameters as a new Gibbs step as suggested in

astro-ph/0310080 rather than going via the Blackwell-Rao Estimator is the l-scaling problem avoided?

We haven't looked into this yet, and my intuition have failed earlier, so I don't think I want to make strong predictions quite yet. However, my gut feeling is that sampling over parameters should help quite a lot, since the dimensionality of the parameter space is so much smaller -- but don't take my word for it; there may be issues I haven't realized yet. So often there are... :-)