[Question] How to discretize the theory function

[steven-murray]

In principal, the observed power spectrum is discretized via

In practice, we will get a discretized window function, W, for which I presume the values have already been averaged in some way. This means we are applying the following approximation:

Where the k' integrals run over small cells in 2D k-space.

Now, the question remains: what to do about the integral over the theory function? Any number of integration rules could be applied -- use the midpoint, use the extremes of each k bin, use quadrature, etc. What should we support, and when does it make a difference?

The point is that getting the theory integral really accurate might not be that much of a win, because we might have made the bigger approximation above by separating the integrals over W and P anyway. The error from this step could possibly be looked into, and could inform us of what to do. In the meantime, we could support multiple options.

[steven-murray]

Tagging in @adeliegorce, @aewallwi and @Stefan-Heimersheim as interested parties.

[Stefan-Heimersheim]

I haven't looked into this for cylindrical power spectra, so far only for spherical power spectra which might be somewhat different (especially because there is no 2D integration then).

In that case I would expect that the integral will indeed not make too much of a difference in the results because the Delta^2 theory power spectra tend to be somewhat flat, at least for the larger values of k. (We are intending to do this for dimensionless power spectra, right?)

Edit: Forgot that P(k) isn't just the theory, but also foregrounds! So saying that Delta^2 is flat isn't necessarily true!

Edit: No, it's P(k) not Delta^2 so ignore my comments about it bring flat!

From the theory perspective, we usually get the theory code to compute the power spectrum at discrete points and then interpolate between them (linear in log P - log k space). This suggests that the integration over P(k) could be done analytically, and be reformatted as a sum over computed P(k) values with weights depending on the interpolation method and the distances between the computed values and the integration boundaries. In this case we would need to accept a set of theory values [k, P(k)] and handle the integration in the code.
Alternatively we can just pass the interpolation function (e.g. generated by scipy.interpolate) to the likelihood code and do a numerical integration internally as the interpolation function should be somewhat fast, although this will probably still introduce a bottleneck in the speed of the likelihood evaluation.

I think for speed and simplicity it would be best to not do this, and use a very simple rule (just use midpoint?). We could also offer one more-complex rule (some simple quadrature rule with a few points per dimension?) to double check that the correction is small.

[steven-murray]

Tagging in @acliu as well.