Tomographic redshift distributions and systematics parameterisation

LSSTDESC / bayesian-pipelines-cosmology

Bayesian Cosmological Inference from Tomographic Maps

MIT License

4 stars 1 forks source link

Tomographic redshift distributions and systematics parameterisation #5

Open EiffL opened 2 years ago

EiffL commented 2 years ago

In this issue we aim to define our baseline choice for photoz distributions.

We should try to look at the specifications of the DC2 SRD choices: https://github.com/LSSTDESC/star-challenge/tree/main/cosmodc2-srd-sample/generation

EiffL commented 2 years ago

Thanks to @elts6570 we now have a demo notebook accessing the photoz for DC2 data \o/ https://github.com/LSSTDESC/bayesian-pipelines-cosmology/blob/main/photo-z/nz_h5_to_fits.ipynb

We are now looking with @aimalz at whether we need to parametrize n(z) systematics in the forward model. For lenses: For sources:

EiffL commented 2 years ago

Based on inputs from @aimalz, we should think about the n(z) provided by DESC PZ as posteriors n(z) with all possible n(z) information available. So, in our forward model, we should use those PZ-provided n(z) as fixed redshift distributions we use in our forward model.

EiffL commented 2 years ago

Here is a link to a notebook where I used the following function to represent an n(z) and implement a systematic shift: https://github.com/DifferentiableUniverseInitiative/jax-cosmo-paper/blob/master/notebooks/DES_Y1_shear_demo.ipynb

To create a redshift distribution object from an existing histogram:
```
import jax_cosmo as jc
```

nz = jc.redshift.kde_nz(z_hist, nz_hist, bw=0.01, gals_per_arcmin2=10)

The bandwidth parameter defines the amount of smoothing of the n(z)

To draw the modeled n(z):
```python
z = linspace(0,2)
plot(z, nz(z))

And to create a new nz object with a systemattic shift:

nz_sys = jc.redshift.systematic_shift(nz, dz)

where dz is the value of the shift

EiffL commented 2 years ago

Thanks @elts6570 for your notebook! I think we can condense the two into your first one, where we could overplot, 1. the original histogram, 2. the model redshift_distribution from the kde, 3. a version of these n(z) with a systematic shift.

I've tried to distill your code into the following blurb:


# Reading the DC2 tomographic bins into redshift distribution objects
with h5py.File("shear_photoz_stack.hdf5") as f:
    group = f["n_of_z"]
    # Read the z grid
    source = group["source"]
    z_shear = source['z'][::]
    # Read the true n(z)
    nz_shear = [jc.redshift.kde_nz(z_shear,  
                                   source[f"bin_{i}"][:], 
                                   bw=0.01, zmax=2.5) for i in range(5)] 

# Create systematically shifted versions of these n(z), 
# all with the same shift of 0.02 for testing
nzs_s_sys = [jc.redshift.systematic_shift(nzi, 0.02, 
                                            zmax=2.5) 
                for i, nzi in enumerate(nz_shear)]

for i in range(5):
  plot(z_shear, nz_shear[i](z_shear), color='C%d'%i, label="Bin %d"%i)
  plot(z_shear, nzs_s_sys[i](z_shear), '--', color='C%d'%i)
legend()
xlim(0,2);

or rsomething like thatt

elts6570 commented 2 years ago

Following up on the discussion about reconstructing pdfs from moments. There is a code that does this. It's called PyMaxEnt and can be found at PyMaxEnt. It does not have a pip package, so I downloaded src/pymaxent.py locally. I pushed the pdf_reconstruction.ipynb notebook on the photo-z branch, where I tried it on one of the DESC n(z)'s. It does not work well. It should be because of degeneracies with other distributions that have similar moments and will be preferred over the (more complicated) n(z)'s due to the nature of the method. At high moments, the method becomes more unstable.

Here are some references if you are interested into looking into the method further: