"Shift" parameterization

[aimalz]

@EiffL has provided a use case for the shift parameterization in the context of the BP TT. @tilmantroester is this consistent with what you had in mind?

[tilmantroester]

Yes - I was writing up and soliciting feedback a detailed use case and context:

This issue concerns the needs of firecrown to marginalise over n(z) uncertainties during sampling. If the described functionality does not yet exist in qp, there is a discussion to be had if the functionality we have in mind should be implemented in qp or in firecrown.

I’m not familiar with the qp terminology, so to avoid confusion with existing terms, I shall call the object that holds the information on the n(z) and its uncertainties just the qp object.

At each sampling step, firecrown ultimately needs a tuple of arrays (z, dndz) that gets passed to CCL to generate Cls with. Naively, one could randomly choose a different realisation of the n(z) at each sampling step to account for the n(z) uncertainty. Certain samplers struggle with such stochastic likelihoods (the likelihood changes even though the sampled parameters do not). For this reason, as well as for backward comparability, it is required to be able to marginalise over the n(z) uncertainty in a parametric fashion. That is, given a (small) set of continuous parameters, deterministically generate realisations of the n(z) such that marginalising over these parameters accounts for the full uncertainty of the n(z). A simplistic example are the current shift parametrisations, where the generated n(z) is just shifted in z by some amount delta_z.

The parameters that govern the generation of the n(z) should ideally be interpreted in terms of the uncertainty of the n(z). So a delta_z = 1 could mean “generate a n(z) that corresponds to the fiducial n(z) shifted to higher z by one sigma, where sigma refers to the uncertainty of the mean of the fiducial n(z).”

The basic API for accounting for n(z) systematics in firecrown at the moment is basically
new_nz = generate_nz(nz_params, nz)
. Here nz_params is a dict of parameters that parametrise how the new n(z) gets generated and nz is an object that holds the information about the n(z), which is either created from the sacc file or from in a previous systematics stage.
For details, have a look at the current shifting implementation here: https://github.com/LSSTDESC/firecrown/blob/master/firecrown/likelihood/gauss_family/statistic/source/weak_lensing.py#L181

For the simple shift parameterisation, we have something like (based on the firecrown implementation)

def generate_nz(nz_params, nz):
    dndz_interp = Akima1DInterpolator(nz.z_array, nz.nz_array)

    new_nz_array = dndz_interp(z_array - nz_params["delta_z"], extrapolate=False)
    new_nz_array[np.isnan(new_nz_array)] = 0.0

    return NZ(new_nz_array, z_array)

What I have in mind to get there is

1. Read n(z) data (for example, n(z) realisations) from a sacc file and generate a qp object at the start of the run.
2. During each MCMC sampling step, the qp object deterministically generates a new n(z), based on some continuos parameters.

In pseudo-code, for the simplistic shifting approach (we would eventually want something more sophisticated):

# During setup at the beginning of the run
qp = QP.from_sacc(sacc_file, tracer_name)

# ...

# During sampling
nz_array, z_array = qp.parametric_sample(mode="shift", delta_z=0.42)

# ...

[aimalz]

Another place where we can experiment with error parameterizations is this code related to the BP TT's STAR efforts. The suggestion in today's telecon was to use the means and standard deviations of a Gaussian mixture model, which is already entirely implemented in qp. (So I don't forget, @jecampagne may have additional code to enhance that example.) @EiffL had some additional thoughts to add in a follow-up comment.

[jecampagne]

@aimalz the code pointing behind "this code" is the code we have used ( @EiffL and me ) to perform a DES Y1 3x2pts exercise. So I'm not sure that I have something different from that.

[tilmantroester]

There are two aspects to this. One is to deterministically sample n(z)'s, based on some parameters. In the BP example, that'd be the jc.redshift.systematic_shift method.

The other aspect is to define priors on these parameters. These priors should be informed by the uncertainty of the n(z). Ideally, the firecrown sampler doesn't need to know about the specifics of these priors and can just sample from a ~fixed distribution, e.g. uniformly between 0 and 1.

I mocked this up here: https://colab.research.google.com/drive/1ie1ZIR-iuAr8owBgh7s5Ap0LPUE1xtOr?usp=sharing
The parameters sampled by firecrown would be u~U(0, 1), which then get transformed internally into normally distributed shifts, with mean 0 and std being the scatter in the mean of the bootstrap n(z) samples.

The question is how much of this functionality can be abstracted away in qp and what should be implemented on the firecrown side. The shift parametrisation is easy enough to deal with but it gets more complicated with more sophisticated approaches.

[aimalz]

This issue is subsumed by #99 so I'm closing this one but linking it there so the discussion is still accessible.