diff --git a/README.md b/README.md index 65f5eea..0d02db1 100644 --- a/README.md +++ b/README.md @@ -30,6 +30,6 @@ As a disclaimer, it will need a major upgrade for flexibility and computational ## License, Contributing etc The code in this repo is available for re-use under the MIT license, which means that you can do whatever you like with it, just don't blame me. -If you end up using any of the code or ideas you find here in your academic research, please cite me as `Malz et al, in preparation\footnote{\texttt{https://github.com/aimalz/chippr}}`. +If you end up using any of the code or ideas you find here in your academic research, please cite me as `Malz & Hogg, 2020 \footnote{[arXiv:2007.12178](https://arxiv.org/abs/2007.12178)\texttt{https://github.com/aimalz/chippr}}`. If you are interested in this project, please do drop me a line via the hyperlinked contact name above, or by [writing me an issue](https://github.com/aimalz/chippr/issues/new). To get started contributing to the `chippr` project, just fork the repo -- pull requests are always welcome! diff --git a/chippr/__init__.py b/chippr/__init__.py index b0f6c29..5e1c823 100644 --- a/chippr/__init__.py +++ b/chippr/__init__.py @@ -1,17 +1,17 @@ -from defaults import * +from .defaults import * -from utils import * -from stat_utils import * -from plot_utils import * +from .utils import * +from .stat_utils import * +from .plot_utils import * -from sim_utils import * -from discrete import * -from gauss import * -from gamma import * -from gmix import * -from mvn import * -from multi_dist import * -from catalog import * +from .sim_utils import * +from .discrete import * +from .gauss import * +from .gamma import * +from .gmix import * +from .mvn import * +from .multi_dist import * +from .catalog import * -from log_z_dens import * -from log_z_dens_plots import * +from .log_z_dens import * +from .log_z_dens_plots import * diff --git a/chippr/catalog.py b/chippr/catalog.py index 2dbda40..d7aaf2f 100644 --- a/chippr/catalog.py +++ b/chippr/catalog.py @@ -5,7 +5,7 @@ import pickle as pkl import matplotlib as mpl -mpl.use('PS') +# mpl.use('PS') import matplotlib.pyplot as plt import chippr @@ -45,7 +45,7 @@ def __init__(self, params={}, vb=True, loc='.', prepend=''): self.params = d.check_sim_params(self.params) if vb: - print self.params + print(self.params) np.random.seed(d.seed) self.cat = {} @@ -427,7 +427,7 @@ def evaluate_lfs(self, pspace, vb=True): lfs = [] for n in self.N_range: points = zip(self.z_fine, [self.samps[n][1]] * self.n_tot) - cur=pspace.pdf(np.array(points)) + cur = pspace.pdf(np.array([p for p in points])) lfs.append(cur) lfs = np.array(lfs) lfs /= np.sum(lfs, axis=-1)[:, np.newaxis] * self.dz_fine diff --git a/chippr/catalog_plots.py b/chippr/catalog_plots.py index 83540c4..1f74e66 100644 --- a/chippr/catalog_plots.py +++ b/chippr/catalog_plots.py @@ -2,7 +2,7 @@ import os import matplotlib as mpl -mpl.use('PS') +# mpl.use('PS') import matplotlib.pyplot as plt from mpl_toolkits.axes_grid1 import make_axes_locatable @@ -31,8 +31,8 @@ def plot_true_histogram(true_samps, n_bins=(10, 50), plot_loc='', prepend='', pl pu.set_up_plot() f = plt.figure(figsize=(5, 5)) sps = f.add_subplot(1, 1, 1) - sps.hist(true_samps, bins=n_bins[1], normed=1, color='k', alpha=0.5, log=True) - sps.hist(true_samps, bins=n_bins[0], normed=1, color='y', alpha=0.5, log=True) + sps.hist(true_samps, bins=n_bins[1], color='k', alpha=0.5, log=True)#, normed=1 + sps.hist(true_samps, bins=n_bins[0], color='y', alpha=0.5, log=True)#, normed=1 sps.set_xlabel(r'$z_{true}$') sps.set_ylabel(r'$n(z_{true})$') f.savefig(os.path.join(plot_loc, prepend+plot_name), bbox_inches='tight', pad_inches = 0, dpi=d.dpi) @@ -105,6 +105,9 @@ def plot_mega_scatter(zs, pfs, z_grid, grid_ends, truth=None, plot_loc='', prepe int_pr: numpy.ndarray, float, optional plot the interim prior with the histograms? """ + print('debug demo') + print((truth, np.shape(zs), np.shape(pfs), np.shape(z_grid), np.shape(grid_ends))) + n = len(zs) zs = zs.T true_zs = zs[0] @@ -154,10 +157,10 @@ def plot_mega_scatter(zs, pfs, z_grid, grid_ends, truth=None, plot_loc='', prepe histy.hist(obs_zs, bins=grid_ends, orientation='horizontal', alpha=0.25, color='k', stacked=False, density=True) if truth is not None: histx.step(truth[0], truth[1], c='k', where='mid') - pu.plot_step(histy, np.pad(truth[1], (1, 1), 'constant', constant_values=(0, 0)), grid_ends, c='k', w=1.5) + pu.plot_step(histy, np.pad(truth[1], (1, 1), 'constant', constant_values=(0, 0)), u.ends_from_mids(truth[0]), c='k', w=1.5) if int_pr is not None: histx.step(int_pr[0], int_pr[1], c='k', alpha=0.5, where='mid') - pu.plot_step(histy, np.pad(int_pr[1], (1, 1), 'constant', constant_values=(0, 0)), grid_ends, c='k', a=0.5, w=1.5) + pu.plot_step(histy, np.pad(int_pr[1], (1, 1), 'constant', constant_values=(0, 0)), u.ends_from_mids(int_pr[0]), c='k', a=0.5, w=1.5) histx.xaxis.set_tick_params(labelbottom=False) histy.yaxis.set_tick_params(labelleft=False) histx.set_yticks([]) diff --git a/chippr/discrete.py b/chippr/discrete.py index c04d5d4..e700487 100644 --- a/chippr/discrete.py +++ b/chippr/discrete.py @@ -30,7 +30,7 @@ def __init__(self, bin_ends, weights): self.bin_ends = bin_ends self.dbins = self.bin_ends[1:] - self.bin_ends[:-1] self.n_bins = len(self.bin_ends)-1 - self.bin_range = range(self.n_bins) + self.bin_range = list(range(self.n_bins)) self.weights = weights # print('dbins: '+str(self.dbins)) diff --git a/chippr/log_z_dens.py b/chippr/log_z_dens.py index 810c2df..3bbb003 100644 --- a/chippr/log_z_dens.py +++ b/chippr/log_z_dens.py @@ -2,11 +2,11 @@ import scipy as sp import os import scipy.optimize as op -import cPickle as cpkl +import pickle as cpkl import emcee import matplotlib as mpl -mpl.use('PS') +# mpl.use('PS') import matplotlib.pyplot as plt import chippr @@ -61,7 +61,7 @@ def __init__(self, catalog, hyperprior, truth=None, loc='.', prepend='', vb=True self.info['log_interim_posteriors'] = self.log_pdfs if vb: - print(str(self.n_bins) + ' bins, ' + str(len(self.log_pdfs)) + ' interim posterior PDFs') + print((str(self.n_bins) + ' bins, ' + str(len(self.log_pdfs)) + ' interim posterior PDFs')) self.hyper_prior = hyperprior @@ -214,12 +214,12 @@ def _objective(log_nz): return -2. * self.evaluate_log_hyper_posterior(log_nz) if vb: - print(self.dir + ' starting at ', start, _objective(start)) + print((self.dir + ' starting at ', start, _objective(start))) res = op.minimize(_objective, start, method="Nelder-Mead", options={"maxfev": 1e5, "maxiter":1e5}) if vb: - print(self.dir + ': ' + str(res)) + print((self.dir + ': ' + str(res))) return res.x def calculate_mmle(self, start, vb=True, no_data=0, no_prior=0): @@ -428,7 +428,7 @@ def distribution(log_nz): full_chain = np.array([[vals[w]] for w in range(self.n_walkers)]) while self.burning_in: if vb: - print('beginning sampling '+str(self.burn_ins)) + print(('beginning sampling '+str(self.burn_ins))) burn_in_mcmc_outputs = self.sample(vals, 10**n_burned) chain = burn_in_mcmc_outputs['chains'] burn_in_mcmc_outputs['chains'] -= u.safe_log(np.sum(np.exp(chain) * self.bin_difs[np.newaxis, np.newaxis, :], axis=2))[:, :, np.newaxis] @@ -489,10 +489,10 @@ def compare(self, vb=True): self.info['stats']['rms']['true_nz' + '__' + key[4:]] = s.calculate_rms(np.exp(self.info['log_tru_nz']), np.exp(self.info['estimators'][key])) self.info['stats']['log_rms']['log_true_nz'+ '__' + key] = s.calculate_rms(self.info['log_tru_nz'], self.info['estimators'][key]) - for i in range(len(self.info['estimators'].keys())): - key_1 = self.info['estimators'].keys()[i] - for j in range(len(self.info['estimators'].keys()[:i])): - key_2 = self.info['estimators'].keys()[j] + for i in range(len(list(self.info['estimators'].keys()))): + key_1 = list(self.info['estimators'].keys())[i] + for j in range(len(list(self.info['estimators'].keys())[:i])): + key_2 = list(self.info['estimators'].keys())[j] # print(((i,j), (key_1, key_2))) self.info['stats']['log_rms'][key_1 + '__' + key_2] = s.calculate_rms(self.info['estimators'][key_1], self.info['estimators'][key_2]) self.info['stats']['rms'][key_1[4:] + '__' + key_2[4:]] = s.calculate_rms(np.exp(self.info['estimators'][key_1]), np.exp(self.info['estimators'][key_2])) @@ -537,11 +537,11 @@ def read(self, read_loc, style='pickle', vb=True): with open(os.path.join(self.res_dir, read_loc), 'rb') as file_location: self.info = cpkl.load(file_location) if vb: - print('The following quantities were read from '+read_loc+' in the '+style+' format:') + print(('The following quantities were read from '+read_loc+' in the '+style+' format:')) for key in self.info: print(key) if 'estimators' in self.info: - print(self.info['estimators'].keys()) + print((list(self.info['estimators'].keys()))) return self.info def write(self, write_loc, style='pickle', vb=True): @@ -560,7 +560,7 @@ def write(self, write_loc, style='pickle', vb=True): with open(os.path.join(self.res_dir, write_loc), 'wb') as file_location: cpkl.dump(self.info, file_location) if vb: - print('The following quantities were written to '+write_loc+' in the '+style+' format:') + print(('The following quantities were written to '+write_loc+' in the '+style+' format:')) for key in self.info: print(key) return diff --git a/chippr/log_z_dens_plots.py b/chippr/log_z_dens_plots.py index 93471f9..190e448 100644 --- a/chippr/log_z_dens_plots.py +++ b/chippr/log_z_dens_plots.py @@ -3,7 +3,7 @@ import scipy as sp import matplotlib as mpl -mpl.use('PS') +# mpl.use('PS') import matplotlib.pyplot as plt from matplotlib import gridspec diff --git a/chippr/mvn.py b/chippr/mvn.py index 69292cf..4e65208 100644 --- a/chippr/mvn.py +++ b/chippr/mvn.py @@ -29,7 +29,7 @@ def __init__(self, mean, var): self.sigma = self.norm_var() self.invvar = self.invert_var() - assert np.linalg.eig(self.var) > 0. + assert(np.all(np.linalg.eig(self.var)[0] > 0.)) self.dist = MGD(self.mean, self.var) diff --git a/chippr/plot_utils.py b/chippr/plot_utils.py index 42c224c..d0a0278 100644 --- a/chippr/plot_utils.py +++ b/chippr/plot_utils.py @@ -1,7 +1,7 @@ import numpy as np import matplotlib as mpl -mpl.use('PS') +# mpl.use('PS') import matplotlib.pyplot as plt import matplotlib.cm as cm diff --git a/chippr/stat_utils.py b/chippr/stat_utils.py index 57d9b09..7491218 100644 --- a/chippr/stat_utils.py +++ b/chippr/stat_utils.py @@ -151,7 +151,7 @@ def multi_parameter_gr_stat(sample): """ dims = np.shape(sample) (n_walkers, n_iterations, n_params) = dims - n_burn_ins = n_iterations / 2 + n_burn_ins = int(n_iterations / 2) chain_ensemble = np.swapaxes(sample, 0, 1) chain_ensemble = chain_ensemble[n_burn_ins:, :] Rs = np.zeros((n_params)) @@ -177,7 +177,7 @@ def gr_test(sample, threshold=d.gr_threshold): True if burning in, False if post-burn in """ gr = multi_parameter_gr_stat(sample) - print('Gelman-Rubin test statistic = '+str(gr)) + print(('Gelman-Rubin test statistic = '+str(gr))) test_result = np.max(gr) > threshold return test_result @@ -198,7 +198,7 @@ def cft(xtimes, lag):#xtimes has ntimes elements autocorrelation time for one time lag for one parameter of one walker """ lent = len(xtimes) - lag - allt = xrange(lent) + allt = range(lent) ans = np.array([xtimes[t+lag] * xtimes[t] for t in allt]) return ans @@ -217,8 +217,8 @@ def cf(xtimes):#xtimes has ntimes elements autocorrelation time over all time lags for one parameter of one walker """ cf0 = np.dot(xtimes, xtimes) - allt = xrange(len(xtimes) / 2) - cf = np.array([sum(cft(xtimes,lag)[len(xtimes) / 2:]) for lag in allt]) / cf0 + allt = range(int(len(xtimes) / 2)) + cf = np.array([sum(cft(xtimes, lag)[int(len(xtimes) / 2):]) for lag in allt]) / cf0 return cf def cfs(x, mode):#xbinstimes has nbins by ntimes elements diff --git a/chippr/utils.py b/chippr/utils.py index 69f0834..e39cf42 100644 --- a/chippr/utils.py +++ b/chippr/utils.py @@ -56,3 +56,39 @@ def ingest(in_info): else: in_dict = in_info return in_dict + +def mids_from_ends(inarr): + """ + Function to make midpoints from grid of endpoints + + Parameters + ---------- + inarr: array, shape=(N) + array of grid endpoints + + Returns + ------- + outarr: array, shape=(N-1) + array of corresponding midpoints + """ + outarr = (inarr[1:] + inarr[:-1]) / 2. + return outarr + +def ends_from_mids(inarr): + """ + Function to make endpoints from grid of midpoints, assuming equal spacing + + Parameters + ---------- + inarr: array, shape=(N) + array of grid midpoints + + Returns + ------- + outarr: array, shape=(N+1) + array of corresponding endpoints + """ + dif = np.mean(inarr[1:] - inarr[:-1]) + int_ends = mids_from_ends(inarr) + outarr = np.concatenate((np.array([int_ends[0] - dif]), int_ends, np.array([int_ends[-1] + dif]))) + return outarr diff --git a/docs/notebooks/demo2.ipynb b/docs/notebooks/demo2.ipynb index 0a4f977..5267974 100644 --- a/docs/notebooks/demo2.ipynb +++ b/docs/notebooks/demo2.ipynb @@ -8,7 +8,9 @@ "\n", "This notebook demonstrates the use of the Cosmological Hierarchical Inference with Probabilistic Photometric Redshifts (CHIPPR) package to estimate population distributions based on a catalog of probability distributions.\n", "\n", - "The package supports two primary objectives: simulation of catalogs and inference of posterior distributions over parameters defining population distributions." + "The package supports two primary objectives: simulation of catalogs and inference of posterior distributions over parameters defining population distributions.\n", + "\n", + "*Note that this notebook dates back to the [Python 2 tagged release](https://github.com/aimalz/chippr/releases/tag/v0.1-beta) and may not work under the latest version.*" ] }, { @@ -451,14 +453,14 @@ "language_info": { "codemirror_mode": { "name": "ipython", - "version": 2 + "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.15rc1" + "pygments_lexer": "ipython3", + "version": "3.8.5" } }, "nbformat": 4, diff --git a/docs/notebooks/demo3.ipynb b/docs/notebooks/demo3.ipynb new file mode 100644 index 0000000..953a9fe --- /dev/null +++ b/docs/notebooks/demo3.ipynb @@ -0,0 +1,499 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# CHIPPR\n", + "\n", + "This notebook demonstrates the use of the Cosmological Hierarchical Inference with Probabilistic Photometric Redshifts (CHIPPR) package to estimate population distributions based on a catalog of probability distributions.\n", + "\n", + "The package supports two primary objectives: simulation of catalogs and inference of posterior distributions over parameters defining population distributions." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import timeit\n", + "# import cProfile, pstats, StringIO\n", + "import os" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib as mpl\n", + "# print(mpl.rcParams.items)\n", + "mpl.use('Agg')\n", + "mpl.rcParams['text.usetex'] = False\n", + "mpl.rcParams['mathtext.rm'] = 'serif'\n", + "mpl.rcParams['font.family'] = 'serif'\n", + "mpl.rcParams['font.serif'] = ['Times New Roman']\n", + "# mpl.rcParams['font.family'] = ['Times New Roman']\n", + "mpl.rcParams['axes.titlesize'] = 25\n", + "mpl.rcParams['axes.labelsize'] = 20\n", + "mpl.rcParams['xtick.labelsize'] = 15\n", + "mpl.rcParams['ytick.labelsize'] = 15\n", + "mpl.rcParams['savefig.dpi'] = 250\n", + "mpl.rcParams['figure.dpi'] = 250\n", + "mpl.rcParams['savefig.format'] = 'pdf'\n", + "mpl.rcParams['savefig.bbox'] = 'tight'\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import chippr\n", + "help(chippr)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Simulation\n", + "\n", + "Many of `chippr`'s modules are used to produce mock catalogs of individual posterior probability distributions." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To create a catalog, we must first define a true redshift distribution function. It may be a Gaussian distribution of the `gauss` class, a Gaussian mixture distribution of the `gmix` class, or a binned discrete distribution of the `discrete` class. In this case, we will consider a mixture of three Gaussian distributions to represent the true redshift distribution." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "true_amps = np.array([0.20, 0.35, 0.55])\n", + "true_means = np.array([0.5, 0.2, 0.75])\n", + "true_sigmas = np.array([0.4, 0.2, 0.1])\n", + "n_mix_comps = len(true_amps)\n", + "\n", + "true_funcs = []\n", + "for c in range(n_mix_comps):\n", + " true_funcs.append(chippr.gauss(true_means[c], true_sigmas[c]**2))\n", + "\n", + "true_nz = chippr.gmix(true_amps, true_funcs, limits=(0., 1.))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`chippr` supports the use of a parameter file to specify various options for the catalog simulator and inference module to turn on and off." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "param_loc = 'params.txt'\n", + "params = chippr.utils.ingest(param_loc)\n", + "params = chippr.defaults.check_sim_params(params)\n", + "print(params)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To make a catalog, we must specify an interim prior redshift distribution $n^{*}(z)$, regardless of what quantity we wish to infer using the catalog. So far, only discrete distributions are supported, but this will soon be changed. The simplest discrete distribution is uniform." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "bin_ends = np.array([params['bin_min'], params['bin_max']])\n", + "weights = np.array([1.])\n", + "\n", + "int_prior = chippr.discrete(bin_ends, weights)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we are ready to make a catalog. To do this we instantiate a `catalog` object and then create a catalog of indiviual posterior distributions based on the true redshift distribution and the interim prior. By default, the catalog generator will make some informative plots along the way. The included plots are a histogram of the true redshifts and a scatterplot of the true redshifts and the centers of the individual Gaussian posteriors. Support for other posterior forms will be added soon. Additionally, the catalog is expressed as normalized binned histogram heights. Support for other parametrizations of the individual posteriors may be added in the future." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(chippr.utils.ends_from_mids(bin_ends))\n", + "print(chippr.utils.mids_from_ends(bin_ends))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "results_loc = os.path.join(os.path.join(os.path.join(os.path.join('..', '..'), 'research'), 'results'), 'demo')\n", + "\n", + "posteriors = chippr.catalog(params=param_loc, loc=results_loc)\n", + "\n", + "output = posteriors.create(true_nz, int_prior, N=params['n_gals'])\n", + "\n", + "data = np.exp(output['log_interim_posteriors'])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can also plot a histogram of the centers of the individual Gaussian posteriors, a binned version of the true redshift distribution, and the $n(z)$ resulting from stacking the individual posteriors." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.hist(posteriors.samps.T[1], bins=100, color=\"k\")#, normed=True\n", + "plt.plot(posteriors.z_coarse, true_nz.evaluate(posteriors.z_coarse), \"r-\")\n", + "plt.plot(posteriors.z_coarse, np.sum(data, axis=0) / 10**params['n_gals'], \"go\")\n", + "plt.xlabel(\"z\")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "It is also informative to see what a few individual likelihoods and binned posteriors look like." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for n, z in enumerate(data[:10]):\n", + " plt.plot(posteriors.z_coarse, data[n], 'ko')\n", + " plt.plot(posteriors.z_fine, posteriors.obs_lfs[n], 'k-')\n", + " plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We finish by saving the data as a plaintext file. Support for more file formats will be added soon." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "saved_location = 'data'\n", + "saved_type = '.txt'\n", + "posteriors.write(loc=saved_location, style=saved_type)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Inference\n", + "\n", + "`chippr` currently contains one inference module to probe the posterior distribution of parameters defining the redshift distribution function." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To perform inference, we must create a catalog object. This may be done by making a new catalog as is done above or by reading in an existing catalog file." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "param_loc = 'params.txt'\n", + "results_loc = os.path.join(os.path.join(os.path.join(os.path.join('..', '..'), 'research'), 'results'), 'demo')\n", + "simulated_posteriors = chippr.catalog(params=param_loc, loc=results_loc)\n", + "\n", + "saved_location = 'data'\n", + "saved_type = '.txt'\n", + "data = simulated_posteriors.read(loc=saved_location, style=saved_type)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The catalog file contains three components: the `bin_ends`, the `log_interim_prior`, and the `log_interim_posteriors`. The bin endpoints can be processed to enable their use in constructing a prior distribution over the parameters determining the redshift distribution function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "zs = data['bin_ends']\n", + "log_nz_intp = data['log_interim_prior']\n", + "log_z_posts = data['log_interim_posteriors']\n", + "\n", + "z_difs = zs[1:]-zs[:-1]\n", + "z_mids = (zs[1:]+zs[:-1])/2.\n", + "n_bins = len(z_mids)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The prior distribution must be a `mvn` object, defined by a mean vector and covariance matrix over the parameters defining the redshift distribution. In this case, it is intuitive to use the definition of the binning strategy to create the prior distribution since the parameters are normalized histogram bin heights, the same parametrization used for the catalog entries themselves." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# prior_sigma = 0.16\n", + "# prior_var = np.eye(n_bins)\n", + "# for b in range(n_bins):\n", + "# prior_var[b] = 1. * np.exp(-0.5 * (z_mids[b] - z_mids) ** 2 / prior_sigma ** 2)\n", + "# l = 1.e-4\n", + "# prior_var = prior_var+l*np.identity(n_bins)\n", + "\n", + "\n", + "a = 1.# / n_bins\n", + "b = 20.#1. / z_difs ** 2\n", + "c = 0.\n", + "prior_var = np.eye(n_bins)\n", + "for k in range(n_bins):\n", + " prior_var[k] = a * np.exp(-0.5 * b * (z_mids[k] - z_mids) ** 2)\n", + "prior_var += c * np.identity(n_bins)\n", + "\n", + "prior_mean = log_nz_intp\n", + "prior = chippr.mvn(prior_mean, prior_var)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We create a `log_z_dens` object from the dictionary of catalog parameters and the prior distribution. We include the optional specification of the true distribution, since it is available in this case. We also include a parameter file that may contain default constants for the inference. Now the `log_z_dens` object can plot some samples from the prior so we can see what it looks like." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nz = chippr.log_z_dens(data, prior, truth=true_nz, vb=True)\n", + "\n", + "prior_samples = prior.sample(7)\n", + "chippr.log_z_dens_plots.plot_ivals(prior_samples, nz.info, nz.plot_dir)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We perform calculations of a few of the simplest estimators of the redshift distribution function $\\hat{n}(z)$. The stacked estimator is defined as $\\hat{n}(z)=\\frac{1}{N}\\sum p(z|\\vec{d},n^{*}(z))$. The marginalized maximum a posteriori estimator is defined as $\\hat{n}(z)=\\hat{n}(\\{argmax[p(z|\\vec{d},n^{*}(z))]\\})$. The marginalized expected value estimator is defined as $\\hat{n}(z)=\\hat{n}(\\{E[p(z|\\vec{d},n^{*}(z))]\\})$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nz_stacked = nz.calculate_stacked()\n", + "nz_mmap = nz.calculate_mmap()\n", + "nz_mexp = nz.calculate_mexp()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The `log_z_dens` object enables easy comparison between estimators using the Kullback-Leibler Divergences (when the true distribution is available) and root-mean-square differences." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We may next calculate the marginalized maximum likelihood estimator (which actually returns the parameters maximizing the posterior probability)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nz_mmle = nz.calculate_mmle(nz_stacked)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If we are very ambitious, we can run an MCMC sampler (currently use of `emcee` is supported, but other samplers may be added in the future) to probe the posterior distribution of the parameter values. To do this, we initialize the sampler with samples from the prior distribution." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "n_ivals = 2 * n_bins\n", + "initial_values = prior.sample(n_ivals)\n", + "\n", + "nz_samps = nz.calculate_samples(initial_values)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nz_stats = nz.compare()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The `log_z_dens` object stores the estimators that have been calculated as well as all metadata associated with the posterior samples. The storage of the metadata and samples will soon be eliminated in favor of saved files, as that information may necessitate a great deal of memory." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nz.info['estimators'].keys()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Currently, the results of all previously calculated estimators (and the true redshift density function, if it was provided) may be plotted automatically." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "nz.plot_estimators()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The `log_z_dens` object supports writing the information associated with the estimators to a file in the `pickle` format, though other formats may be added in the future." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nz.write('nz.p')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here we demonstrate that the written estimators may be loaded from files as well for future use." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nz.info = nz.read('nz.p')\n", + "print(nz)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "chippr (Python 3)", + "language": "python", + "name": "chippr_3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +} diff --git a/research/paper/binned_nz_recovered.png b/research/paper/binned_nz_recovered.png new file mode 100644 index 0000000..1a3e0f2 Binary files /dev/null and b/research/paper/binned_nz_recovered.png differ diff --git a/research/paper/combo_mega_scatter.png b/research/paper/combo_mega_scatter.png new file mode 100644 index 0000000..034ddf4 Binary files /dev/null and b/research/paper/combo_mega_scatter.png differ diff --git a/research/paper/figures/chippr/thesis_neghivarbias-mega_scatter.png b/research/paper/data_bias.png similarity index 100% rename from research/paper/figures/chippr/thesis_neghivarbias-mega_scatter.png rename to research/paper/data_bias.png diff --git a/research/paper/figures/chippr/single_lsst-mega_scatter.png b/research/paper/data_lsst.png similarity index 100% rename from research/paper/figures/chippr/single_lsst-mega_scatter.png rename to research/paper/data_lsst.png diff --git a/research/paper/figures/chippr/thesis_eout-mega_scatter.png b/research/paper/data_outlier_template.png similarity index 100% rename from research/paper/figures/chippr/thesis_eout-mega_scatter.png rename to research/paper/data_outlier_template.png diff --git a/research/paper/figures/chippr/thesis_rout-mega_scatter.png b/research/paper/data_outlier_training.png similarity index 100% rename from research/paper/figures/chippr/thesis_rout-mega_scatter.png rename to research/paper/data_outlier_training.png diff --git a/research/paper/figures/chippr/single_uout-mega_scatter.png b/research/paper/data_outlier_uniform.png similarity index 100% rename from research/paper/figures/chippr/single_uout-mega_scatter.png rename to research/paper/data_outlier_uniform.png diff --git a/research/paper/figures/chippr/single_lsst_tmpr-mega_scatter.png b/research/paper/data_prior_template.png similarity index 100% rename from research/paper/figures/chippr/single_lsst_tmpr-mega_scatter.png rename to research/paper/data_prior_template.png diff --git a/research/paper/figures/chippr/single_lsst_trpr-mega_scatter.png b/research/paper/data_prior_training.png similarity index 100% rename from research/paper/figures/chippr/single_lsst_trpr-mega_scatter.png rename to research/paper/data_prior_training.png diff --git a/research/paper/figures/chippr/thesis_hivarsig-mega_scatter.png b/research/paper/data_scatter_hi.png similarity index 100% rename from research/paper/figures/chippr/thesis_hivarsig-mega_scatter.png rename to research/paper/data_scatter_hi.png diff --git a/research/paper/figures/chippr/single_varsigmas-mega_scatter.png b/research/paper/data_scatter_lo.png similarity index 100% rename from research/paper/figures/chippr/single_varsigmas-mega_scatter.png rename to research/paper/data_scatter_lo.png diff --git a/research/paper/draft.pdf b/research/paper/draft.pdf index a5c3126..9a4a059 100644 Binary files a/research/paper/draft.pdf and b/research/paper/draft.pdf differ diff --git a/research/paper/draft.tex b/research/paper/draft.tex index 49372d1..959a494 100644 --- a/research/paper/draft.tex +++ b/research/paper/draft.tex @@ -1,32 +1,22 @@ -%\RequirePackage[]{lineno} - \documentclass[iop]{emulateapj} \usepackage{tikz} \usepackage{natbib} \usepackage{amsmath} \usepackage{hyperref} -%\usepackage{graphicx} -%\usepackage{lineno} +\usepackage{lineno} \usepackage[percent]{overpic} \usepackage{float} +\usepackage[normalem]{ulem} \usepackage{wrapfig} -%\usetikzlibrary{shapes.geometric, arrows} -%\usetikzlibrary{fit} -% -%\tikzstyle{hyper} = [circle, text centered, draw=black] -%\tikzstyle{param} = [circle, text centered, draw=black] -%\tikzstyle{data} = [circle, text centered, draw=black, line width=2pt] -%\tikzstyle{arrow} = [thick,->,>=stealth] - -%\usepackage{soul} \usepackage[title]{appendix} \newcommand{\todo}[3]{{\color{#2}\emph{#1}: #3}} \newcommand{\aim}[1]{\todo{AIM}{red}{#1}} \newcommand{\dwh}[1]{\todo{Attn. Hogg}{blue}{#1}} \newcommand{\que}[1]{\todo{Question}{cyan}{#1}} +\newcommand{\new}[2]{{\color{red}\sout{#1}}\ {\color{blue}{#2}}}%{#2}% \newcommand{\myemail}{aimalz@nyu.edu} \newcommand{\mathul}[1]{\underline{#1}} @@ -42,50 +32,50 @@ \newcommand{\sdss}{\project{SDSS}} \newcommand{\boss}{\project{BOSS}} \newcommand{\des}{\project{DES}} -\newcommand{\Chippr}{\project{CHIPPR}}% maybe change this to algorithm/italics +\newcommand{\Chippr}{\project{CHIPPR}} \newcommand{\repo}[1]{\texttt{#1}} \newcommand{\qp}{\repo{qp}} \newcommand{\chippr}{\repo{chippr}} \newcommand{\cosmolike}{\repo{CosmoLike}} \newcommand{\emcee}{\repo{emcee}} -\newcommand{\github}{\href{https://github.com}{GitHub}}% maybe remove link and use code font -\newcommand{\python}{\textit{Python}}% establish code font? +\newcommand{\github}{\href{https://github.com}{GitHub}} +\newcommand{\python}{\textit{Python}} -\newcommand{\data}{\ensuremath{\vec{d}}}% could change to bold +\newcommand{\data}{\ensuremath{\vec{d}}} \newcommand{\like}{\mathscr{L}} -\newcommand{\pr}[1]{\ensuremath{\mathrm{p}(#1)}}% could change to Prob or Pr +\newcommand{\pr}[1]{\ensuremath{\mathrm{p}(#1)}} \newcommand{\expect}[1]{\left<#1\right>} \newcommand{\normal}[2]{\mathcal{N} (#1, #2)} -\newcommand{\gvn}{\mid}% could use | or \vert +\newcommand{\gvn}{\mid} \newcommand{\integral}[2]{\ensuremath{\int #1 \mathrm{d} #2}} \newcommand{\sz}{spec-$z$} \newcommand{\Sz}{Spec-$z$} \newcommand{\pz}{photo-$z$} \newcommand{\Pz}{Photo-$z$} -\newcommand{\zpdf}{\pz\ PDF}% could change to posterior -\newcommand{\Zpdf}{\Pz\ PDF}% could change to posterior -\newcommand{\pzpdf}{\pz\ posterior PDF}% could change to implicit posterior -\newcommand{\Pzpdf}{\Pz\ posterior PDF}% could change to implicit posterior +\newcommand{\zpdf}{\pz\ PDF} +\newcommand{\Zpdf}{\Pz\ PDF} +\newcommand{\pzpdf}{\pz\ posterior PDF} +\newcommand{\Pzpdf}{\Pz\ posterior PDF} \newcommand{\pzip}{\pz\ implicit posterior} \newcommand{\nz}{$n(z)$} \newcommand{\Nz}{$N(z)$} \newcommand{\stack}{$\hat{n}(z)$} \newcommand{\ntot}{\ensuremath{N_{\mathrm{tot}}}} -\newcommand{\bvec}[1]{\ensuremath{\boldsymbol{#1}}}% could change to \vec +\newcommand{\bvec}[1]{\ensuremath{\boldsymbol{#1}}} \newcommand{\ndphi}{\bvec{\phi}} -\newcommand{\mmle}{marginalized maximum likelihood estimate}% really marginalized maximum a posteriori estimate, may still change this +\newcommand{\mmle}{marginalized maximum posterior estimate} \begin{document} %\linenumbers \title{How to obtain the redshift distribution from probabilistic redshift estimates} -\author{Alex I. Malz\altaffilmark{1,2}} +\author{Alex I. Malz\altaffilmark{1}} \author{David W. Hogg\altaffilmark{2,3,4,5}} -\email{aimalz@nyu.edu} +\email{aimalz@astro.ruhr-uni-bochum.de} -\altaffiltext{1}{German Centre of Cosmological Lensing, Ruhr-Universit\"{a}t, Universit\"{a}tsstra{\ss}e 150, 44801 Bochum, Germany} +\altaffiltext{1}{Ruhr-University Bochum, German Centre for Cosmological Lensing, Universit\"{a}tsstra{\ss}e 150, 44801 Bochum, Germany} \altaffiltext{2}{Center for Cosmology and Particle Physics, Department of Physics, New York University, 726 Broadway, 9th floor, New York, NY 10003, USA} \altaffiltext{3}{Simons Center for Computational Astrophysics, 162 Fifth Avenue, 7th floor, New York, NY 10010, USA} \altaffiltext{4}{Center for Data Science, New York University, 60 Fifth Avenue, 7th floor, New York, NY 10003, USA} @@ -99,7 +89,7 @@ We present the \texttt{chippr} prototype code %and use it to forecast constraints in the space of cosmological parameters , noting that the mathematically justifiable approach incurs computational expense. -The \textsc{CHIPPR} approach is applicable to any one-point statistic of any random variable, provided the prior probability density used to produce the posteriors is explicitly known; +The \textsc{CHIPPR} approach is applicable to any one-point statistic of any random variable, provided the prior probability density used to produce the posteriors is explicitly known; if the prior is implicit, as may be the case for popular photo-$z$ techniques, then the resulting posterior PDFs cannot be used for scientific inference. We therefore recommend that the photo-$z$ community focus on developing methodologies that enable the recovery of photo-$z$ likelihoods with support over all redshifts, either directly or via a known prior probability density. \end{abstract} @@ -108,28 +98,28 @@ \maketitle -%\aim{TODO: Add test case panel labels to all figs of \chippr\ sims/outputs. Take higher-z samples in mock \pzpdf\ plots relative to outlier populations. Include $\bar{z}$ in results plots.} +%\aim{TODO: Add test case panel labels to all figs of \chippr\ sims/outputs. Take higher-z samples in mock \pzpdf\ plots relative to outlier populations. Include $\bar{z}$ in results plots and cite DES/KiDS papers that use n(z) bias as primary metric.} \section{Introduction} \label{sec:intro} % what are photo-zs? -Photometric redshift (\pz) estimation has been a staple of studies of galaxy evolution, large-scale structure, and cosmology since its conception half a century ago \citep{baum_photoelectric_1962}. -An extremely coarse spectrum in the form of photometry in a handful of broadband filters can be an effective substitute for the time- and photon-intensive process of obtaining a spectroscopic redshift (\sz), a procedure that may only be applied to relatively bright galaxies. +Photometric redshift (\pz) estimation has been a staple of studies of galaxy evolution, large-scale structure, and cosmology since its conception half a century ago \citep{baum_photoelectric_1962}. +An extremely coarse spectrum in the form of photometry in a handful of broadband filters can be an effective substitute for the time- and photon-intensive process of obtaining a spectroscopic redshift (\sz), a procedure that may only be applied to relatively bright galaxies. Once the photometric colors are calibrated against either a library of spectral energy distribution (SED) templates or a data set of spectra for galaxies with known redshifts, a correspondence between photometric colors and redshifts may be constructed, forming a trustworthy basis for \pz\ estimation or testing. % why do we need photo-zs? -Calculations of correlation functions of cosmic shears and galaxy positions that constrain the cosmological parameters require large numbers of high-confidence redshifts of surveyed galaxies. -Many more \pz s may be obtained in the time it would take to observe a smaller number of \sz s, and \pz s may be measured for galaxies too dim for accurate \sz\ confirmation, permitting the compilation of large catalogs of galaxies spanning a broad range of redshifts and luminosities. -\Pz s have thus enabled the era of precision cosmology, heralded by weak gravitational lensing tomography and baryon acoustic oscillation peak measurements. +Calculations of correlation functions of cosmic shears and galaxy positions that constrain the cosmological parameters require large numbers of high-confidence redshifts of surveyed galaxies. +Many more \pz s may be obtained in the time it would take to observe a smaller number of \sz s, and \pz s may be measured for galaxies too dim for accurate \sz\ confirmation, permitting the compilation of large catalogs of galaxies spanning a broad range of redshifts and luminosities. +\Pz s have thus enabled the era of precision cosmology, heralded by weak gravitational lensing tomography and baryon acoustic oscillation peak measurements. % what's wrong with photo-zs? -However, \pz s are susceptible to inaccuracy and imprecision in the form of their inherent noisiness resulting from the coarseness of photometric filters, catastrophic errors in which galaxies of one SED at one redshift are mistaken for galaxies of another SED at a different redshift, and systematics introduced by observational techniques, data reduction processes, and training or template set limitations. +However, \pz s are susceptible to inaccuracy and imprecision in the form of their inherent noisiness resulting from the decreasing signal-to-noise ratio with increasing redshift, coarseness of photometric filters, catastrophic errors in which galaxies of one SED at one redshift are mistaken for galaxies of another SED at a different redshift, and systematics introduced by observational techniques, data reduction processes, and training or template set limitations. Figure~\ref{fig:pedagogical_scatter} is an adaptation of the ubiquitous plots of \pz\ vs. \sz\ illustrating the assumptions underlying \pz\ estimation in general, that \sz s are a good approximation to true redshifts and \pz s represent special non-linear projections of observed photometry to a scalar variable that approximates the true redshift. \begin{figure} \begin{center} - \includegraphics[width=0.45\textwidth]{figures/jain05.png} + \includegraphics[width=0.4\textwidth]{jain15.png} \caption{ A generic probability space (darker in areas of higher probability density) of true redshift ($x$-axis) and a nonlinear projection of photometric data ($y$-axis), with vertical cuts and marginals (orange) indicating the construction of likelihoods and horizontal cuts and marginals (blue) indicating the construction of posteriors, with a theoretically perfect \pz\ estimate on the diagonal (yellow) for reference. The data points were extracted using WebPlotDigitizer \citep{rohatgi_webplotdigitizer_2019} from a \sz\ vs. \pz\ plot in \citet{jain_whole_2015}. @@ -149,12 +139,12 @@ \section{Introduction} Once propagated through the calculations of correlation functions of cosmic shear and galaxy positions, \pz\ errors are a dominant contributor to the total uncertainties reported on cosmological parameters \citep{abruzzo_impact_2019}. As progress has been made on the influence of other sources of systematic error, the uncertainties associated with \pz s have come to dominate the error budget of cosmological parameter estimates made by current surveys such as \des\ \citep{hoyle_dark_2018}, \project{HSC} \citep{tanaka_photometric_2018}, and \project{KiDS} \citep{hildebrandt_kids-450:_2017}. Based on the goals of a photometric galaxy survey, limits can be placed on the tolerance to these effects. -For example, the Science Requirements Document \citep{mandelbaum_weak_2017} states \lsst's requirements for the main cosmological sample, reproduced in Table~\ref{tab:lsstsrd}. +For example, the Science Requirements Document \citep{mandelbaum_weak_2017} states \lsst's requirements \new{for}{on \pz\ error tolerances for constraining cosmology with} the main cosmological sample, reproduced in Table~\ref{tab:lsstsrd}. \begin{table} \begin{center} \caption{\Pz\ requirements for \lsst\ cosmology\\ - \citep{mandelbaum_weak_2017}.} + \citep{mandelbaum_weak_2017}\new{.}{}} \begin{tabular}{ll} Number of galaxies & $\approx 10^{7}$\\ Root-mean-square error & $< 0.02 (1 + z)$\\ @@ -170,19 +160,22 @@ \section{Introduction} The presence of galaxies whose SEDs are not represented by the template library tends to lead to catastrophic outliers distributed like the horizontally oriented population of \Fig{fig:pedagogical_scatter}. For data-driven approaches, training sets that are incomplete in redshift coverage tend to result in catastrophic outliers like the vertically oriented population of \Fig{fig:pedagogical_scatter}. The approaches of using a training set versus a template library are related to one another by \citet{budavari_unified_2009}. -Sophisticated Bayesian techniques and machine learning methods have been employed to improve precision \citep{carliles_random_2010} and accuracy \citep{sadeh_annz2:_2016}, while other advances have focused on identifying and removing catastrophic outliers when using \pz s for inference \citep{gorecki_new_2014}. +Sophisticated Bayesian techniques and machine learning methods have been employed to improve precision \citep{carliles_random_2010} and accuracy \citep{sadeh_annz2:_2016}, while other advances have focused on identifying and removing catastrophic outliers when using \pz s for inference \citep{gorecki_new_2014}. % PDFs are a better way to improve photo-zs The probability density function (PDF) in redshift space for each galaxy, commonly written as $\pr{z}$, is an alternative to the MLE (with or without presumed Gaussian error bars) \citep{koo_photometric_1999}. This option is favorable because it contains more potentially useful information about the uncertainty on each galaxy's redshift, incorporating our understanding of precision, accuracy, and systematic error. -However, denoting \zpdf s as ``$\pr{z}$'' is an abuse of notation, as it does not adequately convey what information is being used to constrain the redshift $z$; \zpdf s are \textit{posterior} PDFs, conditioned on the photometric data and prior knowledge. +However, denoting \zpdf s as ``$\pr{z}$'' is an abuse of notation, as it does not adequately convey what information is being used to constrain the redshift $z$; +\zpdf s are \textit{posterior} PDFs, conditioned on the photometric data and prior knowledge. In terms of \Fig{fig:pedagogical_scatter}, \zpdf s are horizontal cuts, probabilities of redshift conditioned on a specific value of data, i.e. posteriors $\pr{z \gvn \data}$, which constrain redshifts, whereas vertical cuts through this space are probabilities of data conditioned on a specific redshift, i.e. likelihoods $\pr{\data \gvn z}$, from which photometric data is actually drawn. % photo-z PDFs are established -\Pzpdf s have been produced by completed surveys \citep{hildebrandt_cfhtlens:_2012, sheldon_photometric_2012} and will be produced by ongoing and upcoming surveys \citep{abell_lsst_2009, carrasco_kind_exhausting_2014, bonnett_redshift_2016, masters_mapping_2015}. -\Pzpdf s are not without their own shortcomings, however, including the resources necessary to calculate and record them for large galaxy surveys \citep{carrasco_kind_sparse_2014, malz_approximating_2018} and the divergent results of each method used to derive them \citep{hildebrandt_phat:_2010, dahlen_critical_2013, sanchez_clustering_2013, bonnett_redshift_2016, tanaka_photometric_2018}. +\Pzpdf s have been produced by completed surveys \citep{hildebrandt_cfhtlens:_2012, sheldon_photometric_2012} and will be produced by ongoing and upcoming surveys \citep{abell_lsst_2009, carrasco_kind_exhausting_2014, bonnett_redshift_2016, masters_mapping_2015}. +\Pzpdf s are not without their own shortcomings, however, including the resources necessary to calculate and record them for large galaxy surveys \citep{carrasco_kind_sparse_2014, malz_approximating_2018} and the divergent results of each method used to derive them \citep{hildebrandt_phat:_2010, dahlen_critical_2013, sanchez_clustering_2013, bonnett_redshift_2016, tanaka_photometric_2018}. Though the matter is outside the scope of this paper, reviews of various methods have been presented in the literature \citep{sheldon_photometric_2012, ball_robust_2008, carrasco_kind_tpz:_2013, carrasco_kind_exhausting_2014, schmidt_evaluation_2020}. -The most concerning weakness of \pzpdf s, however, is their usage in the literature, which is at best inconsistent and at worst incorrect. +The most concerning weakness of \pzpdf s, however, is their usage in the literature, which is at best inconsistent and at worst incorrect. + +{\color{white}{\Fig{fig:pgm}}} % photo-z PDFs are most often reduced to point estimates Though their potential to improve estimates of physical parameters is tremendous, \pzpdf s have been applied only to a limited extent, most often by reduction to familiar point estimates. @@ -192,7 +185,7 @@ \section{Introduction} % photo-z PDFs may also be used to define cuts Regardless of how it is done, any procedure that reduces \pzpdf s to point estimates discards valuable information about the uncertainty on redshift. -\Pzpdf s have also been used to form selection criteria of samples from galaxy surveys without propagation through the calculations of physical parameters \citep{van_breukelen_reliable_2009, viironen_high_2015}. +\Pzpdf s have also been used to form selection criteria of samples from galaxy surveys without propagation through the calculations of physical parameters \citep{van_breukelen_reliable_2009, viironen_high_2015}. Probability cuts on Bayesian quantities are not uncommon \citep{leung_bayesian_2017, dipompeo_quasar_2015}, but that procedure does not fully take advantage of all information contained in a probability distribution for parameter inference. The most prevalent application of \pzpdf s that preserves their information content is the estimation of the \textit{redshift distribution function \Nz}, or, interchangably, its normalized cousin the \textit{redshift density function \nz}. @@ -201,7 +194,7 @@ \section{Introduction} As a key input to the traditional calculation of the power spectra of weak gravitational lensing and large-scale structure, the accuracy and precision to which \Nz\ is estimated can strongly impact our constraints on the parameters of cosmological models \citep{bonnett_using_2015, masters_mapping_2015, viironen_high_2015, asorey_galaxy_2016, bonnett_redshift_2016, yang_calibrating_2018}, so it is unsurprising that this last application dominates the canonical bias requirement of Table~\ref{tab:lsstsrd}. %\aim{TODO: Say why \Nz\ matters for cosmology, what precision is needed for \Nz\ for future weak lensing surveys, motivate how well we need to know \Nz.} Even with \pz s adhering to the \lsst\ requirements of \Tab{tab:lsstsrd}, the degree to which constraints on the cosmological parameters can advance is limited by the accuracy and precision to which \nz\ is known \citep{abruzzo_impact_2019}. -% Say what precision the mass function is needed (in, say cluster studies) for precision cosmology. +%\aim{TODO: Say what precision the mass function is needed (in, say cluster studies) for precision cosmology.} % how people get n(z) from photo-z PDFs Though it is traditional to estimate \nz\ from \pz\ point estimates \citep{abruzzo_impact_2019}, it has become more common to use \pzpdf s directly to calculate the conceptually simple but mathematically inconsistent \citep{hogg_data_2012} \textit{stacked estimator} $\hat{n}(z)$ of the redshift density function \citep{lima_estimating_2008} @@ -210,16 +203,17 @@ \section{Introduction} \hat{n}(z) &= \frac{1}{J} \sum_{j = 0}^{J} \pr{z}_{j} \end{align} for a sample of $J$ galaxies $j$, or, equivalently, the redshift distribution function $\hat{N}(z) = J \hat{n}(z)$, by effectively averaging the \pzpdf s. -This summation procedure has been used extensively in cosmological analyses with photometric galaxy samples \citep{mandelbaum_precision_2008, benjamin_cfhtlens_2013, kelly_weighing_2014}. +\new{T}{Despite its problems \citep{malz_how_2021}, t}his summation procedure has been used extensively in cosmological analyses with photometric galaxy samples \citep{mandelbaum_precision_2008, benjamin_cfhtlens_2013, kelly_weighing_2014}. %\aim{TODO: Continue adding equations/citations for the use of \Nz\ in cosmology.} % what this paper is about -Despite the growing prevalence of \pzpdf\ production, no implementation of inference using \pzpdf s has yet been presented with a mathematically consistent methodology. -This paper challenges the logically invalid yet pervasive analysis procedure of stacking \pzpdf s by presenting and validating a hierarchical Bayesian technique for the use of \pzpdf s\ in the inference of \nz, yielding a method applicable to arbitrary one-point statistics relevant to cosmology, large-scale structure, and galaxy evolution; future work will extend this methodology to higher-order statistics. +Despite the growing prevalence of \pzpdf\ production, no implementation of inference using \pzpdf s has yet been presented with a mathematically consistent methodology. +This paper challenges the logically invalid yet pervasive analysis procedure of stacking \pzpdf s by presenting and validating a hierarchical Bayesian technique for the use of \pzpdf s\ in the inference of \nz, yielding a method applicable to arbitrary one-point statistics relevant to cosmology, large-scale structure, and galaxy evolution; +future work will extend this methodology to higher-order statistics. We aim to develop a clear methodology guiding the use of \pzpdf s in inference so they may be utilized effectively by the cosmology community. -Though others have approached the problem before \citep{leistedt_hierarchical_2016, leistedt_hierarchical_2018}, the method presented here differs in that it makes use of any existing catalog of \pzpdf s, rather than requiring a simultaneous derivation of the \pzpdf s and the redshift distribution, making it preferable to ongoing surveys for which there may be inertia preventing a complete restructuring of the analysis pipeline. +Though others have approached the problem before \citep{leistedt_hierarchical_2016, leistedt_hierarchical_2019}, the method presented here differs in that it makes use of any existing catalog of \pzpdf s, rather than requiring a simultaneous derivation of the \pzpdf s and the redshift distribution, making it preferable to ongoing surveys for which there may be inertia preventing a complete restructuring of the analysis pipeline. -In Section~\ref{sec:meth}, we present the \Chippr\ model for characterizing the full posterior probability landscape of \Nz\ using \pzpdf s. +In Section~\ref{sec:meth}, we present the \Chippr\ model for characterizing the full posterior probability landscape of \Nz\ using \pzpdf s. In Section~\ref{sec:application}, we present the \chippr\ implementation of the \Chippr\ model and the experimental set up by which we validate it, including the forward modeling of mock \pzpdf s. In Section~\ref{sec:alldata}, we present a number of informative test cases and compare the results of \chippr\ with alternative approaches. In Section~\ref{sec:results}, we stress-test the \Chippr\ model under nontraditional conditions. @@ -229,18 +223,16 @@ \section{Introduction} \section{Model} \label{sec:meth} -Consider a survey of $J$ galaxies $j$, each with photometric data $\data_{j}$; -thus the entire survey over some solid angle produces the ensemble of photometric magnitudes (or colors) and their associated observational errors $\{\data_{j}\}$. -Each galaxy $j$ has a redshift parameter $z_{j}$ that we would like to learn. +Consider a survey of $J$ galaxies $j$, each with photometric data $\data_{j}$; +thus the entire survey over some solid angle produces the ensemble of photometric magnitudes (or colors) and their associated observational errors $\{\data_{j}\}$. +Each galaxy $j$ has a redshift parameter $z_{j}$ that we would like to learn. The distribution of the ensemble of redshift parameters $\{z_{j}\}$ may be described by the hyperparameters defining the redshift distribution function \nz\ that we would like to quantify. -The redshift distribution function \nz\ is the number of galaxies per unit redshift, effectively defining the evolution in the number of galaxies convolved with the selection function of the sample \citep{menard_clustering-based_2013}. +The redshift distribution function \nz\ is the number of galaxies per unit redshift, effectively defining the evolution in the number of galaxies convolved with the selection function of the sample \citep{menard_clustering-based_2013}. -In \Sect{sec:forward}, we establish a forward model encapsulating the causal relationship between \nz\ and photometry $\data$. +In \Sect{sec:forward}, we establish a forward model encapsulating the causal relationship between \nz\ and photometry $\data$. In \Sect{sec:prob}, we present the directed acyclic graph of this probabilistic generative model and interpret the corresponding mathematical expression, whose full derivation may be found in the Appendix. In \Sect{sec:limitations}, we summarize the necessary assumptions of the model. -%\aim{TODO: check for consistent notation of \pzip s vs. \pzpdf s.} - \subsection{Forward Model} \label{sec:forward} @@ -251,7 +243,8 @@ \subsection{Forward Model} \int_{-\infty}^{\infty}\ n(z)\ dz\ \equiv\ \frac{1}{J}\ \int_{-\infty}^{\infty}\ \sum_{j=1}^{J}\ \delta(z_{j},\ z)\ dz = 1 \end{equation} of finding a galaxy $j$ in a catalog of $J$ galaxies having a redshift $z$. -We believe that galaxy redshifts are indeed sampled, or drawn, from \nz, making it a probability density over redshift; this fact can also be confirmed by dimensional analysis of \Eq{eqn:nz}, as suggested in \citet{hogg_data_2012}. +We believe that galaxy redshifts are indeed sampled, or drawn, from \nz, making it a probability density over redshift; +this fact can also be confirmed by dimensional analysis of \Eq{eqn:nz}, as suggested in \citet{hogg_data_2012}. We may without loss of generality impose a parameterization \begin{equation} @@ -260,7 +253,7 @@ \subsection{Forward Model} \end{equation} in terms of some parameter vector $\ndphi$. At this point, the parameter vector is quite general and may represent coefficients in a high-order polynomial as a function of redshift, a set of means and variances defining Gaussians that sum to the desired distribution, a set of histogram heights that describe a binned version of the redshift distribution function, etc. -Upon doing so, we may rewrite \Eq{eqn:fz} as +Upon doing so, we may rewrite \Eq{eqn:fz} as \begin{equation} \label{eqn:pz} z_{j}\ \sim\ \pr{z \gvn \ndphi}\ \equiv\ f(z; \ndphi), @@ -275,9 +268,9 @@ \subsection{Forward Model} n(z)\ =\ \integral{\pr{z, \data}}{\data} \end{equation} over the dimension of data in that joint probability space. -Note that galaxies may have different observational data despite sharing the same redshift, and that galaxies at different redshifts may have identical photometry; +Note that galaxies may have different observational data despite sharing the same redshift, and that galaxies at different redshifts may have identical photometry; the space $\pr{z, \data}$ need not be one-to-one. -We assume a stronger version of statistical independence here, that draws $(z_{j}, \data_{j})$ are independent of draws $(z_{j'}, \data_{j'})$ in this space; +We assume a stronger version of statistical independence here, that draws $(z_{j}, \data_{j})$ are independent of draws $(z_{j'}, \data_{j'})$ in this space; the data and redshift of each galaxy are independent of those of other galaxies. However, this problem has additional causal structure that we can acknowledge. @@ -309,10 +302,10 @@ \subsection{Probabilistic Model} \begin{figure} \begin{center} - \includegraphics[height=0.25\textheight]{figures/chippr/pgm.png} + \includegraphics[height=0.25\textheight]{pgm.png} \caption{The directed acyclic graph of the CHIPPR model, where circles indicate random variables and arrows indicate causal relationships. - The redshift distribution \nz\ parameterized by hyperparameters $\ndphi$ exists independent of the survey of $J$ galaxies, indicated as a box. - The redshifts $\{z_{j}\}$ of all galaxies in the survey are latent variables independently drawn from the redshift distribution, which is a function of $\ndphi$. + The redshift distribution \nz\ parameterized by hyperparameters $\ndphi$ exists independent of the survey of $J$ galaxies, indicated as a box. + The redshifts $\{z_{j}\}$ of all galaxies in the survey are latent variables independently drawn from the redshift distribution, which is a function of $\ndphi$. The photometric data $\data_{j}$ for each galaxy is drawn from a function of its redshift $z_{j}$ and observed, indicated by a shaded circle. } \label{fig:pgm} @@ -321,29 +314,24 @@ \subsection{Probabilistic Model} The problem facing cosmologists is to determine the true value of $\ndphi$ from observing the photometry $\{\data_{j}\}$ of a large sample of $J$ galaxies $j$. To self-consistently propagate the uncertainty in the inference of redshift, however, it is more appropriate to estimate the posterior $\pr{\ndphi \gvn \{\data_{j}\}}$ over all possible values of $\ndphi$ conditioned on all the observed data $\{\data_{j}\}$ available in a generic catalog. -In order to use the DAG of \Fig{fig:pgm} to derive an expression for $\pr{\ndphi \gvn \{\data_{j}\}}$ in terms of \pzpdf s, we must introduce two more concepts, confusingly named the \textit{implicit prior} and the \textit{prior probability density} (\textit{prior PDF}), elaborated upon below. +In order to use the DAG of \Fig{fig:pgm} to derive an expression for $\pr{\ndphi \gvn \{\data_{j}\}}$ in terms of \pzpdf s, we must introduce two more concepts, confusingly named the \new{\textit{implicit prior} and the \textit{prior probability density} (\textit{prior PDF})}{implicit prior and the (hyper)prior probability density function}, elaborated upon below. When we constrain the redshift of a galaxy using its observed photometric data $\data_{j}$, we are effectively estimating a posterior $\pr{z \gvn \data_{j}}$, the probability of an unknown quantity conditioned on the quantity we have in hand, i.e the photometric data. This posterior is effectively a marginalization with respect to redshift at a given value of $\data = \data_{j}$ of the \textit{empirical frequency distribution} $\pr{z, \data \gvn \ndphi^{\dagger}}$, the joint probability density corresponding to the true redshift distribution parameterized by $\ndphi^{\dagger}$, which exists in nature but need not be known. -%\que{Propagate new $\pr{z, \data \gvn \ndphi^{\dagger}}$ notation through appendix?} As the hyperparameters $\ndphi^{\dagger}$ of the true redshift distribution are in general unknown, the investigator seeking to estimate a posterior $\pr{z \gvn \data_{j}}$ must have a model $\phi^{*}$ for the general relationship between redshifts and photometry, whether empirical, as is the case for machine learning \pzpdf\ methods, or analytic, as is the case for template-based \pzpdf\ methods. If we were to marginalize over the photometry in $\pr{\data, z}$, we would obtain a one-dimensional PDF $\pr{z \gvn \ndphi^{*}}$ over redshift, which can by definition be parameterized by the same functional form as \nz, for some $\ndphi^{*}$ specific to the estimation procedure that may or may not bear any relation to the hyperparameters $\ndphi^{\dagger}$ of the true \nz. -Rather, $\ndphi^{*}$ is a consequence of the generative model for how photometry results from redshift, including the influence of intrinsic galaxy spectra and instrumental effects. +Rather, $\ndphi^{*}$ is a consequence of the generative model for how photometry results from redshift, including the influence of intrinsic galaxy spectra and instrumental effects. -We call $\pr{z \gvn \ndphi^{*}}$ the \textit{implicit prior}, as it is rarely explicitly known nor chosen by the researcher\footnote{For template-based methods, the implicit prior is often an explicitly known input to the algorithm, engineered as an initial guess for the true $\ndphi$, with an aim for a realistic choice guided by an earlier spectroscopic survey. +We call $\pr{z \gvn \ndphi^{*}}$ the \textit{implicit prior}, as it is rarely explicitly known nor chosen by the researcher\footnote{For template-based methods, the implicit prior is often an explicitly known input to the algorithm, engineered as an initial guess for the true $\ndphi$, with an aim for a realistic choice guided by an earlier spectroscopic survey. (See \citet{benitez_bayesian_2000} for more detail.) It may thus be more appropriate to call it an \textit{interim prior}, but we will use the former term throughout this paper for generality.} Because the implicit prior is unavoidable and almost inherently not uninformative, the \pzpdf s reported by any method must be \textit{implicit posteriors} ${\pr{z \gvn \data, \ndphi^{*}}}$ weighted by the implicit prior. -%Posteriors differ from likelihoods by way of a prior distribution, so we cannot simply assume that the available data products are \pz\ posteriors $\pr{z \gvn \data_{j}}$. -%Rather, we have a catalog of implicit-prior weighted \pz\ posteriors $\pr{z \gvn \data_{j}, \ndphi^{*}}$. -%There must have been some interim prior probability distribution $p(z|\vec{\theta}^{*})$ defined in terms of the interim prior parameter values (hereafter the interim prior) $\vec{\theta}^{*}$ explicitly chosen or implicitly made to perform the calculation of the probabilistic photo-$z$s. -%If it is implicit, it may not be representable in the parametrization we have chosen, and furthermore it may not be known at all; a method that produces interim photo-$z$ posteriors of this kind is not suitable for inference. -%However, so long as the implicit prior is known, hierarchical inference is possible. - -The prior probability density $\pr{\ndphi}$ is a more familiar concept in astronomy; to progress, we will have to choose a prior probability density over all possible values of the hyperparameters $\ndphi$. -This prior need not be excessively proscriptive; for example, it may be chosen to enforce smoothness at physically motivated scales in redshift without imposing any particular region as over- or under-dense. +The \textit{prior probability density} $\pr{\ndphi}$ is a more familiar concept in astronomy; +to progress, we will have to choose a prior probability density over all possible values of the hyperparameters $\ndphi$. +This prior need not be excessively proscriptive; +for example, it may be chosen to enforce smoothness at physically motivated scales in redshift without imposing any particular region as over- or under-dense. With inputs of the \pzip\ catalog $\{\pr{z \gvn \data, \ndphi^{*}}\}$, the implicit prior $\pr{z \gvn \ndphi^{*}}$, and the prior PDF $\pr{\ndphi}$, we thus aim to obtain the posterior probability $\pr{\ndphi \gvn \{\data_{j}\}}$ of the redshift density function given all the photometric data. By performing the derivation of the Appendix, we arrive at the desired expression @@ -351,7 +339,7 @@ \subsection{Probabilistic Model} \label{eqn:fullpost} \pr{\ndphi \gvn \{\data_{j}\}} \propto \pr{\ndphi} \integral{\prod_{j=1}^{J} \frac{\pr{z \gvn \data_{j}, \ndphi^{*}} \pr{z \gvn \ndphi}}{\pr{z \gvn \ndphi^{*}}}}{z}, \end{equation} -which is the very heart of \Chippr, also given as \Eq{eqn:final}. +which is the very heart of \Chippr, also given as \Eq{eqn:final} and in the form of interpretive dance \citep{malz_probabilistic_2019}. This in effect replaces the implicit prior with the sampled model hyperparameters, thereby converting the \pzip s into likelihoods in order to obtain unbiased posteriors. \subsection{Model Limitations} @@ -360,49 +348,49 @@ \subsection{Model Limitations} Finally, we explicitly review the assumptions made by this approach, which are as follows: \begin{enumerate} \item Photometric measurements of galaxies are statistically independent Poisson draws from the set of all galaxies such that \Eq{eqn:indiedat} and \Eq{eqn:indie} hold. - \item We take the reported \pzip s to be accurate, free of model misspecification; + \item We take the reported \pzip s to be accurate, free of model misspecification; draws thereof must not be inconsistent with the distribution of photometry and redshifts. Furthermore, we must be given the implicit prior $\ndphi^{*}$ used to produce the \pzip s. \item We must assume a hyperprior distribution $\pr{\ndphi}$ constraining the underlying probability distribution of the hyperparameters, which is informed by our prior beliefs about the true redshift distribution function. \end{enumerate} -These assumptions have known limitations. -First, the photometric data are not a set of independent measurements; +These assumptions have known limitations. +First, the photometric data are not a set of independent measurements; the data are correlated not only by the conditions of the experiment under which they were observed (instrument and observing conditions) but also by redshift covariances resulting from physical processes governing underlying galaxy spectra and their relation to the redshift distribution function. -Second, the reported \pzip s may not be trustworthy; -there is not yet agreement on the best technique to obtain \pzpdf s, and the implicit prior may not be appropriate or even known to us as consumers of \pzip s. +Second, the reported \pzip s may not be trustworthy; +there is not yet agreement on the best technique to obtain \pzpdf s, and the implicit prior may not be appropriate or even known to us as consumers of \pzip s. Third, the hyperprior may be quite arbitrary and poorly motivated if the underlying physics is complex, and it can only be appropriate if our prior beliefs about \nz\ are accurate. -Furthermore, in Section~\ref{sec:prob}, we have made an assumption of \textit{support}, meaning the model $\pr{z, \data \gvn \ndphi}$ has mutual coverage with the parameter values that real galaxies can take. -In other words, any probability distribution over the $(z, \data)$ space must be nonzero where real galaxies can exist. +Furthermore, in Section~\ref{sec:prob}, we have made an assumption of \textit{support}, meaning the model $\pr{z, \data \gvn \ndphi}$ has mutual coverage with the parameter values that real galaxies can take. +In other words, any probability distribution over the $(z, \data)$ space must be nonzero where real galaxies can exist. Additionally, the hyperprior $\pr{\ndphi}$ must be nonzero at the hyperparameters $\ndphi^{\dagger}$ of the true redshift density function \nz. This assumption cannot be violated under the experimental design of Section~\ref{sec:forward}, but it is not generically guaranteed when performing inference on real data; thus the chosen $\pr{z, \data \gvn \ndphi^{*}}$ and $\pr{\ndphi}$ must be sufficiently general as to not rule out plausible areas of parameter space. +%{\color{white}{\Fig{fig:flowchart}, \Fig{fig:mega_scatter}}} \section{Methods \& Data} \label{sec:application} Here we describe the method by which we demonstrate the \Chippr\ model. In \Sect{sec:exp}, we outline the implementation of the \chippr\ code. -%In \Sect{sec:sheldon}, we introduce the alternative \nz\ estimators against which \Chippr\ is compared. -%In \Sect{sec:diag}, we present the quantitative metrics by which the \nz\ estimators are compared. In \Sect{sec:mock}, we outline the procedure for emulating mock \pzip s. \subsection{Implementation} \label{sec:exp} -We implement the \Chippr\ model in code in order to perform tests of its validity and to compare its performance to that of traditional alternatives. +We implement the \Chippr\ model in \new{}{a prototype} code in order to perform tests of its validity and to compare its performance to that of traditional alternatives. In \Sect{sec:mcmc}, we describe the publicly available \chippr\ library. In \Sect{sec:sheldon}, we introduce the alternative approaches evaluated for comparison with \Chippr. In \Sect{sec:diag}, we describe the diagnostic criteria by which we assess estimators of \nz. -%In \Sect{app:acorr}, we outline how \chippr\ can be used to sample the full log-posterior distribution $\ln[\pr{\ndphi \gvn \{\data_{j}\}}]$. \subsubsection{Code} \label{sec:mcmc} -\chippr\ is a \python\ 2 library\footnote{\url{https://github.com/aimalz/chippr}} that includes an implementation of the \Chippr\ model as well as an extensive suite of tools for comparing \Chippr\ to other approaches. +\chippr\ is a \new{}{public} \python\ 2 library\footnote{\url{https://github.com/aimalz/chippr}} that includes an implementation of the \Chippr\ model as well as an extensive suite of tools for comparing \Chippr\ to other approaches. +\new{}{For the sake of reproducibility, this paper presents the prototype-level code used for the experiments discussed herein. +However, a modern, flexible implementation appropriate for at-scale application to future data sets is currently under development\footnote{A working prototype of a \python 3-based implementation built on the \texttt{qp} back-end presented in \citet{malz_approximating_2018} may be found at \url{https://github.com/LSSTDESC/qit}, which is anticipated to evolve in sophistication over time.}.} -Though there are plans for future expansion to more flexible parameterizations, the current version of \chippr\ uses a log-space piecewise constant parameterization +\new{Though there are plans for future expansion to more flexible parameterizations, the version of }{}\chippr\ uses a log-space piecewise constant parameterization \begin{equation} \label{eqn:logstepfunc} f(z; \ndphi) = \exp[\phi^{k}]\ \mathrm{if}\ z^{k} < z < z^{k+1} @@ -416,29 +404,24 @@ \subsubsection{Code} Thus each $\pr{z \gvn \data_{j}} = f(z; \ndphi_{j})$ has parameters $\ndphi_{j}$ that are defined in the same basis as those of \nz. To infer the full log-posterior distribution $\ln[\pr{\ndphi \gvn \{\data_{j}\}}]$, one must provide a plaintext file with $K+1$ redshift bin endpoints $\{z_{k}\}$, the parameters $\ndphi^{*}$ of the implicit log-prior, and the parameters $\{\ndphi_{j}\}$ of the log-posteriors $\ln[\pr{z \gvn \data_{j}, \ndphi^{*})}$. -The \emcee \citep{foreman-mackey_emcee_2013} implementation of ensemble sampling is used to sample the full log-posterior of \Eq{eqn:final}. +The \emcee \citep{foreman-mackey_emcee_2013} implementation of \new{}{MCMC} ensemble sampling is used to sample the full log-posterior of \Eq{eqn:final}. +\new{}{As a consequence of the Central Limit Theorem applied to a Markov Chain, the sample distribution satisfies a statistical guarantee of normality \citep{brooks_handbook_2011}, meaning the output samples $\{\ndphi_{i}\}$ of \chippr\ naturally define a notion of standard deviation interpretable as conventional error bars\footnote{Additionally, the proofs are tractable by hand in the special cases of simplified \pzpdf s introduced in \citet{malz_how_2021}.}.} + \chippr\ accepts a configuration file of user-specified parameters, among them the number $W$ of walkers. At each iteration $i$ and for each walker, a proposal distribution $\hat{\ndphi}_{i}$ is drawn from the log-prior distribution and evaluated for acceptance to or rejection from the full log-posterior distribution. -%Two threshold conditions are defined, one designating all previous samples to be ignored as as products of a burn-in phase and another indicating when a sufficient number of post-burn samples have been accepted. -%In this case, the first threshold (described in \Sect{app:acorr}) is defined in terms of sub-runs of $10^{3}$ accepted samples, and the second is defined as an accumulation of $10^{4}$ samples. -%Though previous versions used \texttt{HDF5} for the primary I/O format due to its efficiency for large quantities of data, it was abandoned in favor of \texttt{pickle} in the working release due to the instability of the \python\ implementation of the format on high-performance computing systems. - -The resulting output -% is a set of ordered \texttt{pickle} files % enumerated by $\rho$ each containing the state information of $10^{3}$ accepted samples, with the last ten files including only samples taken after the completion of the burn-in phase -% after each sub-run. -%The state information -includes $\frac{I_{0}}{s}$ accepted samples $\ndphi_{i}$ for a pre-specified chain thinning factor $s$ and their full posterior probabilities $\pr{\ndphi_{i} \gvn \{\data_{j}\}}$, as well as the autocorrelation times and acceptance fractions calculated for each element of $\ndphi$, divided into separate files before and after the completion of the burn-in phase, as defined by the Gelman-Rubin statistic \citep{gelman_inference_1992}. +%In \Sect{app:acorr}, we outline how \chippr\ can be used to sample the full log-posterior distribution $\ln[\pr{\ndphi \gvn \{\data_{j}\}}]$. +The resulting output includes $\frac{I_{0}}{s}$ accepted samples $\ndphi_{i}$ for a pre-specified chain thinning factor $s$ and their full posterior probabilities $\pr{\ndphi_{i} \gvn \{\data_{j}\}}$, as well as the autocorrelation times and acceptance fractions calculated for each element of $\ndphi$, divided into separate files before and after the completion of the burn-in phase, as defined by the Gelman-Rubin statistic $R$ \citep{gelman_inference_1992}. \subsubsection{Alternative approaches for comparison} \label{sec:sheldon} -In this study, we compare the results of \Eq{eqn:fullpost} to those of the two most common approaches to estimating \nz\ from a catalog of \pzip s: +In this study, we compare the results of \Eq{eqn:fullpost} to those of the two most common approaches to estimating \nz\ from a catalog of \pzip s: the distribution $n(z_{\mathrm{max}})$ of the redshifts at maximum posterior probability \begin{equation} \label{eqn:mmap} f^{MMAP}(z; \hat{\ndphi}) = \sum_{j=1}^{J}\ \delta(z, \mathrm{mode}[\pr{z \gvn \data_{j}, \ndphi^{*}}]) \end{equation} -(i.e. the distribution of modes of the \pzip s) and the stacked estimator of \Eq{eqn:stacked}, which can be rewritten as +(i.e. the distribution of modes of the \pzip s) and the stacked estimator of \Eq{eqn:stacked}, which can be rewritten as \begin{equation} \label{eqn:stacked} f^{stack}(z; \hat{\ndphi}) = \sum_{j=1}^{J}\ \pr{z \gvn \data_{j}, \ndphi^{*}} @@ -446,57 +429,41 @@ \subsubsection{Alternative approaches for comparison} in terms of the \pzip s we have. These two approaches have been compared to one another by \citet{hildebrandt_cfhtlens:_2012}, \citet{benjamin_cfhtlens_2013}, and \citet{asorey_galaxy_2016} in the past but not to \Chippr. -Point estimation converts the implicit \pz\ posteriors $\pr{z \gvn \data_{j}, \ndphi^{*}}$ into delta functions with all probability at a single estimated redshift. +Point estimation converts the implicit \pz\ posteriors $\pr{z \gvn \data_{j}, \ndphi^{*}}$ into delta functions with all probability at a single estimated redshift. Some variants of point estimation choose this single redshift to be that of maximum a posteriori probability $\mathrm{mode}[\pr{z \gvn \data_{j}, \ndphi^{*}}]$ or the expected value of redshift $\langle z \rangle = \integral{z \pr{z \gvn \data_{j}, \ndphi^{*}}}{z}$. \citet{tanaka_photometric_2018} directs attention to deriving an optimal point estimate reduction of a \pzpdf, but since the purpose of this paper is to compare against the most established alternative estimators of \nz, its use will be postponed until a future study. Stacking these modified \pzip s leads to the marginalized maximum a posteriori (MMAP) estimator and the marginalized expected value (MExp) estimator, though only the former is included in this study since the latter has fallen out of favor in recent years\footnote{And for good reason! Consider a bimodal \pzpdf; its expected value may very well fall in a region of very low probability, yielding a less probable point estimate than the point at which either peak achieves its maximum.}. -It is worth discussing the relationship between point estimation and stacking. -When the point estimator of redshift is equal to the true redshift, stacking delta function \pzpdf s will indeed lead to an accurate recovery of the true redshift distribution function. -However, stacking is in general applied indiscriminately to broader \pzpdf s and imperfect point estimators of redshift. +It is worth discussing the relationship between point estimation and stacking. +When the point estimator of redshift is equal to the true redshift, stacking delta function \pzpdf s will indeed lead to an accurate recovery of the true redshift distribution function. +However, stacking is in general applied indiscriminately to broader \pzpdf s and imperfect point estimators of redshift. It is for these reasons that alternatives are considered here. -A final estimator of the hyperparameters is the maximum marginalized likelihood estimator (MMLE), the value of $\ndphi$ maximizing the log posterior given by \Eq{eqn:final} using any optimization code. -The MMLE can be obtained in substantially less time than enough samples to characterize the full log-posterior distribution of \nz. -However, the MMLE yields only a point estimate of \nz\ rather than characterizing the full log-posterior on $\ndphi$, and it does not escape the dependence on the choice of hyperprior distribution. +A final estimator of the hyperparameters is the maximum marginalized likelihood estimator (\mmle), the value of $\ndphi$ maximizing the log posterior given by \Eq{eqn:final} using any optimization code. +The \mmle can be obtained in substantially less time than enough samples to characterize the full log-posterior distribution of \nz. +However, the \mmle yields only a point estimate of \nz\ rather than characterizing the full log-posterior on $\ndphi$, and it does not escape the dependence on the choice of hyperprior distribution. Furthermore, derivatives will not in general be available for the full posterior distribution, restricting optimization methods used, and, as is true for any optimization code, there is a risk of numerical instability. -%\begin{equation} -%\label{eq:mmle} -%\ln[p(\{\vec{d}_{j}\}|\vec{\theta})] \propto -\int\ f_{\vec{\theta}}(z)\ -%dz+\sum_{j=1}^{J}\ln\left[\int\ -%\exp\left[\ln[p(z_{j}|\vec{d}_{j},\vec{\theta}^{*})]+\ln[f_{\vec{\theta}}(z)]-\l -%n[f_{\vec{\theta}^{*}}(z)]\right]\ dz\right], -%\end{equation} -%accessible with any optimization code. \subsubsection{Performance metrics} \label{sec:diag} -The results of the computation described in \Sect{sec:exp} are evaluated for accuracy on the basis of some quantitative measures. -Beyond visual inspection of samples, we calculate summary statistics to quantitatively compare different estimators' precision and accuracy. -Since MCMC samples of hyperparameters are Gaussian distributions, we can quantify the breadth of the distribution for each hyperparameter using the standard deviation regardless of whether the true values are known. +The results of the computation described in \Sect{sec:exp} are evaluated for accuracy on the basis of some quantitative measures. +Beyond visual inspection of samples, we calculate summary statistics to quantitatively compare different estimators' precision and accuracy. +Since MCMC samples of these hyperparameters are Gaussian distributions, we can quantify the breadth of the distribution for each hyperparameter using the standard deviation regardless of whether the true values are known. -In simulated cases where the true parameter values are known, we calculate the Kullback-Leibler divergence (KLD), given by +In simulated cases where the true parameter values are known, we calculate the Kullback-Leibler divergence (KLD), given by \begin{equation} \label{eqn:kl} -KL_{\ndphi,\ndphi^{\ddagger}} = \integral{\pr{z \gvn \ndphi} \ln \left[ \frac{\pr{z \gvn \ndphi}}{\pr{z \gvn \ndphi^{\dagger}}} \right]}{z} , +KL\new{}{D}_{\ndphi,\ndphi^{\ddagger}} = \integral{\pr{z \gvn \ndphi} \ln \left[ \frac{\pr{z \gvn \ndphi}}{\pr{z \gvn \ndphi^{\dagger}}} \right]}{z} , \end{equation} -which measures a distance from parameter values $\ndphi$ to true parameter values $\ndphi^{\dagger}$. +which measures a distance from parameter values $\ndphi$ to true parameter values $\ndphi^{\dagger}$. The KLD is a measure of the information loss, in units of nats, due to using $\ndphi$ to approximate the true $\ndphi^{\dagger}$ when it is known. A detailed exploration of the KLD may be found in the Appendix to \citet{malz_approximating_2018}. -%We note that $KL_{\ndphi,\ndphi^{\ddagger}} \neq KL_{\ndphi^{\ddagger},\ndphi}$ and is only interpretable when there is a notion that $\ndphi^{\dagger}$ is closer to the truth than $\ndphi$. -%In simulated tests, $\ndphi^{\dagger}$ is the true value and $\ndphi'$ is that produced by one of the methods in question. - -%\aim{TODO: Also include the mean redshift/$\Delta_{z}$ for each estimator of \nz, cite the DES/KiDS/HSC papers that motivate this.} - -\subsection{Validation on mock data} -\label{sec:mock} \begin{figure*} \begin{center} - \includegraphics[width=0.7\textwidth]{figures/chippr/flowchart.pdf} + \includegraphics[width=0.7\textwidth]{flowchart.pdf} \caption{ - % \que{Change ``true joint model'' to true forward model, condition on $\phi^{\dagger}$, add footnote explaining distinction between dagger, star, prime, etc., check that derivation matches notation} A flow chart illustrating the forward model used to generate mock data in the validation of \Chippr, as described in \Sect{sec:forward}. Ovals indicate a quantity that must be chosen in order to generate the data, rectangles indicate an operation we perform, and rounded rectangles indicate a quantity created by the forward model. Arrows indicate the inputs and outputs of each operation performed to simulate mock \pzip\ catalogs. @@ -505,6 +472,22 @@ \subsection{Validation on mock data} \end{center} \end{figure*} +\subsection{Validation on mock data} +\label{sec:mock} + +\begin{figure}%[H] + \begin{center} + \includegraphics[width=0.4\textwidth]{data_lsst.png} + \caption{ + The joint probability space of true and estimated redshift for the three concerning \pz\ systematics at the level of the \lsst\ requirements: + intrinsic scatter, uniformly distributed catastrophic outliers, and bias. + The main panel shows samples (black points) in the space of mock data and redshift, akin to the standard scatterplots of true and estimated redshift, the $z_{\mathrm{spec}} = z_{\mathrm{phot}}$ diagonal (gray line), and posterior probabilities evaluated at the given estimated redshift (colored step functions). + The insets show marginal histograms (light gray) in each dimension, that can be compared with the true \nz\ used to make the figure (black) to see the effect of these systematics, as well as the implicit prior (dark gray). + } + \label{fig:mega_scatter} + \end{center} +\end{figure} + We compare the results of \Chippr\ to those of stacking and the histogram of \pzip\ maxima (modes) on mock data in the form of catalogs of emulated \pzip s generated via the forward model discussed in \Sect{sec:forward}. \Fig{fig:flowchart} illustrates the implementation of the forward model, defined by the much simpler \Fig{fig:pgm}, used for validating the method presented here. The irony of a simple model and complex validation procedure is not lost on the authors. @@ -521,88 +504,98 @@ \subsection{Validation on mock data} \label{eqn:gamma} n^{\dagger}(z) = \frac{1}{2 c_{z}} \left(\frac{z}{c_{z}}\right)^{2}\ \exp\left[-\frac{z}{c_{z}}\right] \end{equation} -with $c_{z} = 0.3$, because it has been used in forecasting studies for \des\ and \lsst. -%\aim{I learned this from talking to people and don't know of a published source that talks about the nitty gritty of the internal validation tests performed before there was data.} +with $c_{z} = 0.3$, because it has been used in forecasting studies for \lsst \new{}{\citep{chang_effective_2013, the_lsst_dark_energy_science_collaboration_lsst_2018}}. The mock data emulates the three sources of error of highest concern to the \pz\ community that are explored in detail later in this section: intrinsic scatter (\Sect{sec:scatter}), catastrophic outliers (\Sect{sec:outliers}), and canonical bias (\Sect{sec:bias}). \Fig{fig:mega_scatter} illustrates these three effects simultaneously at the tolerance of \lsst\ for demonstrative purposes, harking back to Figure~\ref{fig:pedagogical_scatter}. -\begin{figure}%[H] - \begin{center} - \includegraphics[width=0.45\textwidth]{figures/chippr/single_lsst-mega_scatter.png} - \caption{ - The joint probability space of true and estimated redshift for the three concerning \pz\ systematics at the level of the \lsst\ requirements: - intrinsic scatter, uniformly distributed catastrophic outliers, and bias. - The main panel shows samples (black points) in the space of mock data and redshift, akin to the standard scatterplots of true and estimated redshift, the $z_{\mathrm{spec}} = z_{\mathrm{phot}}$ diagonal (gray line), and posterior probabilities evaluated at the given estimated redshift (colored step functions). - The insets show marginal histograms (light gray) in each dimension, that can be compared with the true \nz\ used to make the figure (black) to see the effect of these systematics, as well as the implicit prior (dark gray). -% \aim{TODO: Include slice further up so we can see outliers. -% Label panel.} - } - \label{fig:mega_scatter} - \end{center} -\end{figure} - The hyperprior distribution chosen for these tests is a multivariate normal distribution with mean $\vec{\mu}$ equal to the implicit prior $\ndphi^{*}$ and covariance \begin{equation} \label{eqn:priorcov} -\Sigma_{k,k'} = q\ \exp[-\frac{e}{2}\ (\bar{z}_{k}-\bar{z}_{k'})^{2}]\ +\ t\delta(k,k') +\Sigma_{k,k'} = q\ \exp\left[-\frac{\new{e}{\varepsilon}}{2}\ (\bar{z}_{k}-\bar{z}_{k'})^{2}\right]\ +\ t\delta(k,k') \end{equation} -inspired by one used in Gaussian processes, where $k$ and $k'$ are indices ranging from $1$ to $K$ and $q=1.0$, $e=100.0$, and $t=q\cdot10^{-5}$ are constants chosen to permit draws from this prior distribution to produce shapes similar to that of a true $\tilde{\ndphi}$. +inspired by one used in Gaussian processes, where $k$ and $k'$ are indices ranging from $1$ to $K$ and $q = 1.0$, $\new{e}{\varepsilon} = 100.0$, and $t = q \cdot 10^{-5}$ are constants chosen to permit draws from this prior distribution to produce shapes similar to that of a true $\tilde{\ndphi}$. We adapt the full log-posterior of \Eq{eqn:final} to the chosen binning of redshift space. -%An example of such samples from the prior are shown in \Fig{fig:prior}. -%\que{Should I add back the figure of prior samples?} - -The sampler is initialized with $W=100$ walkers each with a value chosen from a Gaussian distribution of identity covariance around a sample from the hyperprior distribution. +The sampler is initialized with $W = 100$ walkers each with a value chosen from a Gaussian distribution of identity covariance around a sample from the hyperprior distribution. +\new{}{The computational expense is dominated by the sampler's burn-in phase and driven predominantly by the number of parameters $K$; +for the tests presented in this study, with $J = 10^{4}$ galaxies, $K = 35$ parameters, $W = 100$ walkers, a Gelman-Rubin convergence threshold of $R = 1.5$, a thinning factor of $s = 1$, and $\frac{I_{0}}{s} = 10^{3}$ accepted post-burn samples, the single-threaded computation time on a 2017-era laptop is on the order of a few hours.} +%{\color{white}{\Fig{fig:pzs-scatter} \Fig{fig:results-scatter} \Fig{fig:uniform-outliers } \Fig{fig:nonuniform-outliers-data} \Fig{fig:nonuniform-outliers-results} \Fig{fig:bias}}} \section{Results} \label{sec:alldata} Here, we compare the results of the \Chippr\ methodology with those of established \nz\ estimators under the three traditional measures of \pz\ uncertainty one at a time: \Sect{sec:scatter} concerns the redshift-dependent intrinsic scatter, \Sect{sec:outliers} concerns realistically complex catastrophic outlier populations, and \Sect{sec:bias} concerns the canonical bias in the mean redshift. -\subsection{Intrinsic scatter} -\label{sec:scatter} - -%Several factors contribute to photometric redshifts' intrinsic scatter. -%Distant galaxies are dimmer compared to galaxies of identical luminosity that are closer, driving up photometric errors in flux-limited surveys. -%The nature of the galaxy sample at higher redshifts also changes, meaning the generation of the photometric redshift posterior based on an a locally-calibrated SED template library or spectroscopically-confirmed training set is more likely to be inappropriate, leading to broader features. -%In general, the galaxies that could not have been observed spectroscopically will have different and noisier photo-$z$ likelihoods than those that could fall into a spectroscopic training set (or spectroscopically derived template library). -%This effect may be stronger for high-redshift galaxies. - -\Fig{fig:pzs-scatter} shows some examples of \pzpdf s generated with only the systematic of intrinsic scatter, at the level of the \lsst\ requirements on the left and twice that on the right. -One can see that the histogram of redshift estimates is broader than that of true redshifts, and that the effect is substantially more pronounced by just doubling the intrinsic scatter from the level of the \lsst\ requirements. +\begin{figure*} + \begin{center} + \includegraphics[width=0.4\textwidth]{data_scatter_lo.png} + \includegraphics[width=0.4\textwidth]{data_scatter_hi.png} + \caption{ + Examples of mock \pzpdf s generated with intrinsic scatter at the \lsst\ requirements (left) and twice the \lsst\ requirements (right), including samples from the probability space of true and observed redshift (black points), \pzpdf s (colored step functions), and the true redshifts of the example \pzpdf s (colored vertical lines). + A histogram (light gray) of points in each dimension is shown in the respective inset, with the true redshift distribution (black) and implicit prior (dark gray). + } + \label{fig:pzs-scatter} + \end{center} +\end{figure*} \begin{figure*} \begin{center} - \includegraphics[width=0.45\textwidth]{figures/chippr/single_varsigmas-mega_scatter.png} - \includegraphics[width=0.45\textwidth]{figures/chippr/thesis_hivarsig-mega_scatter.png} - \caption{ - Examples of mock \pzpdf s generated with intrinsic scatter at the \lsst\ requirements (left) and twice the \lsst\ requirements (right), including samples from the probability space of true and observed redshift (black points), \pzpdf s (colored step functions), and the true redshifts of the example \pzpdf s (colored vertical lines). - A histogram (light gray) of points in each dimension is shown in the respective inset, with the true redshift distribution (black) and implicit prior (dark gray). -% \aim{TODO: Label panels. -% Show the mean redshift for each estimator.} - } - \label{fig:pzs-scatter} + \includegraphics[width=0.4\textwidth]{results_scatter_lo.png} + \includegraphics[width=0.4\textwidth]{results_scatter_hi.png} + \caption{ + The results of \Chippr\ (samples in light blue and optimization in dark blue) and the alternative approaches (the stacked estimator in red and the histogram of modes in yellow) on \pzpdf s with intrinsic scatter of the \lsst\ requirements (left) and twice that (right), with the true redshift density (black curve) and implicit prior (gray curve). + \Chippr\ is robust to intrinsic scatter, but the alternatives suffer from overly broad \nz\ estimates that worsen with increasing intrinsic scatter. + } + \label{fig:results-scatter} \end{center} \end{figure*} +\subsection{Intrinsic scatter} +\label{sec:scatter} + +\Fig{fig:pzs-scatter} shows some examples of \pzpdf s generated with only the systematic of intrinsic scatter, at the level of the \lsst\ requirements on the left and twice that on the right. +One can see that the histogram of redshift estimates is broader than that of true redshifts, and that the effect is substantially more pronounced by just doubling the intrinsic scatter from the level of the \lsst\ requirements. + \Fig{fig:results-scatter} shows the \nz\ recovered by \Chippr\ and the alternative approaches. As expected, the estimates of \nz\ based on the modes of the \pzpdf s and stacking are broader than the marginalized maximum likelihood estimator from \chippr, with more broadening as the intrinsic scatter increases. \Chippr's \mmle\ is robust to intrinsic scatter and is unaffected by increased intrinsic scatter, though the \Chippr\ posterior distribution on the redshift distribution is itself broader for the higher intrinsic scatter case than for the \lsst\ requirements. The broadening of the alternative estimators corresponds to a loss of 3-4 times as many nats of information about \nz\ for the \lsst\ requirements relative to the \mmle\ of \Chippr. +\begin{figure} + \begin{center} + \includegraphics[width=0.4\textwidth]{data_outlier_uniform.png}\\ + \includegraphics[width=0.4\textwidth]{results_outlier_uniform.png} + \caption{ + Top: Examples of \pzpdf s with a uniformly distributed catastrophic outlier population at the level of the \lsst\ requirements, including samples from the probability space of true and observed redshift (black points), \pzpdf s (colored step functions), and the true redshifts of the example \pzpdf s (colored vertical lines), with marginal histograms (light gray) for each dimension with the true redshift distribution (black) and implicit prior (dark gray) in the insets. + Bottom: The results of \Chippr\ (samples in light blue, optimization in dark blue) and the alternative approaches (the stacked estimator in red, the histogram of modes in yellow) on \pzpdf s with uniformly distributed catastrophic outliers, with the true redshift density (black curve) and implicit prior (gray curve). + The presence of the catastrophic outlier population broadens the histogram of modes and stacked estimator of the redshift distribution, but the result of \Chippr\ is unbiased. + } + \label{fig:uniform-outliers} + \end{center} +\end{figure} + \begin{figure*} \begin{center} - \includegraphics[width=0.45\textwidth]{figures/chippr/single_varsigmas_log_estimators.png} - \includegraphics[width=0.45\textwidth]{figures/chippr/thesis_hivarsig_log_estimators.png} - \caption{ - The results of \Chippr\ (samples in light blue and optimization in dark blue) and the alternative approaches (the stacked estimator in red and the histogram of modes in yellow) on \pzpdf s with intrinsic scatter of the \lsst\ requirements (left) and twice that (right), with the true redshift density (black curve) and implicit prior (gray curve). - \Chippr\ is robust to intrinsic scatter, but the alternatives suffer from overly broad \nz\ estimates that worsen with increasing intrinsic scatter. -% \aim{TODO: Label panels. -% Show the mean redshift for each estimator.} - } - \label{fig:results-scatter} + \includegraphics[width=0.4\textwidth]{data_outlier_template.png} + \includegraphics[width=0.4\textwidth]{data_outlier_training.png} + \caption{ + Examples of \pzpdf s with a catastrophic outlier population like that seen in template-fitting \pzpdf\ codes (left) and machine learning \pzpdf\ codes (right), including samples from the probability space of true and observed redshift (black points), \pzpdf s (colored step functions), and the true redshifts of the example \pzpdf s (colored vertical lines), with marginal histograms (light gray) for each dimension with the true redshift distribution (black) and implicit prior (dark gray) in the insets. + } + \label{fig:nonuniform-outliers-data} + \end{center} +\end{figure*} + +\begin{figure*} + \begin{center} + \includegraphics[width=0.4\textwidth]{results_outlier_template.png} + \includegraphics[width=0.4\textwidth]{results_outlier_training.png} + \caption{ + The results of \Chippr\ (samples in light blue and optimization in dark blue) and the alternative approaches (the stacked estimator in red, the histogram of modes in yellow) on \pzpdf s with catastrophic outliers like those seen in template-fitting \pzpdf\ codes (left) and machine learning \pzpdf\ codes (right) to the \lsst\ requirements, with the true redshift density (black curve) and implicit prior (gray curve). + Though the histogram of modes is most sensitive to a catastrophic outlier population, the stacked estimator also overestimates \nz\ under (machine learning-like outliers) and beyond (template fitting-like outliers). + } + \label{fig:nonuniform-outliers-results} \end{center} \end{figure*} @@ -618,23 +611,6 @@ \subsection{Catastrophic outliers} The intrinsic scatter of the tests in this section does not increase with redshift as indicated in Table~\ref{tab:lsstsrd} in order to isolate the effect of outliers, and is instead held at a constant $\sigma_{z} = 0.02$. -\begin{figure} - \begin{center} - \includegraphics[width=0.45\textwidth]{figures/chippr/single_uout-mega_scatter.png}\\ - \includegraphics[width=0.45\textwidth]{figures/chippr/single_uout_log_estimators.png} - \caption{ - Top: Examples of \pzpdf s with a uniformly distributed catastrophic outlier population at the level of the \lsst\ requirements, including samples from the probability space of true and observed redshift (black points), \pzpdf s (colored step functions), and the true redshifts of the example \pzpdf s (colored vertical lines), with marginal histograms (light gray) for each dimension with the true redshift distribution (black) and implicit prior (dark gray) in the insets. -% \aim{TODO: Include slice further up. -% Label panels.} - Bottom: The results of \Chippr\ (samples in light blue, optimization in dark blue) and the alternative approaches (the stacked estimator in red, the histogram of modes in yellow) on \pzpdf s with uniformly distributed catastrophic outliers, with the true redshift density (black curve) and implicit prior (gray curve). -% \aim{TODO: Label panels. -% Show the mean redshift for each estimator.} - The presence of the catastrophic outlier population broadens the histogram of modes and stacked estimator of the redshift distribution, but the result of \Chippr\ is unbiased. - } - \label{fig:uniform-outliers} - \end{center} -\end{figure} - \Fig{fig:uniform-outliers} shows that at the level of the \lsst\ requirements, the alternative estimators are overly broad, whereas \Chippr's \mmle\ yields an unbiased estimate of \nz. Further, the result of stacking is even broader than that of the histogram of modes, corresponding to ten times the information loss of \Chippr's \mmle, making it worse than the most naive reduction of \pzpdf s to point estimates. @@ -646,38 +622,24 @@ \subsection{Catastrophic outliers} It is less straightforward to emulate catastrophic outliers like those anticipated of a machine learning code, those that are truly multimodal. The testing conditions here, illustrated in the right panel of \Fig{fig:nonuniform-outliers-data}, gives $10\%$ of galaxies at the redshift affected by outliers an observed redshift that is uniformly distributed relative to the true redshift, meaning that far fewer than $10\%$ of all galaxies in the sample are catastrophic outliers. -\begin{figure*} - \begin{center} - \includegraphics[width=0.45\textwidth]{figures/chippr/thesis_eout-mega_scatter.png} - \includegraphics[width=0.45\textwidth]{figures/chippr/thesis_rout-mega_scatter.png} - \caption{ - Examples of \pzpdf s with a catastrophic outlier population like that seen in template-fitting \pzpdf\ codes (left) and machine learning \pzpdf\ codes (right), including samples from the probability space of true and observed redshift (black points), \pzpdf s (colored step functions), and the true redshifts of the example \pzpdf s (colored vertical lines), with marginal histograms (light gray) for each dimension with the true redshift distribution (black) and implicit prior (dark gray) in the insets. -% \aim{TODO: Label panels. -% Include slices higher up.} - } - \label{fig:nonuniform-outliers-data} - \end{center} -\end{figure*} - The results of \Chippr\ and the alternative estimators of \nz\ are presented in \Fig{fig:nonuniform-outliers-results}. The most striking feature is that the histogram of modes is highly sensitive to both outlier populations, producing a severe overestimate in the case of an outlier population like those seen in template-fitting codes and a severe underestimate in the case of an outlier population like those seen in machine learning codes, corresponding to a twenty-fold loss of information compared to the \Chippr\ \mmle\ in both cases. The effect on the stacked estimator of \nz\ is more subtle though still concerning. In the case of outliers like those resulting from template-fitting, the stacked estimator is overly broad even without realistic intrinsic scatter, resulting in ten times the information loss compared to the \Chippr\ \mmle, and in the case of outliers like those resulting from machine learning, the stacked estimator features an overestimate at the redshift affected by the outlier population, resulting in about five times the information loss as the \Chippr\ \mmle. The \Chippr\ \mmle, however, appears unbiased and withstands these effects, and the breadth of the distribution of samples of \nz\ is invariant. -\begin{figure*} +\begin{figure} \begin{center} - \includegraphics[width=0.45\textwidth]{figures/chippr/thesis_eout_log_estimators.png} - \includegraphics[width=0.45\textwidth]{figures/chippr/thesis_rout_log_estimators.png} - \caption{ - The results of \Chippr\ (samples in light blue and optimization in dark blue) and the alternative approaches (the stacked estimator in red, the histogram of modes in yellow) on \pzpdf s with catastrophic outliers like those seen in template-fitting \pzpdf\ codes (left) and machine learning \pzpdf\ codes (right) to the \lsst\ requirements, with the true redshift density (black curve) and implicit prior (gray curve). - Though the histogram of modes is most sensitive to a catastrophic outlier population, the stacked estimator also overestimates \nz\ under (machine learning-like outliers) and beyond (template fitting-like outliers). -% \aim{TODO: Label panels. -% Show the mean redshift for each estimator.} - } - \label{fig:nonuniform-outliers-results} + \includegraphics[width=0.4\textwidth]{data_bias.png}\\ + \includegraphics[width=0.4\textwidth]{results_bias.png} + \caption{ + Top: Examples of \pzpdf s with ten times the bias of the \lsst\ requirements, including samples from the probability space of true and observed redshift (black points), \pzpdf s (colored step functions), and the true redshifts of the example \pzpdf s (colored vertical lines), with marginal histograms (light gray) for each dimension with the true redshift distribution (black) and implicit prior (dark gray) in the insets. + Bottom: The results of \Chippr\ (samples in light blue, optimization in dark blue) and the alternative approaches (the stacked estimator in red, the histogram of modes in yellow) on \pzpdf s with ten times the bias of the \lsst\ requirements, with the true redshift density (black curve) and implicit prior (gray curve). + The impact of bias at even ten times the level of the \lsst\ requirements is almost imperceptible on all estimators, though the \Chippr\ \mmle\ minimizes the information loss regardless. + } + \label{fig:bias} \end{center} -\end{figure*} +\end{figure} \subsection{Canonical bias} \label{sec:bias} @@ -690,252 +652,196 @@ \subsection{Canonical bias} Consider that if the canonical bias were included in the framework of Figure~\ref{fig:pedagogical_scatter}, it could be trivially modeled out as a simple linear transformation of $z_{\mathrm{phot}} \to z_{\mathrm{phot}} - \Delta_{z} (1 + z_{\mathrm{phot}})$ of the $(z_{\mathrm{spec}}, z_{\mathrm{phot}})$ space. Regardless, for completeness, a test at ten times the canonical bias of the \lsst\ requirements, with no redshift-dependent intrinsic scatter nor catastrophic outliers, is provided in \Fig{fig:bias}. -\begin{figure} - \begin{center} - \includegraphics[width=0.45\textwidth]{figures/chippr/thesis_neghivarbias-mega_scatter.png}\\ - \includegraphics[width=0.45\textwidth]{figures/chippr/thesis_neghivarbias_log_estimators.png} - \caption{ - Top: Examples of \pzpdf s with ten times the bias of the \lsst\ requirements, including samples from the probability space of true and observed redshift (black points), \pzpdf s (colored step functions), and the true redshifts of the example \pzpdf s (colored vertical lines), with marginal histograms (light gray) for each dimension with the true redshift distribution (black) and implicit prior (dark gray) in the insets. -% \aim{TODO: Include slice further up. -% Label panels.} - Bottom: The results of \Chippr\ (samples in light blue, optimization in dark blue) and the alternative approaches (the stacked estimator in red, the histogram of modes in yellow) on \pzpdf s with ten times the bias of the \lsst\ requirements, with the true redshift density (black curve) and implicit prior (gray curve). -% \aim{TODO: Label panel. -% Show the mean redshift for each estimator.} - The impact of bias at even ten times the level of the \lsst\ requirements is almost imperceptible on all estimators, though the \Chippr\ \mmle\ minimizes the information loss regardless. - } - \label{fig:bias} - \end{center} -\end{figure} - As expected based on self-consistency of the forward-modeled \pzpdf s, \Chippr\ is immune to linear bias of the form of $\Delta_{z}$. Furthermore, the alternative estimators are only weakly affected, with information loss two and four times greater than that of the \Chippr\ \mmle\ for the histogram of modes and stacked estimator respectively. (This general robustness may suggest that the canonical bias may not be the most relevant measure of performance of estimators of \nz.) +%{\color{white}{\Fig{fig:pzs-priors} \Fig{fig:results-priors} \Fig{fig:mischaracterized}}} \section{Discussion} \label{sec:results} The experiments of \Sect{sec:alldata} quantify the influence on each estimator of \nz\ due to each of the canonical types of \pz\ error one at a time in isolation. -Now, we stress-test \Chippr\ by exploring the impact of the implicit prior, which has thus far not received much attention in the literature. -%two realistically complex cases, one in which the \nz\ estimates are made tomographically as in a modern cosmological analysis (\Sect{sec:lsstdemo}) and one +Now, we \new{}{forecast the influence of these \nz\ estimators on cosmological parameter constraints in a tomographic analysis and} stress-test \Chippr\ by exploring the impact of the implicit prior, which has thus far not received attention in the literature. +%(\Sect{sec:lsstdemo}) \Sect{sec:interim} demonstrates the sensitivity of \nz\ estimation methods to realistically complex implicit priors, and \Sect{sec:violations} demonstrates the consequences of mischaracterization of the implicit prior used to generate the \pzip\ catalog. +\new{}{\Sect{sec:lsstdemo} presents the impact of \nz\ estimation methods on forecasted cosmological parameter constraints.} These results provide compelling motivation for the \pz\ community to prioritize the study of implicit priors of existing and developing \pzpdf\ techniques. -%\que{Add back the results of the LSST requirements here?} - -%\subsection{LSST Requirements} -%\label{sec:lsstdemo} -% -%It is of interest to explore the impact of incorrectly estimated \nz\ on the cosmological inference to answer the question of how wrong we will be in our understanding of the universe if we incorrectly constrain \nz. -%To test the impact of these uncertainties, we simulate mock data with all three effects with which \lsst\ is concerned at the levels of Table~\ref{tab:lsstsrd} and propagate the results of \Chippr\ and the other estimators to a Fisher matrix forecast using \cosmolike\ \citep{krause_cosmolike_2017}, a publicly available cosmological forecasting code. -% -%\begin{figure} -% \begin{center} -% \includegraphics[width=0.45\textwidth]{figures/chippr/cosmolike_inputs.png} -% \caption{ -% The \lsst-like tomographic binning and true redshift distribution, where the truth (solid) is a PDF evaluated on a fine grid of $350$ redshifts $0.0101 < z < 3.5001$, and the binned (dashed) and drawn (dotted) \nz\ are piecewise constant functions evaluated in $35$ evenly spaced bins, for four different tomographic bins (colors). -% } -% \label{fig:tomobins} -% \end{center} -%\end{figure} -% -%\dwh{We consider as ground truth a set of known \nz\ corresponding to each of four hypothetical samples of galaxies and the corresponding cosmological parameter covariance matrix. -%The \nz\ of each galaxy subsample emulates that anticipated of galaxies binned by a redshift point estimate, as is common in tomographic redshift analyses, though our experimental procedure is agnostic to how the samples are identified. -%The cosmological parameter covariance matrices are those used for \desc\ forecasting with the ground truth \nz\ in the same four bins.} -%The true \nz\ in each pre-defined bin is already provided in the form of an evaluation of a function on a fine grid of $350$ redshifts $0.0101 < z < 3.5001$. -% -%First, we bin them down to a piecewise constant parameterization with a manageable $35$ hyperparameters for \chippr's sampling capabilities. -%Next, we draw $10^{4}$ true redshifts from the binned true \nz\ for each tomographic bin. -%The original, binned, and drawn \nz\ are shown in \Fig{fig:tomobins}. -%We emulate \pzpdf s for the $10^{4}$ true redshifts drawn from the true \nz\ in each bin using the procedure of \Fig{fig:flowchart} with all three effects of Table~\ref{tab:lsstsrd} -%at their given levels. -%Illustrations of this process are provided in \Fig{fig:per-bin-scatter}. -% -%\begin{figure*} -% \begin{center} -% \includegraphics[width=0.24\textwidth]{figures/chippr/0single_lsst_mega_scatter.png} -% \includegraphics[width=0.24\textwidth]{figures/chippr/1single_lsst_mega_scatter.png} -% \includegraphics[width=0.24\textwidth]{figures/chippr/2single_lsst_mega_scatter.png} -% \includegraphics[width=0.24\textwidth]{figures/chippr/3single_lsst_mega_scatter.png} -% \caption{As in \Fig{fig:mega_scatter}, with a different tomographic bin in each panel and the three effects of intrinsic scatter, uniformly distributed catastrophic outliers, and bias at the levels of the \lsst\ SRD, given in Table~\ref{tab:lsstsrd}. -% \aim{TODO: Make this one big plot instead of four little ones to eliminate repeated insets and legend. -% Enlarge axis labels. -% Label panels. -% Also, add watermark of ``mock data'' in UL corner, same for ``results of inference'' on other kind of plot. -% Show the mean redshift for each estimator of N(z), cite the DES papers that motivate this.} -% } -% \label{fig:per-bin-scatter} -% \end{center} -%\end{figure*} -% -%\que{Is the distinction between binned samples determined by some observational property and probabilistic redshift distributions sufficiently clear?} -% -%We then make a point estimate of \nz\ using \chippr's \mmle\ optimization option as well as the alternative methods on the \pzpdf\ catalog for each tomographic bin, shown in \Fig{fig:per-bin-ests}, because \cosmolike\ produces cosmology constraints from a single \nz\ result, rather than samples from the full posterior probability density of possible \nz. -%Note that \Fig{fig:per-bin-ests} is shown in linear rather than log probability units, unlike all other plots in this paper, to better show the behavior at low probability. -%The excessive breadth of the alternative estimators can be seen quite plainly. -% -%\begin{figure*} -% \begin{center} -% \includegraphics[width=0.24\textwidth]{figures/chippr/0single_lsst_lin_estimators.png} -% \includegraphics[width=0.24\textwidth]{figures/chippr/1single_lsst_lin_estimators.png} -% \includegraphics[width=0.24\textwidth]{figures/chippr/2single_lsst_lin_estimators.png} -% \includegraphics[width=0.24\textwidth]{figures/chippr/3single_lsst_lin_estimators.png} -% \caption{ -% The \chippr-derived and other estimators of \nz\ in each tomographic bin, with the true \nz\ (black), the implicit prior (gray), stacked estimator (red), histogram of modes (yellow), and \Chippr\ \mmle\ (blue). -% The result of stacking is far too broad for \lsst-like \pzpdf s, even moreso than the simplistic histogram of modes. -% \aim{TODO: Make this one big plot instead of four little ones to eliminate repeated axis labels and legend. -% Label figure/panels. -% Also, add watermark of ``results of inference'' in UL corner, same for ``mock data'' on other kind of plot. -% Show the mean redshift for each estimator of N(z), cite the DES papers that motivate this.} -% } -% \label{fig:per-bin-ests} -% \end{center} -%\end{figure*} -% -%We then use the different estimators of \nz\ in a cosmological forecasting procedure with \cosmolike, constraining $\Omega_{m}$, $\Omega_{b}$, $w_{a}$, $w_{0}$, $n_{s}$, $S_{8}$, and $H_{0}$. -%Though there are also slight differences in the angle of the error ellipses, the most striking effect is the broadening of the contours under the alternative estimators relative to \Chippr, which are almost indistinguishable from those derived by using the true redshift distribution in each bin. -%The stacked estimator is significantly worse than the \Chippr\ \mmle\ for all parameters except $\Omega_{b}$ and $H_{0}$. -%Stacking, however, outperforms the histogram of modes for all parameters except $\Omega_{m}$ and $S_{8}$, for which their constraints are quite similar. -%Though the true values of the parameters themselves were not accessible with the Fisher matrix-based framework, we calculate the seven-dimensional KLD for the three \nz\ estimators relative to the constraints derived from the true \nz, showing that \Chippr\ preserves information $200-800$ times better than the alternatives, with the histogram of modes doing about four times better than stacking. -%\que{So a Fisher matrix analysis inherently can't provide the bias, because that requires data. -%Because we have no data, only posteriors conditioned on hypothetical data, this isn't possible. -%Do you think it's sufficient to provide the bias on the moments of \nz, since that's what everyone uses anyway?} -% -%\begin{figure*} -% \begin{center} -% \includegraphics[width=0.9\textwidth]{figures/chippr/final_plot.png} -% \caption{ -% \que{Are the contours any easier to see now?} -% The result of propagating the estimators of \nz\ by stacking (red), the histogram of modes (yellow), \Chippr\ (blue), and the true \nz\ (black) of \Fig{fig:per-bin-ests} to a subset of cosmological parameters. -% For all parameters considered, \Chippr\ yields contours no broader than those corresponding to the true \nz, whereas for most parameters, stacking and the histogram of modes yield broader contours. -% \aim{TODO: Fix formatting of axis labels. -% Standardize ticklabels.} -% } -% \label{fig:cornerplot} -% \end{center} -%\end{figure*} -% -%\que{Does this discussion adequately quantify how good \Nz\ has to be and how wrong we'll be if we estimate it wrong? -% Suggestions for how to better establish context would be appreciated.} - -\subsection{Realistically complex implicit prior} -\label{sec:interim} - -\chippr\ can handle any implicit prior with support over the redshift range where \nz\ is defined, but some archetypes of implicit prior are more likely to be encountered in the wilds of \pzip\ codes. -Ideally, an uninformative implicit prior would be used, although it may be complicated to compute from the covariances of the raw data. -Template-fitting codes have an explicit prior input formed by redshifting a small number of templates, leading to a highly nonuniform but physically-motivated interim prior. -%Another potential method for selecting an interim prior with support over the entire redshift range expected of the photometric survey is to sum two or more $N(z)$ distributions obtained from reliable photometric surveys in the past. -%This is just as problematic as using a biased spectroscopically derived $N(z)$ as the interim prior because the sum of redshift distributions for two or more surveys does not reflect our beliefs about the true distribution for a single survey even though it provides support over the same redshift range. -%To simulate this case, we choose an interim prior with more weight at high and low redshifts than for mid-range redshifts. -Machine learning approaches tend to be trained on previously observed data sets that are biased towards low redshift, which biases the implicit prior towards low redshift. -Some efforts have been made to modify an observationally informed implicit prior so that it is more representative of the photometric data for which redshifts are desired \citep{sheldon_photometric_2012}, but, unless it is equal to the true \nz, it will propagate to the results of traditional \nz\ estimation methods. -%Because low-redshift galaxies are more likely to be bright enough to be observed by such a survey, $N(z)$ determined from that sample may be heavily biased to low redshift galaxies. -%By contrast, the galaxies that were unobserved in such a survey are more likely be dimmer, making them more likely to be at higher redshifts. -%Since the interim prior is not compatible with our beliefs about the true redshift distribution, the resulting interim redshift posteriors will be inappropriate. - -\Fig{fig:pzs-priors} shows examples of \pzip s with a low-redshift favoring implicit prior emulating that of a machine learning approach to \pz\ estimation (left panel) and a more complex interim prior emulating that of a template-fitting \pz\ method (right panel). -One can see that the \pzip s take different shapes from one another even though the marginal histograms of the points are identical. -The machine learning-like implicit prior has been modified to have nonzero value at high-redshift because the implicit prior must be strictly positive definite for the \Chippr\ model to be valid. - \begin{figure*} \begin{center} - \includegraphics[width=0.45\textwidth]{figures/chippr/single_lsst_trpr-mega_scatter.png} - \includegraphics[width=0.45\textwidth]{figures/chippr/single_lsst_tmpr-mega_scatter.png} + \includegraphics[width=0.4\textwidth]{data_prior_training.png} + \includegraphics[width=0.4\textwidth]{data_prior_template.png} \caption{ Examples of mock \pzip s generated with a machine learning-like implicit prior (left) and a template-fitting-like implicit prior (right), including samples from the probability space of true and observed redshift (black points), \pzip s (colored step functions), the true redshifts of the example \pzip s (colored vertical lines). A histogram (light gray) of points in each dimension is shown in the respective inset, with the true redshift distribution (black) and implicit prior (dark gray). -% \aim{TODO: Label panels. -% Include slices higher up.} } \label{fig:pzs-priors} \end{center} \end{figure*} -\Fig{fig:results-priors} shows the performance of \Chippr\ and the traditional methods on \pzip s generated with nontrivial implicit priors. -In both cases, the \Chippr\ \mmle\ effectively recovers the true redshift distribution, and the distribution of \nz\ parameter values reflects higher uncertainty where the implicit prior undergoes large changes in derivative. -The alternatives, on the other hand, are biased by the implicit prior except where it is flat, in the case of high redshifts for the machine learning-like implicit prior, resulting in over $1,000$ times the information loss on \nz\ for the machine learning-like implicit prior and some $5-20$ times the information loss for the template fitting-like implicit prior, relative to the \Chippr\ \mmle. - \begin{figure*} \begin{center} - \includegraphics[width=0.45\textwidth]{figures/chippr/single_lsst_trpr_log_estimators.png} - \includegraphics[width=0.45\textwidth]{figures/chippr/single_lsst_tmpr_log_estimators.png} + \includegraphics[width=0.4\textwidth]{results_prior_training.png} + \includegraphics[width=0.4\textwidth]{results_prior_template.png} \caption {The results of \Chippr\ (samples in light blue and optimization in dark blue) and the alternative approaches (the stacked estimator in red and the histogram of modes in yellow) on \pzip s with an implicit prior like that of machine learning \pzip\ approaches (left) and an implicit prior like that of template-fitting \pzip\ codes (right), with the true redshift density (black curve) and implicit prior (gray curve). \Chippr\ is robust to a nontrivial implicit prior, but the alternatives are biased toward the implicit prior. -% \aim{TODO: Label panels. -% Show the mean redshift for each estimator.} } \label{fig:results-priors} \end{center} \end{figure*} -The main implication of the response of \nz\ estimates to a nontrivial implicit prior is that the implicit prior must be accounted for when using \pzip\ catalogs. +\subsection{Realistically complex implicit prior} +\label{sec:interim} -\subsection{Violations of the model} -\label{sec:violations} +\chippr\ can handle any implicit prior with support over the redshift range where \nz\ is defined, but some archetypes of implicit prior are more likely to be encountered in the wilds of \pzip\ codes. +Ideally, an uninformative implicit prior would be used, although it may be complicated to compute from the covariances of the raw data. +Template fitting codes have an explicit prior input formed by redshifting a number of templates, leading to a highly nonuniform but physically-motivated interim prior. +Machine learning approaches tend to be trained on one of more previously observed data sets that include only galaxies for which spectroscopy is accessible, typically biasing the implicit prior towards atypically bright and/or low redshift populations. +Some efforts have been made to modify an observationally informed implicit prior so that it is more representative of the photometric data for which redshifts are desired \citep{sheldon_photometric_2012}, but, unless it is equal to the true \nz, it will propagate to the results of traditional \nz\ estimation methods. -In this test, the \pzip s are made to the \lsst\ requirements but the implicit prior used for the inference is not the same as the implicit prior used for generating the data. -\Pzpdf\ codes do not generally provide their implicit prior, with the exception of some template-fitting techniques for which it is a known input. -If we naively used the \pzip\ catalog produced by a generic machine learning or template-fitting code and assumed a flat implicit prior, we would observe the contents of \Fig{fig:mischaracterized}. +\Fig{fig:pzs-priors} shows examples of \pzip s with a low-redshift favoring implicit prior emulating that of a machine learning approach to \pz\ estimation (left panel) and a more complex interim prior emulating that of a template-fitting \pz\ method (right panel). +One can see that the \pzip s take different shapes from one another even though the marginal histograms of the points are identical. +The machine learning-like implicit prior has been modified to have nonzero value at high-redshift because the implicit prior must be strictly positive definite for the \Chippr\ model to be valid. + +\Fig{fig:results-priors} shows the performance of \Chippr\ and the traditional methods on \pzip s generated with nontrivial implicit priors. +In both cases, the \Chippr\ \mmle\ effectively recovers the true redshift distribution, and the distribution of \nz\ parameter values reflects higher uncertainty where the implicit prior undergoes large changes in derivative. +The alternatives, on the other hand, are biased by the implicit prior except where it is flat, in the case of high redshifts for the machine learning-like implicit prior, resulting in over $1,000$ times the information loss on \nz\ for the machine learning-like implicit prior and some $5-20$ times the information loss for the template fitting-like implicit prior, relative to the \Chippr\ \mmle. + +The main implication of the response of \nz\ estimates to a nontrivial implicit prior is that the implicit prior must be accounted for when using \pzip\ catalogs. \begin{figure*} \begin{center} - \includegraphics[width=0.45\textwidth]{figures/chippr/single_lsst_trpr_wrong_log_estimators.png} - \includegraphics[width=0.45\textwidth]{figures/chippr/single_lsst_tmpr_wrong_log_estimators.png} + \includegraphics[width=0.4\textwidth]{results_wrongprior_training.png} + \includegraphics[width=0.4\textwidth]{results_wrongprior_template.png} \caption {The results of \Chippr\ (samples in light blue, optimization in dark blue) and the alternative approaches (the stacked estimator in red, the histogram of modes in yellow) when run with an incorrectly specified implicit prior (gray curve). The data upon which each panel's results are based are provided in Figure~\ref{fig:pzs-priors}, where the left corresponds to the sort of implicit prior anticipated of machine learning approaches and the right corresponds to an implicit prior like that of a template-fitting code. Here, \Chippr\ has been provided with a uniform implicit prior rather than those used to produce the mock \pzip s, and its performance is notably worse than when it is provided an accurate implicit prior, as in Figure~\ref{fig:results-priors}. -% \aim{Label panels. -% Show the mean redshift for each estimator.} When the incorrect implicit prior is provided to \chippr, even Bayesian inference cannot recover the true \nz. } \label{fig:mischaracterized} \end{center} \end{figure*} +\subsection{Violations of the model} +\label{sec:violations} + +In this test, the \pzip s are made to the \lsst\ requirements but the implicit prior used for the inference is not the same as the implicit prior used for generating the data. +\Pzpdf\ codes do not generally provide their implicit prior, with the exception of some template-fitting techniques for which it is a known input. +If we naively used the \pzip\ catalog produced by a generic machine learning or template-fitting code and assumed a flat implicit prior, we would observe the contents of \Fig{fig:mischaracterized}. + The results of using a mischaracterized implicit prior are disastrous, causing every estimator, including \Chippr, to be strongly biased. The stacked estimator and histogram of modes don't make use of the implicit prior so do no worse than when the implicit prior is accurately provided, but \Chippr\ is sensitive to prior misspecification, which violates the model upon which it is based. It is thus crucial that \pzip\ methods always characterize and provide the implicit prior. -% removed ancient investigation on real data, would be hard to redo with the new \chippr\ code on this timescale given that I haven't even looked at the data format in 4 years +\new{}{ + \begin{figure} + \begin{center} + \includegraphics[width=0.45\textwidth]{combo_mega_scatter.png} + \caption{\new{}{ + As in \Fig{fig:mega_scatter}, for each of four \desc-like tomographic bins (colors). + Main panel: $(z_{spec}, z_{phot})$ pairs frawn from a model with \lsst-like systematic errors. + Top panel: The \lsst-like tomographic \nz\ model (thin, opaque line), after a reduction in the number of parameters of \nz\ (dashed line), and histogram of $z_{spec}$ values drawn from that model (thick, semitransparent line). + Side panel: Histograms of $z_{spec}$ (thin, opaque line) and $z_{phot}$ (semitransparent bars). + }} + \label{fig:per-bin-scatter} + \end{center} + \end{figure} + + \begin{figure} + \begin{center} + \includegraphics[width=0.4\textwidth]{binned_nz_recovered.png} + \caption{\new{}{ + The \chippr-derived and other estimators of \nz\ in each tomographic bin (panels), with the true \nz\ (solid black), stacked estimator (thin dashed red), histogram of modes (dotted yellow), and \Chippr\ \mmle\ (thick dashed blue). + The result of stacking is far too broad for \lsst-like \pzpdf s, even moreso than the simplistic histogram of modes, while the \Chippr\ \mmle\ more closely follows the amplitude of the true \nz. + }} + \label{fig:per-bin-ests} + \end{center} + \end{figure} + + \begin{figure*} + \begin{center} + \includegraphics[width=0.9\textwidth]{final_plot.png} + \caption{\new{}{ + The result of propagating the estimators of \nz\ by stacking (red horizontal hatching), the histogram of modes (yellow vertical hatching), \Chippr\ (blue shaded), and the true \nz\ (black outline) of \Fig{fig:per-bin-ests} to a subset of cosmological parameters. + For all parameters considered, \Chippr\ yields contours no broader than those corresponding to the true \nz, whereas stacking and the histogram of modes yield broader contours. + }} + \label{fig:cornerplot} + \end{center} + \end{figure*} + + \subsection{Propagation to cosmology} + \label{sec:lsstdemo} + + It is of interest to explore the impact of incorrectly estimated \nz\ on derived constraints on the cosmological parameters to answer the question of how wrong we will be in our understanding of the universe if we incorrectly constrain \nz. + To test the impact of these uncertainties, we emulate mock data with all three effects with which \lsst\ is concerned at the levels of Table~\ref{tab:lsstsrd} and propagate the results of \Chippr\ and the other estimators of \nz\ through a Fisher matrix forecast using \cosmolike\ \citep{krause_cosmolike_2017}, a publicly available cosmological forecasting code. + + %\begin{figure} + % \begin{center} + % \includegraphics[width=0.45\textwidth]{combo_mega_scatter.png} + % \caption{\new{}{ + % The \lsst-like tomographic binning and true redshift distribution, where the truth (solid) is a PDF evaluated on a fine grid of $350$ redshifts $0.0101 < z < 3.5001$, and the binned (dashed) and drawn (dotted) \nz\ are piecewise constant functions evaluated in $35$ evenly spaced bins, for four different tomographic bins (colors). + % }} + % \label{fig:tomobins} + % \end{center} + %\end{figure} + + We consider as ground truth a set of known \nz\ corresponding to each of four hypothetical samples of galaxies and the corresponding cosmological parameter covariance matrix. + The \nz\ of each galaxy subsample emulates that anticipated of galaxies binned by a redshift point estimate, as is common in tomographic redshift analyses, though our experimental procedure is agnostic to how the samples are identified. + The cosmological parameter covariance matrices are those used for \desc\ forecasting with the ground truth \nz\ in the same four bins. +% \aim{TODO: put all relevant files on GitHub.} + The true \nz\ in each pre-defined bin is provided in the form of an evaluation of the function on a fine grid of $350$ redshifts $0.0101 < z < 3.5001$. + + First, we bin down the true redshift distribution to a piecewise constant parameterization with a manageable $35$ hyperparameters for \chippr's sampling capabilities. + Next, we draw $10^{4}$ true redshifts from the binned true \nz\ for each tomographic bin. + The original, binned, and drawn \nz\ are shown in the top panel of \Fig{fig:per-bin-scatter}. + We emulate \pzpdf s for the $10^{4}$ true redshifts drawn from the true \nz\ in each bin using the procedure of \Fig{fig:flowchart} with all three effects of Table~\ref{tab:lsstsrd} at their given levels. + The drawn $(z_{spec}, z_{phot})$ pairs are shown in the main panel of \Fig{fig:per-bin-scatter}, and the side panel shows their marginal histograms. + %Illustrations of this process are provided in \Fig{fig:per-bin-scatter}. + + We then make a point estimate of \nz\ using \chippr's \mmle\ optimization option as well as the alternative methods on the \pzpdf\ catalog for each tomographic bin, shown in \Fig{fig:per-bin-ests}, because \cosmolike\ produces cosmology constraints from a single \nz\ result, rather than samples from the full posterior probability density of possible \nz. + %Note that \Fig{fig:per-bin-ests} is shown in linear rather than log probability units, unlike all other plots in this paper, to better show the behavior at low probability. + The excessive breadth of the alternative estimators, especially in the regions of low probability at high redshift, can be seen quite plainly. + + %\subsubsection{Impact on cosmology constraints} + %\label{sec:cosmo} + + We then use the different estimators of \nz\ in a cosmological forecasting procedure with \cosmolike, constraining $\Omega_{m}$, $\Omega_{b}$, $w_{a}$, $w_{0}$, $n_{s}$, $S_{8}$, and $H_{0}$, shown in \Fig{fig:cornerplot}. + Though there are also slight differences in the angle of the error ellipses, the most striking effect is the broadening of the contours under the alternative estimators relative to \Chippr, which are almost indistinguishable from those derived by using the true redshift distribution in each bin. + The stacked estimator is significantly worse than the \Chippr\ \mmle\ for all parameters except $\Omega_{b}$ and $H_{0}$. + Stacking, however, outperforms the histogram of modes for all parameters except $\Omega_{m}$ and $S_{8}$, for which their constraints are quite similar. + Though the Fisher matrix approach cannot be used to evaluate the bias of point estimates of the cosmological parameters relative to their true values, we instead calculate the seven-dimensional KLD for the three \nz\ estimators relative to the constraints derived from the true \nz, finding that \Chippr\ preserves information $200-800$ times better than the alternatives, and even the histogram of modes outperforms stacking by a factor of four. + +} \section{Conclusion} \label{sec:con} +%\aim{TODO: Break up conclusion into subsections of \Sect{sec:alldata}/\Sect{sec:results}?} -%\que{TODO: Break up conclusion into subsections of \Sect{sec:alldata}/\Sect{sec:results}?} - -This study derives and demonstrates a mathematically consistent inference of a one-point statistic, the redshift density function \nz, based on an arbitrary catalog of \pzpdf s. -The fully Bayesian \Chippr\ model, based in the fundamental laws of probability, begins with a probabilistic graphical model corresponding to equations for the full posterior distribution over the parameters for \nz. -The \Chippr\ model is implemented in the publicly available \chippr\ code. -The method is implemented in the publicly available \chippr\ code and validated on mock data. +This study derives and demonstrates a mathematically consistent inference of a one-point statistic, the redshift density function \nz, based on an arbitrary catalog of \pzpdf s. +The fully Bayesian \Chippr\ model\new{, based in the fundamental laws of probability,}{} begins with a probabilistic graphical model corresponding to equations for the full posterior distribution over the parameters for \nz. +The \Chippr\ model is implemented in the publicly available \chippr\ \new{}{prototype} code and validated on mock data. % at the level of \nz\ as well as in terms of constraining power on the cosmological parameters. -%In the tests on simulated data performed here, the full posterior distribution over the hyperparameters defining $N(z)$ derived by this method is consistent with the true redshift distribution function, making the mean of sampled values an excellent point estimator of $N(z)$. -%The information contained in the full posterior distribution's shape convey the traditional error bar information without having to explicitly propagate any error estimates. -%The results of those tests is summarized below and in \Tab{tab:kld}, where lower values indicate a closer match between the true $N(z)$ and the estimator. -%Tests were also performed on subsets of BOSS DR10 data with results consistent with those of simulations. - Using a flexible, self-consistent forward model of the relationship between true and estimated redshifts, capable of encapsulating the complexity of observed redshift-photometry relations (e.g. \Fig{fig:pedagogical_scatter}), we emulate the canonical \pz\ error statistics, intrinsic scatter (\Sect{sec:scatter}), catastrophic outliers (\Sect{sec:outliers}), and canonical bias (\Sect{sec:bias}) one at a time. -Though these test cases may appear overly simplistic, they enable rigorous quantification of the relative performance of each \nz\ estimation techniques under the controlled conditions of each type of error in isolation, at levels equal to and beyond those of \lsst. -%\aim{TODO: point out that Fig 1 is uglier than 4, 5, 7, 8, 10, 11; these tests are simplistic, reality is more complex, but because \Chippr is provably correct, it can handle complexity of the real thing if p(z) is correct and captures complexity} +Though these test cases may appear overly simplistic, they enable rigorous quantification of the \new{relative}{} performance of each \new{}{of two established} \nz\ estimation techniques \new{}{relative to \Chippr} under the controlled conditions of each type of error in isolation, at levels equal to and beyond those of \lsst. Based on our tests, the following statements about the \Chippr\ methodology may be made with confidence: \begin{itemize} -\item \Chippr\ outperforms traditional estimators of \nz\ under realistically complex conditions, even at pessimistic levels relative to future survey requirements on the traditional \pz\ error statistics, as demonstrated both by eye and according to KLD values corresponding to $10\%$ the information loss of alternative methods. -%\aim{TODO: Refer to quantitative results on KLD of \Nz.} +\item \Chippr\ outperforms traditional estimators of \nz\ under realistically complex \new{conditions}{systematics}, even at pessimistic levels relative to future survey requirements on the \new{traditional}{standard} \pz\ error statistics, \new{as demonstrated both by eye and according to KLD values corresponding to}{achieving $\sim$} $10\%$ the information loss of alternative methods. \item Both the \Chippr\ \mmle\ and the mean of \chippr\ samples are good point estimators of \nz, whereas the histogram of modes is very sensitive to outliers and the stacked estimator is always excessively broad. -\item The error bars on the posterior distribution over \nz\ hyperparameters are interpretable and arise naturally under \Chippr, unlike those that may be assumed for the conventional point estimators. +\item The error bars \new{on}{of} the posterior distribution over \nz\ hyperparameters are interpretable and arise naturally under \Chippr, unlike those that \new{may be}{are sometimes} assumed for the conventional point estimators. +\item \new{}{When the \nz\ estimates for LSST-like data are propagated through a cosmological forecast, \Chippr's information loss in the $(\Omega_{m}, \Omega_{b}, w_{a}, w_{0}, n_{s}, S_{8}, H_{0})$ parameter space is 2-3 times lower than that corresponding to traditional \nz\ estimators.} \end{itemize} -Not only is \Chippr\ the only mathematically correct approach to the problem, it also recovers the true values of the hyperparameters defining \nz\ better than popular alternatives, as measured by the loss of information in \nz. -% \ and the size of error ellipses in the space of cosmological parameters. -However, the mathematically valid approach to inference with probabilistic data products incurs nontrivial computational expense, motivating future work to optimize the implementation. +Not only is \Chippr\ the only mathematically correct approach to the problem\new{}{ of obtaining \nz\ from an external catalog of \pzpdf s}, it also recovers the true values of \new{}{both} the hyperparameters defining \nz \new{}{and the cosmological parameters} better than popular alternatives, as measured by the \new{loss of information in \nz}{KLD}. +However, the mathematically valid approach to inference with probabilistic data products incurs nontrivial computational expense, motivating \new{future}{ongoing} work to optimize the implementation. -Additionally, this work highlights a crucial and almost entirely overlooked complication to the usage of \pzpdf s, namely the implicit prior, motivating the following recommendations: +Additionally, this work highlights a crucial and almost entirely overlooked complication to the usage of \pzpdf s, namely the implicit prior, motivating the following recommendations: \begin{itemize} \item In the presence of a nontrivial implicit prior corresponding to the specifics of the architecture of the method by which \pzpdf s are obtained, established methods cannot recover \nz; a principled hierarchical inference such as \Chippr\ is the only way to recover \nz\ from \pzpdf s. -\item %\Chippr, however, is sensitive to misspecification of the implicit prior; -Neither \Chippr\ nor traditional alternatives can recover \nz\ in the presence of a misspecified implicit prior; +\item Neither \Chippr\ nor traditional alternatives can recover \nz\ in the presence of a misspecified implicit prior; the implicit prior used to produce the \pzpdf\ catalog must be known and provided to \Chippr\ in order to recover the true \nz. \end{itemize} Given the significance of the implicit prior \citep{schmidt_evaluation_2020}, it is therefore imperative that those developing codes to obtain \pzpdf s provide a way to isolate the implicit prior and that those publishing \pzpdf\ catalogs provide the implicit prior to users. @@ -945,51 +851,29 @@ \section{Conclusion} In this case, it may not be possible to apply \Chippr\ without marginalizing over additional variables $\psi$ for the SEDs. In other words, obtaining the implicit prior from a template fitting code may be challenging or even require consideration of higher-dimensional PDFs such as $\pr{z, \mathrm{SED} \gvn \psi^{*}}$. -The situation appears more dire for data-driven techniques, whose training sets may not straightforwardly translate into an implicit prior. +The situation is more dire for data-driven techniques, whose training sets may not straightforwardly translate into a recoverable implicit prior. For example, some training set galaxies may contribute to the \pzpdf s more than others, resulting in different effective weights when factoring into, say, a histogram of training set redshifts as the implicit prior. Additionally, the weights may be stochastic, depending on the random seed used to initialize non-deterministic methods, precluding reproducibility. It is thus unclear whether the implicit prior can be meaningfully obtained from such methods at all. -%\aim{TODO: new paragraph for what can go wrong, -% call to community for what to do about it, -% what aspects of implicit prior are knowable and not knowable? -% if testable how so? outside the scope of paper, data producers be warned! -% a likelihood is better -- give us that if you can! -% focus on methods that acknowledge probabilistic structure of problem.} - -A thorough investigation of the degree to which the implicit prior can be meaningfully obtained is outside this paper but should be a priority for all consumers of \pzpdf s. -As an alternative, however, we must point out that if likelihoods were available rather than posteriors, the trouble with the implicit prior would be avoided altogether. -We thus encourage the community of those making \pzpdf s to consider developing such methods so that the resulting data products may be correctly used in scientific inference more generically. - -%\aim{TODO: Claim that implications for tomographic binning are severe, if that can be motivated by mean \Nz shifts.} -% by Section~\ref{sec:lsstdemo}.} - -%\subsection*{Recommendations for future work} - -%The following conclusions and recommendations can be made with confidence: +A thorough investigation of the degree to which the implicit prior can be meaningfully obtained is outside \new{}{the scope of} this paper but should be a priority for all consumers of \pzpdf s. +Alternatively, the trouble with the implicit prior would be avoided altogether if likelihoods were produced rather than posteriors. +We thus encourage the community of those making \pzpdf s to consider developing methods yielding likelihoods rather than posteriors so that the resulting data products may be correctly used in scientific inference more generically. -%\begin{enumerate} -% -% % \item The marginalized maximum likelihood estimator is an excellent estimator for strongly featured redshift distribution function with simple, clean photo-$z$ posteriors; stacking smooths features more than sampling and photo-$z$ point estimation. -% \item When the implicit prior is known to be a poor match to the data, only the results of \Chippr\ are satisfactory estimators of the redshift distribution function because they are the only methods that can account for the bias induced on the \pzpdf\ catalog by the method that produces it; this is the most compelling case for the sampler because of the ubiquity of inappropriate interim priors. -%\end{enumerate} - -By showing that \Chippr\ is effective in recovering the true redshift distribution function and posterior distributions on its parameters from catalogs of \pzpdf s, this work supports the production of \pzpdf s by upcoming photometric surveys such as \lsst\ to enable more accurate inference of the cosmological parameters. -We discourage researchers from co-adding \pzpdf s or converting them into point estimates of redshift and instead recommend the use of Bayesian probability to guide the usage of \pzpdf s. -We emphasize to those who produce \pzpdf s from data that it is essential to release the implicit prior used in generating this data product in order for any valid inference to be conducted by consumers of this information. +By showing that \Chippr\ is effective in recovering the true redshift distribution function and posterior distributions on its parameters from catalogs of \pzpdf s, this work supports the production of \pzpdf s by upcoming photometric surveys such as \lsst\ to enable more accurate inference of the cosmological parameters. +We discourage researchers from co-adding \pzpdf s or converting them into point estimates of redshift and instead recommend the use of Bayesian probability to guide the usage of \pzpdf s. +We emphasize to those who produce \pzpdf s from data that it is essential to \new{}{determine and} release the implicit prior used in generating this data product in order for any valid inference to be conducted by consumers of this information. Methodologies for obtaining \pzpdf s must therefore be designed such that there is a known implicit prior, i.e. one that is not implicit at all, so that likelihoods may be recovered. -The technique herein developed is applicable with minimal modification to other one-point statistics of redshift to which we will apply this method in the future, such as the redshift-dependent luminosity function and weak lensing mean distance ratio. +\new{}{While the \chippr\ prototype code is publicly available, follow-up work is currently under way to build an at-scale implementation of the \Chippr\ algorithm that will be applied to current and future data sets.} +The \new{}{\Chippr} technique herein developed is applicable with minimal modification to other one-point statistics of redshift to which we will apply this method in the future, such as the redshift-dependent luminosity function and weak lensing mean distance ratio. Future work will also include the extension of this fully probabilistic approach to higher-order statistics of redshift such as the two-point correlation function. \begin{acknowledgements} AIM acknowledges support from the Max Planck Society and the Alexander von Humboldt Foundation in the framework of the Max Planck-Humboldt Research Award endowed by the Federal Ministry of Education and Research. During the completion of this work, AIM was supported by National Science Foundation grant AST-1517237 and the U.S. Department of Energy, Office of Science, Office of Workforce Development for Teachers and Scientists, Office of Science Graduate Student Research (SCGSR) program, administered by the Oak Ridge Institute for Science and Education for the DOE under contract number DE‐SC0014664. The authors thank Phil Marshall for advice on relevant examples, Elisabeth Krause for assistance with the \cosmolike\ code, Mohammadjavad Vakili for statistical insights, Geoffrey Ryan for programming advice, and Boris Leistedt for other helpful comments in the development of \Chippr. -% The authors also acknowledge -% \aim{TODO: Circulate draft to GCCL; submit to ApJ; circulate to Dan Foreman-Mackey, Boris Leistedt, Kate Storey-Fisher; post to arXiv; circulate to Johann Cohen-Tanugi, Will Hartley, Alan Heavens, Mike Jarvis, Francois Lanusse, Ann Lee, Rachel Mandelbaum, Phil Marshall, Chris Morrison, Jeff Newman, Sam Schmidt, Anze Slosar, Josh Speagle, others for feedback.} This work was completed with generous nutritional support from the Center for Computational Astrophysics. -% \aim{TODO: Thank thesis readers.} \end{acknowledgements} %\aim{TODO: add software citation section.} @@ -999,20 +883,11 @@ \section{Conclusion} %\renewcommand{\thesubsection}{\Alph{subsection}} \numberwithin{equation}{section} -\section{Derivation} -\label{app:math} - -%We begin by parametrizing $N(z)$ in terms of $\vec{\theta}$, comprising some set of hyperparameters that define the form $N(z)$ may take in whatever basis we choose. -%We define a function $f_{\vec{\theta}}(z)=N(z)$ that transforms these hyperparameters into the redshift distribution function $N(z)$. -%Because -%\begin{equation} -%\eqlabel{eq:definition} -%N(z) \propto p(z \gvn \vec{\theta}), -%\end{equation} -%we may discontinue discussion of $N(z)$ in favor of the likelihood $p(z|\vec{\theta})$. +\section{}%Derivation} +%\label{app:math} -We perform the derivation of \Eq{eqn:fullpost} using log-probabilities. -What we wish to estimate is then the full log-posterior probability distribution (hereafter the full log-posterior) of the hyperparameters $\ndphi$ given the catalog of photometry $\{\data_{j}\}$. +We perform the derivation of \Eq{eqn:fullpost} using log-probabilities. +What we wish to estimate is \new{then}{} the full log-posterior probability distribution (hereafter the full log-posterior) of the hyperparameters $\ndphi$ given the catalog of photometry $\{\data_{j}\}$. By Bayes' Rule, the full log-posterior \begin{equation} @@ -1022,30 +897,30 @@ \section{Derivation} may be expressed in terms of the full log-likelihood probability distribution (hereafter the full log-likelihood) $\ln[\pr{\{\data_{j}\} \gvn \ndphi}]$ by way of a hyperprior log-probability distribution (hereafter the hyperprior) $\ln[\pr{\ndphi}]$ over the hyperparameters and the log-evidence probability of the data $\ln[\pr{\{\data_{j}\}}]$. However, the evidence is rarely known, so we probe the full log-posterior modulo an unknown constant of proportionality. -The full log-likelihood may be expanded in terms of a marginalization over the redshifts as parameters, as in +The full log-likelihood may be expanded in terms of a marginalization over the redshifts as parameters, as in \begin{equation} \label{eqn:marginalize} \ln[\pr{\{\data_{j}\} \gvn \ndphi}] = \ln\left[\integral{\pr{\{\data_{j}\} \gvn \{z_{j}\}} \pr{\{z_{j}\} \gvn \ndphi}}{\{z_{j}\}}\right]. \end{equation} -We shall make two assumptions of independence in order to make the problem tractable; their limitations are be discussed below. -First, we take $\ln[\pr{\{\data_{j}\} \gvn \{z_{j}\}}]$ to be the sum of $J$ individual log-likelihood distribution functions $\ln[\pr{\data_{j} \gvn z_{j}}]$, as in +We shall make two assumptions of independence in order to make the problem tractable; +their limitations are be discussed below. +First, we take $\ln[\pr{\{\data_{j}\} \gvn \{z_{j}\}}]$ to be the sum of $J$ individual log-likelihood distribution functions $\ln[\pr{\data_{j} \gvn z_{j}}]$, as in \begin{equation} \label{eqn:indiedat} \ln[\pr{\{\data_{j}\} \gvn \{z_{j}\}}] = \sum_{j=1}^{J}\ \ln[\pr{\data_{j} \gvn z_{j}}], \end{equation} a result of the definition of probabilistic independence encoded by the box in \Fig{fig:pgm}. -Second, we shall assume the true redshifts $\{z_{j}\}$ are $J$ independent draws from the true $\pr{z \gvn \ndphi}$. -Additionally, $J$ itself is a Poisson random variable. -The combination of these assumptions is given by +Second, we shall assume the true redshifts $\{z_{j}\}$ are $J$ independent draws from the true $\pr{z \gvn \ndphi}$. +Additionally, $J$ itself is a Poisson random variable. +The combination of these assumptions is given by \begin{equation} \label{eqn:indie} \ln[\pr{\{z_{j}\} \gvn \ndphi}] = -\integral{f(z; \ndphi)}{z} + \sum_{j=1}^{J}\ \ln[\pr{z_{j} \gvn \ndphi}]. \end{equation} -%It is important to note that the integral $\integral{n(z)}{z} N(z)\ dz$ is not constrained to equal the variable defining the Poisson distribution but instead $J$ by \Eq{eq:definition}, which can be thought of as another parameter. The derivation differs when $J$ is not known, say, when we want to learn about a distribution in nature rather than a distribution specific to data in hand, but for a photometric galaxy catalog where the desired quantity is $n(z)$ for the galaxies entering a larger cosmology calculation, it is a fixed quantity. -A detailed discussion of this matter may be found in \citet{foreman-mackey_exoplanet_2014}. -Applying Bayes' Rule, we may combine terms to obtain +A detailed discussion of this matter may be found in \citet{foreman-mackey_exoplanet_2014}. +Applying Bayes' Rule, we may combine terms to obtain \begin{align} \begin{split} \label{eqn:posterior} @@ -1053,16 +928,11 @@ \section{Derivation} \end{split} \end{align} -%\Eq{eq:posterior} contains two quantities that merit further discussion, the prior distribution $p(\vec{\theta})$ discussed further in \Sect{sec:exp} and the photo-$z$ log-likelihoods $\ln[p(\vec{d}_{j}|z_{j})]$ that have not been mentioned since \Eq{eq:marginalize}. -%Though photo-$z$ log-likelihoods would be desirable for use in these equations, they are not generally the product of either empirical and data-driven methods for obtaining photo-$z$ probability distributions. -%Though probabilistic photo-$z$s are typically reported as generic probability distributions $p(z_{j})$, the methods that produce them may be understood to always yield posteriors, probability distributions conditioned on the data we believe to be true. -%If they were not based in this assumption, they would require a sum over an infinite space of possible datasets. - Since we only have access to \pzip s, we must be able to write the full log-posterior in terms of log \pzip s rather than the log-likelihoods of \Eq{eqn:posterior}. -To do so, we will need an explicit statement of this implicit prior $\ndphi^{*}$ for whatever method is chosen to produce the \pzip s. +To do so, we will need an explicit statement of this implicit prior $\ndphi^{*}$ for whatever method is chosen to produce the \pzip s. -To perform the necessary transformation from likelihoods to posteriors, we follow the reasoning of \citet{foreman-mackey_exoplanet_2014}. -Let us consider the probability of the parameters conditioned on the data and an interim prior and rewrite the problematic likelihood of \Eq{eqn:posterior} as +To perform the necessary transformation from likelihoods to posteriors, we follow the reasoning of \citet{foreman-mackey_exoplanet_2014}. +Let us consider the probability of the parameters conditioned on the data and an interim prior and rewrite the problematic likelihood of \Eq{eqn:posterior} as \begin{align} \label{eqn:trick} \begin{split} @@ -1070,14 +940,14 @@ \section{Derivation} \end{split} \end{align} -Once the implicit prior $\ndphi^{*}$ is explicitly introduced, we may expand the last term in \Eq{eqn:trick} according to Bayes' Rule to get +Once the implicit prior $\ndphi^{*}$ is explicitly introduced, we may expand the last term in \Eq{eqn:trick} according to Bayes' Rule to get \begin{align} \begin{split} \label{eqn:expand} \ln[\pr{\data_{j} \gvn z}] = & \ln[\pr{\data_{j} \gvn z}] + \ln[\pr{z \gvn \data_{j}, \ndphi^{*}}] + \ln[\pr{\data_{j} \gvn \ndphi^{*}}] - \ln[\pr{z \gvn \ndphi^{*}}] - \ln[\pr{\data_{j} \gvn z, \ndphi^{*}}]. \end{split} \end{align} -Because there is no direct dependence of the data upon the hyperparameters, we may again expand the term $\ln[\pr{\data_{j} \gvn z, \ndphi^{*}}]$ to obtain +Because there is no direct dependence of the data upon the hyperparameters, we may again expand the term $\ln[\pr{\data_{j} \gvn z, \ndphi^{*}}]$ to obtain \begin{align} \begin{split} \label{eqn:indterm} @@ -1089,7 +959,7 @@ \section{Derivation} \label{eqn:cancel} \ln[\pr{\data_{j} \gvn z}] = \ln[\pr{z \gvn \data_{j}, \ndphi^{*}}] - \ln[\pr{z \gvn \ndphi^{*}}]. \end{equation} -We put this all together to get the full log-posterior probability distribution of +We put this all together to get the full log-posterior probability distribution of \begin{align} \begin{split} \label{eqn:final} @@ -1098,54 +968,42 @@ \section{Derivation} \end{align} which is equivalent to that of \citet{hogg_inferring_2010}, though the context differs. -The argument of the integral in the log-posterior of \Eq{eqn:final} depends solely on knowable quantities (and those we must explicitly assume) and can be calculated for a given sample of log \pzip s $\{\ln[\pr{z \gvn \data_{j}, \ndphi^{*}}]\}$ and the implicit prior $\pr{z \gvn \ndphi^{*}}$ with which they were obtained, noting the relation of +The argument of the integral in the log-posterior of \Eq{eqn:final} depends solely on knowable quantities (and those we must explicitly assume) and can be calculated for a given sample of log \pzip s $\{\ln[\pr{z \gvn \data_{j}, \ndphi^{*}}]\}$ and the implicit prior $\pr{z \gvn \ndphi^{*}}$ with which they were obtained, noting the relation of \begin{equation} \label{eqn:params} \pr{z \gvn \ndphi} = \frac{f(z; \ndphi)}{\integral{f(z; \ndphi)}{z}}. \end{equation} -Since we cannot know constant of proportionality, we sample the desired full log-posterior $\ln[\pr{\ndphi \gvn \{\data_{j}\}}]$ using Monte Carlo-Markov chain (MCMC) methods. - -% -%\begin{align} -%\begin{split} -%\eqlabel{eqn:fullpost} -%\ln[\pr{\ndphi \gvn \{\data_{j}\}}] & \propto \ln[\pr{\ndphi}] + \ln \left[\integral{\exp \left[\sum_{j=1}^{J} \left(\ln[\pr{z \gvn \data_{j}, \ndphi^{*}}] + \ln[\pr{z \gvn \ndphi}] - \ln[\pr{z \gvn \ndphi^{*}}]\right)\right]}{z}\right] , -%\end{split} -%\end{align} +Since we cannot know constant of proportionality, we sample the desired full log-posterior $\ln[\pr{\ndphi \gvn \{\data_{j}\}}]$ using Monte Carlo-Markov chain (MCMC) methods. %\section{Convergence Criteria} %\label{app:acorr} % -%\que{Cut convergence criteria section?} -% -%In addition to qualitative visual inspection of the chains, two quantities that probe the convergence of the sampler are used in this study, the autocorrelation time and the Gelman-Rubin convergence criterion. +%In addition to qualitative visual inspection of the chains, two quantities that probe the convergence of the sampler are used in this study, the autocorrelation time and the Gelman-Rubin convergence criterion. %%\Fig{fig:chains} shows the %evolution of the values of one parameter of one walker over the course of all %iterations of the sampler. % %%\begin{figure} %%%\includegraphics[width=0.5\textwidth]{figs/null/chain0.pdf} -%%\caption{This figure shows the evolution of one walker's parameter values for -%%one element of the parameter vector $\vec{\theta}$ as a function of iteration +%%\caption{This figure shows the evolution of one walker's parameter values for +%%one element of the parameter vector $\vec{\theta}$ as a function of iteration %%number, demonstrating the completion of the burn-in phase.} %%\label{fig:chains} %%\end{figure} % -%The autocorrelation time is effectively a measure of the efficiency of the method and can be described as the expected number of iterations necessary to accept a new sample independent of the current accepted sample. -%A sampler that converges faster will have a smaller autocorrelation time, and smaller autocorrelation times are preferable because it means fewer iterations are wasted on non-independent samples when independent samples are desired. -%See \citet{foreman-mackey_emcee_2013} for a more complete exploration of the autocorrelation time. -%In all tests discussed here, autocorrelation times across walkers and parameters were approximately 20, meaning two samples 20 or more iterations apart were independent, a satisfactory level of efficiency. -%Low autocorrelation times are a necessary but not always sufficient convergence condition, as the autocorrelation times calculated for tests in this paper were constant across all sub-runs, even those that were obviously burning in. +%The autocorrelation time is effectively a measure of the efficiency of the method and can be described as the expected number of iterations necessary to accept a new sample independent of the current accepted sample. +%A sampler that converges faster will have a smaller autocorrelation time, and smaller autocorrelation times are preferable because it means fewer iterations are wasted on non-independent samples when independent samples are desired. +%See \citet{foreman-mackey_emcee_2013} for a more complete exploration of the autocorrelation time. +%In all tests discussed here, autocorrelation times across walkers and parameters were approximately 20, meaning two samples 20 or more iterations apart were independent, a satisfactory level of efficiency. +%Low autocorrelation times are a necessary but not always sufficient convergence condition, as the autocorrelation times calculated for tests in this paper were constant across all sub-runs, even those that were obviously burning in. % %The Gelman-Rubin statistic %\begin{equation} %\label{eqn:gr} %R_{k} = \sqrt{\frac{(1 - \frac{2}{I_{0}}) w_{k} + \frac{2}{I_{0}} b_{k}}{w_{k}}}, %\end{equation} -%a weighted sum of the mean $w_{k}$ of the variances within individual walkers' chains and the variance $b_{k}$ between chains of different walkers $m$, is calculated over each sub-run $i$ to determine the duration of the burn-in period. -%Convergence is achieved when the statistic approaches unity. +%a weighted sum of the mean $w_{k}$ of the variances within individual walkers' chains and the variance $b_{k}$ between chains of different walkers $m$, is calculated over each sub-run $i$ to determine the duration of the burn-in period. +%Convergence is achieved when the statistic approaches unity. \bibliographystyle{apj} -\bibliography{draft} - -%\aim{TODO: find way to cite Dance Your Ph.D. video} +\bibliography{ms} \end{document} diff --git a/research/paper/figures/chippr/final_plot.png b/research/paper/figures/chippr/final_plot.png deleted file mode 100644 index 0d84000..0000000 Binary files a/research/paper/figures/chippr/final_plot.png and /dev/null differ diff --git a/research/paper/final_plot.png b/research/paper/final_plot.png new file mode 100644 index 0000000..d7e826a Binary files /dev/null and b/research/paper/final_plot.png differ diff --git a/research/paper/figures/chippr/flowchart.pdf b/research/paper/flowchart.pdf similarity index 100% rename from research/paper/figures/chippr/flowchart.pdf rename to research/paper/flowchart.pdf diff --git a/research/paper/figures/jain05.png b/research/paper/jain15.png similarity index 100% rename from research/paper/figures/jain05.png rename to research/paper/jain15.png diff --git a/research/paper/draft.bib b/research/paper/ms.bib similarity index 77% rename from research/paper/draft.bib rename to research/paper/ms.bib index c3426cc..e1ff053 100644 --- a/research/paper/draft.bib +++ b/research/paper/ms.bib @@ -4,6 +4,8 @@ @article{carrasco_kind_sparse_2014 volume = {441}, issn = {0035-8711}, shorttitle = {Sparse representation of photometric redshift probability density functions}, + url = {https://academic.oup.com/mnras/article/441/4/3550/1229381/Sparse-representation-of-photometric-redshift}, + doi = {10.1093/mnras/stu827}, number = {4}, journal = {Mon Not R Astron Soc}, author = {Carrasco Kind, Matias and Brunner, Robert J.}, @@ -16,6 +18,8 @@ @article{benitez_bayesian_2000 title = {Bayesian {Photometric} {Redshift} {Estimation}}, volume = {536}, issn = {0004-637X}, + url = {http://stacks.iop.org/0004-637X/536/i=2/a=571}, + doi = {10.1086/308947}, language = {en}, number = {2}, journal = {ApJ}, @@ -28,6 +32,8 @@ @article{krause_cosmolike_2017 title = {cosmolike {\textendash} cosmological likelihood analyses for photometric galaxy surveys}, volume = {470}, issn = {0035-8711}, + url = {https://academic.oup.com/mnras/article/470/2/2100/3850221/cosmolike-cosmological-likelihood-analyses-for}, + doi = {10.1093/mnras/stx1261}, number = {2}, journal = {Mon Not R Astron Soc}, author = {Krause, Elisabeth and Eifler, Tim}, @@ -38,6 +44,9 @@ @article{krause_cosmolike_2017 @article{abell_lsst_2009, title = {{LSST} {Science} {Book}, {Version} 2.0}, + volume = {0912.0201}, + url = {http://arxiv.org/abs/0912.0201}, + journal = {arXiv}, author = {Abell, Paul A. and Allison, Julius and Anderson, Scott F. and Andrew, John R. and Angel, J. Roger P. and Armus, Lee and Arnett, David and Asztalos, S. J. and Axelrod, Tim S. and Bailey, Stephen and Ballantyne, D. R. and Bankert, Justin R. and Barkhouse, Wayne A. and Barr, Jeffrey D. and Barrientos, L. Felipe and Barth, Aaron J. and Bartlett, James G. and Becker, Andrew C. and Becla, Jacek and Beers, Timothy C. and Bernstein, Joseph P. and Biswas, Rahul and Blanton, Michael R. and Bloom, Joshua S. and Bochanski, John J. and Boeshaar, Pat and Borne, Kirk D. and Bradac, Marusa and Brandt, W. N. and Bridge, Carrie R. and Brown, Michael E. and Brunner, Robert J. and Bullock, James S. and Burgasser, Adam J. and Burge, James H. and Burke, David L. and Cargile, Phillip A. and Chandrasekharan, Srinivasan and Chartas, George and Chesley, Steven R. and Chu, You-Hua and Cinabro, David and Claire, Mark W. and Claver, Charles F. and Clowe, Douglas and Connolly, A. J. and Cook, Kem H. and Cooke, Jeff and Cooray, Asantha and Covey, Kevin R. and Culliton, Christopher S. and de Jong, Roelof and de Vries, Willem H. and Debattista, Victor P. and Delgado, Francisco and Dell'Antonio, Ian P. and Dhital, Saurav and Di Stefano, Rosanne and Dickinson, Mark and Dilday, Benjamin and Djorgovski, S. G. and Dobler, Gregory and Donalek, Ciro and Dubois-Felsmann, Gregory and Durech, Josef and Eliasdottir, Ardis and Eracleous, Michael and Eyer, Laurent and Falco, Emilio E. and Fan, Xiaohui and Fassnacht, Christopher D. and Ferguson, Harry C. and Fernandez, Yanga R. and Fields, Brian D. and Finkbeiner, Douglas and Figueroa, Eduardo E. and Fox, Derek B. and Francke, Harold and Frank, James S. and Frieman, Josh and Fromenteau, Sebastien and Furqan, Muhammad and Galaz, Gaspar and Gal-Yam, A. and Garnavich, Peter and Gawiser, Eric and Geary, John and Gee, Perry and Gibson, Robert R. and Gilmore, Kirk and Grace, Emily A. and Green, Richard F. and Gressler, William J. and Grillmair, Carl J. and Habib, Salman and Haggerty, J. S. and Hamuy, Mario and Harris, Alan W. and Hawley, Suzanne L. and Heavens, Alan F. and Hebb, Leslie and Henry, Todd J. and Hileman, Edward and Hilton, Eric J. and Hoadley, Keri and Holberg, J. B. and Holman, Matt J. and Howell, Steve B. and Infante, Leopoldo and Ivezi{\'c}, {\v Z}eljko and Jacoby, Suzanne H. and Jain, Bhuvnesh and R and Jedicke and Jee, M. James and Jernigan, J. Garrett and Jha, Saurabh W. and Johnston, Kathryn V. and Jones, R. Lynne and Juric, Mario and Kaasalainen, Mikko and Styliani and Kafka and Kahn, Steven M. and Kaib, Nathan A. and Kalirai, Jason and Kantor, Jeff and Kasliwal, Mansi M. and Keeton, Charles R. and Kessler, Richard and Knezevic, Zoran and Kowalski, Adam and Krabbendam, Victor L. and Krughoff, K. Simon and Kulkarni, Shrinivas and Kuhlman, Stephen and Lacy, Mark and Lepine, Sebastien and Liang, Ming and Lien, Amy and Lira, Paulina and Long, Knox S. and Lorenz, Suzanne and Lotz, Jennifer M. and Lupton, R. H. and Lutz, Julie and Macri, Lucas M. and Mahabal, Ashish A. and Mandelbaum, Rachel and Marshall, Phil and May, Morgan and McGehee, Peregrine M. and Meadows, Brian T. and Meert, Alan and Milani, Andrea and Miller, Christopher J. and Miller, Michelle and Mills, David and Minniti, Dante and Monet, David and Mukadam, Anjum S. and Nakar, Ehud and Neill, Douglas R. and Newman, Jeffrey A. and Nikolaev, Sergei and Nordby, Martin and O'Connor, Paul and Oguri, Masamune and Oliver, John and Olivier, Scot S. and Olsen, Julia K. and Olsen, Knut and Olszewski, Edward W. and Oluseyi, Hakeem and Padilla, Nelson D. and Parker, Alex and Pepper, Joshua and Peterson, John R. and Petry, Catherine and Pinto, Philip A. and Pizagno, James L. and Popescu, Bogdan and Prsa, Andrej and Radcka, Veljko and Raddick, M. Jordan and Rasmussen, Andrew and Rau, Arne and Rho, Jeonghee and Rhoads, James E. and Richards, Gordon T. and Ridgway, Stephen T. and Robertson, Brant E. and Roskar, Rok and Saha, Abhijit and Sarajedini, Ata and Scannapieco, Evan and Schalk, Terry and Schindler, Rafe and Schmidt, Samuel and Schmidt, Sarah and Schneider, Donald P. and Schumacher, German and Scranton, Ryan and Sebag, Jacques and Seppala, Lynn G. and Shemmer, Ohad and Simon, Joshua D. and Sivertz, M. and Smith, Howard A. and Smith, J. Allyn and Smith, Nathan and Spitz, Anna H. and Stanford, Adam and Stassun, Keivan G. and Strader, Jay and Strauss, Michael A. and Stubbs, Christopher W. and Sweeney, Donald W. and Szalay, Alex and Szkody, Paula and Takada, Masahiro and Thorman, Paul and Trilling, David E. and Trimble, Virginia and Tyson, Anthony and Van Berg, Richard and Berk, Daniel Vanden and VanderPlas, Jake and Verde, Licia and Vrsnak, Bojan and Walkowicz, Lucianne M. and Wandelt, Benjamin D. and Wang, Sheng and Wang, Yun and Warner, Michael and Wechsler, Risa H. and West, Andrew A. and Wiecha, Oliver and Williams, Benjamin F. and Willman, Beth and Wittman, David and Wolff, Sidney C. and Wood-Vasey, W. Michael and Wozniak, Przemek and Young, Patrick and Zentner, Andrew and Zhan, Hu}, month = dec, year = {2009}, @@ -48,6 +57,8 @@ @article{sadeh_annz2:_2016 volume = {128}, issn = {1538-3873}, shorttitle = {{ANNz2}}, + url = {http://stacks.iop.org/1538-3873/128/i=968/a=104502}, + doi = {10.1088/1538-3873/128/968/104502}, language = {en}, number = {968}, journal = {PASP}, @@ -71,6 +82,8 @@ @article{hildebrandt_kids-450:_2017 volume = {465}, issn = {0035-8711}, shorttitle = {{KiDS}-450}, + url = {http://adsabs.harvard.edu/abs/2017MNRAS.465.1454H}, + doi = {10.1093/mnras/stw2805}, journal = {Mon Not R Astron Soc}, author = {Hildebrandt, H. and Viola, M. and Heymans, C. and Joudaki, S. and Kuijken, K. and Blake, C. and Erben, T. and Joachimi, B. and Klaes, D. and Miller, L. and Morrison, C. B. and Nakajima, R. and Verdoes Kleijn, G. and Amon, A. and Choi, A. and Covone, G. and de Jong, J. T. A. and Dvornik, A. and Fenech Conti, I. and Grado, A. and Harnois-D{\'e}raps, J. and Herbonnet, R. and Hoekstra, H. and K{\"o}hlinger, F. and McFarland, J. and Mead, A. and Merten, J. and Napolitano, N. and Peacock, J. A. and Radovich, M. and Schneider, P. and Simon, P. and Valentijn, E. A. and van den Busch, J. L. and van Uitert, E. and Van Waerbeke, L.}, month = feb, @@ -83,6 +96,8 @@ @article{herbel_redshift_2017 volume = {2017}, issn = {1475-7516}, shorttitle = {The redshift distribution of cosmological samples}, + url = {http://stacks.iop.org/1475-7516/2017/i=08/a=035}, + doi = {10.1088/1475-7516/2017/08/035}, language = {en}, number = {08}, journal = {J. Cosmol. Astropart. Phys.}, @@ -95,21 +110,26 @@ @article{coe_fisher_2009 title = {Fisher {Matrices} and {Confidence} {Ellipses}: {A} {Quick}-{Start} {Guide} and {Software}}, volume = {0906.4123}, shorttitle = {Fisher {Matrices} and {Confidence} {Ellipses}}, + url = {http://arxiv.org/abs/0906.4123}, journal = {arXiv}, author = {Coe, Dan}, month = jun, year = {2009}, + note = {arXiv: 0906.4123}, } @article{hearin_general_2012, title = {General {Requirements} on {Matter} {Power} {Spectrum} {Predictions} for {Cosmology} with {Weak} {Lensing} {Tomography}}, volume = {2012}, issn = {1475-7516}, + url = {http://arxiv.org/abs/1111.0052}, + doi = {10.1088/1475-7516/2012/04/034}, number = {04}, journal = {Journal of Cosmology and Astroparticle Physics}, author = {Hearin, Andrew P. and Zentner, Andrew R. and Ma, Zhaoming}, month = apr, year = {2012}, + note = {arXiv: 1111.0052}, pages = {034--034}, } @@ -117,6 +137,8 @@ @article{schmidt_improved_2013 title = {Improved photometric redshifts via enhanced estimates of system response, galaxy templates and magnitude priors}, volume = {431}, issn = {0035-8711}, + url = {https://academic.oup.com/mnras/article/431/3/2766/1080084}, + doi = {10.1093/mnras/stt373}, language = {en}, number = {3}, journal = {Mon Not R Astron Soc}, @@ -129,6 +151,7 @@ @article{schmidt_improved_2013 @article{mandelbaum_weak_2017, title = {Weak lensing for precision cosmology}, volume = {1710.03235}, + url = {http://arxiv.org/abs/1710.03235}, journal = {arXiv}, author = {Mandelbaum, Rachel}, month = oct, @@ -139,6 +162,8 @@ @article{tanaka_photometric_2018 title = {Photometric redshifts for {Hyper} {Suprime}-{Cam} {Subaru} {Strategic} {Program} {Data} {Release} 1}, volume = {70}, issn = {0004-6264}, + url = {https://academic.oup.com/pasj/article/70/SP1/S9/4494086}, + doi = {10.1093/pasj/psx077}, language = {en}, number = {SP1}, journal = {Publ Astron Soc Jpn Nihon Tenmon Gakkai}, @@ -152,8 +177,10 @@ @article{asorey_galaxy_2016 title = {Galaxy clustering with photometric surveys using {PDF} redshift information}, volume = {459}, issn = {0035-8711, 1365-2966}, + url = {http://arxiv.org/abs/1601.00357}, + doi = {10.1093/mnras/stw721}, number = {2}, - journal = {Monthly Notices of the Royal Astronomical Society}, + journal = {Mon Not R Astron Soc}, author = {Asorey, J. and Kind, M. Carrasco and Sevilla-Noarbe, I. and Brunner, R. J. and Thaler, J.}, month = jun, year = {2016}, @@ -164,6 +191,8 @@ @article{sheldon_photometric_2012 title = {Photometric {Redshift} {Probability} {Distributions} for {Galaxies} in the {SDSS} {DR8}}, volume = {201}, issn = {0067-0049}, + url = {http://stacks.iop.org/0067-0049/201/i=2/a=32}, + doi = {10.1088/0067-0049/201/2/32}, language = {en}, number = {2}, journal = {ApJS}, @@ -176,6 +205,8 @@ @article{mandelbaum_precision_2008 title = {Precision photometric redshift calibration for galaxy{\textendash}galaxy weak lensing}, volume = {386}, issn = {0035-8711}, + url = {https://academic.oup.com/mnras/article/386/2/781/1056461}, + doi = {10.1111/j.1365-2966.2008.12947.x}, language = {en}, number = {2}, journal = {Mon Not R Astron Soc}, @@ -187,16 +218,20 @@ @article{mandelbaum_precision_2008 @article{laureijs_euclid_2011, title = {Euclid {Definition} {Study} {Report}}, + url = {http://arxiv.org/abs/1110.3193}, journal = {ESA/SRE(2011)12}, author = {Laureijs, R. and Amiaux, J. and Arduini, S. and Augu{\`e}res, J.-L. and Brinchmann, J. and Cole, R. and Cropper, M. and Dabin, C. and Duvet, L. and Ealet, A. and Garilli, B. and Gondoin, P. and Guzzo, L. and Hoar, J. and Hoekstra, H. and Holmes, R. and Kitching, T. and Maciaszek, T. and Mellier, Y. and Pasian, F. and Percival, W. and Rhodes, J. and Criado, G. Saavedra and Sauvage, M. and Scaramella, R. and Valenziano, L. and Warren, S. and Bender, R. and Castander, F. and Cimatti, A. and F{\`e}vre, O. Le and Kurki-Suonio, H. and Levi, M. and Lilje, P. and Meylan, G. and Nichol, R. and Pedersen, K. and Popa, V. and Lopez, R. Rebolo and Rix, H.-W. and Rottgering, H. and Zeilinger, W. and Grupp, F. and Hudelot, P. and Massey, R. and Meneghetti, M. and Miller, L. and Paltani, S. and Paulin-Henriksson, S. and Pires, S. and Saxton, C. and Schrabback, T. and Seidel, G. and Walsh, J. and Aghanim, N. and Amendola, L. and Bartlett, J. and Baccigalupi, C. and Beaulieu, J.-P. and Benabed, K. and Cuby, J.-G. and Elbaz, D. and Fosalba, P. and Gavazzi, G. and Helmi, A. and Hook, I. and Irwin, M. and Kneib, J.-P. and Kunz, M. and Mannucci, F. and Moscardini, L. and Tao, C. and Teyssier, R. and Weller, J. and Zamorani, G. and Osorio, M. R. Zapatero and Boulade, O. and Foumond, J. J. and Di Giorgio, A. and Guttridge, P. and James, A. and Kemp, M. and Martignac, J. and Spencer, A. and Walton, D. and Bl{\"u}mchen, T. and Bonoli, C. and Bortoletto, F. and Cerna, C. and Corcione, L. and Fabron, C. and Jahnke, K. and Ligori, S. and Madrid, F. and Martin, L. and Morgante, G. and Pamplona, T. and Prieto, E. and Riva, M. and Toledo, R. and Trifoglio, M. and Zerbi, F. and Abdalla, F. and Douspis, M. and Grenet, C. and Borgani, S. and Bouwens, R. and Courbin, F. and Delouis, J.-M. and Dubath, P. and Fontana, A. and Frailis, M. and Grazian, A. and Koppenh{\"o}fer, J. and Mansutti, O. and Melchior, M. and Mignoli, M. and Mohr, J. and Neissner, C. and Noddle, K. and Poncet, M. and Scodeggio, M. and Serrano, S. and Shane, N. and Starck, J.-L. and Surace, C. and Taylor, A. and Verdoes-Kleijn, G. and Vuerli, C. and Williams, O. R. and Zacchei, A. and Altieri, B. and Sanz, I. Escudero and Kohley, R. and Oosterbroek, T. and Astier, P. and Bacon, D. and Bardelli, S. and Baugh, C. and Bellagamba, F. and Benoist, C. and Bianchi, D. and Biviano, A. and Branchini, E. and Carbone, C. and Cardone, V. and Clements, D. and Colombi, S. and Conselice, C. and Cresci, G. and Deacon, N. and Dunlop, J. and Fedeli, C. and Fontanot, F. and Franzetti, P. and Giocoli, C. and Garcia-Bellido, J. and Gow, J. and Heavens, A. and Hewett, P. and Heymans, C. and Holland, A. and Huang, Z. and Ilbert, O. and Joachimi, B. and Jennins, E. and Kerins, E. and Kiessling, A. and Kirk, D. and Kotak, R. and Krause, O. and Lahav, O. and van Leeuwen, F. and Lesgourgues, J. and Lombardi, M. and Magliocchetti, M. and Maguire, K. and Majerotto, E. and Maoli, R. and Marulli, F. and Maurogordato, S. and McCracken, H. and McLure, R. and Melchiorri, A. and Merson, A. and Moresco, M. and Nonino, M. and Norberg, P. and Peacock, J. and Pello, R. and Penny, M. and Pettorino, V. and Di Porto, C. and Pozzetti, L. and Quercellini, C. and Radovich, M. and Rassat, A. and Roche, N. and Ronayette, S. and Rossetti, E. and Sartoris, B. and Schneider, P. and Semboloni, E. and Serjeant, S. and Simpson, F. and Skordis, C. and Smadja, G. and Smartt, S. and Spano, P. and Spiro, S. and Sullivan, M. and Tilquin, A. and Trotta, R. and Verde, L. and Wang, Y. and Williger, G. and Zhao, G. and Zoubian, J. and Zucca, E.}, month = oct, year = {2011}, + note = {arXiv:1110.3193}, } @article{malz_approximating_2018, title = {Approximating {Photo}- z {PDFs} for {Large} {Surveys}}, volume = {156}, issn = {1538-3881}, + url = {http://stacks.iop.org/1538-3881/156/i=1/a=35}, + doi = {10.3847/1538-3881/aac6b5}, language = {en}, number = {1}, journal = {AJ}, @@ -209,6 +244,8 @@ @article{leung_bayesian_2017 title = {Bayesian {Redshift} {Classification} of {Emission}-line {Galaxies} with {Photometric} {Equivalent} {Widths}}, volume = {843}, issn = {0004-637X}, + url = {http://stacks.iop.org/0004-637X/843/i=2/a=130}, + doi = {10.3847/1538-4357/aa71af}, language = {en}, number = {2}, journal = {ApJ}, @@ -222,6 +259,8 @@ @article{masters_mapping_2015 volume = {813}, issn = {0004-637X}, shorttitle = {Mapping the {Galaxy} {Color}{\textendash}{Redshift} {Relation}}, + url = {http://stacks.iop.org/0004-637X/813/i=1/a=53}, + doi = {10.1088/0004-637X/813/1/53}, language = {en}, number = {1}, journal = {ApJ}, @@ -234,20 +273,24 @@ @article{ball_robust_2008 title = {Robust {Machine} {Learning} {Applied} to {Astronomical} {Data} {Sets}. {III}. {Probabilistic} {Photometric} {Redshifts} for {Galaxies} and {Quasars} in the {SDSS} and {GALEX}}, volume = {683}, issn = {0004-637X}, + url = {http://iopscience.iop.org/article/10.1086/589646}, + doi = {10.1086/589646}, language = {en}, number = {1}, - journal = {The Astrophysical Journal}, + journal = {ApJ}, author = {Ball, Nicholas M. and Brunner, Robert J. and Myers, Adam D. and Strand, Natalie E. and Alberts, Stacey L. and Tcheng, David}, month = aug, year = {2008}, pages = {12--21} } -@article{baum_photoelectric_1962, +@incollection{baum_photoelectric_1962, series = {Problems of {Extra}-{Galactic} {Research}}, title = {Photoelectric {Magnitudes} and {Red}-{Shifts}}, volume = {15}, - journal = {Proceedings from IAU Symposium}, + url = {http://adsabs.harvard.edu/abs/1962IAUS...15..390B}, + booktitle = {Proceedings from {IAU} {Symposium}}, + publisher = {Macmillan Press}, author = {Baum, W. A.}, year = {1962}, pages = {390} @@ -257,8 +300,10 @@ @article{bonnett_using_2015 title = {Using neural networks to estimate redshift distributions. {An} application to {CFHTLenS}}, volume = {449}, issn = {0035-8711}, + url = {http://mnras.oxfordjournals.org/content/449/1/1043}, + doi = {10.1093/mnras/stv230}, number = {1}, - journal = {Monthly Notices of the Royal Astronomical Society}, + journal = {Mon Not R Astron Soc}, author = {Bonnett, C.}, month = mar, year = {2015}, @@ -269,9 +314,11 @@ @article{budavari_unified_2009 title = {A {UNIFIED} {FRAMEWORK} {FOR} {PHOTOMETRIC} {REDSHIFTS}}, volume = {695}, issn = {0004-637X}, + url = {http://iopscience.iop.org/article/10.1088/0004-637X/695/1/747}, + doi = {10.1088/0004-637X/695/1/747}, language = {en}, number = {1}, - journal = {The Astrophysical Journal}, + journal = {ApJ}, author = {Budav{\'a}ri, Tam{\'a}s}, month = apr, year = {2009}, @@ -282,9 +329,11 @@ @article{carliles_random_2010 title = {{RANDOM} {FORESTS} {FOR} {PHOTOMETRIC} {REDSHIFTS}}, volume = {712}, issn = {0004-637X}, + url = {http://iopscience.iop.org/article/10.1088/0004-637X/712/1/511}, + doi = {10.1088/0004-637X/712/1/511}, language = {en}, number = {1}, - journal = {The Astrophysical Journal}, + journal = {ApJ}, author = {Carliles, Samuel and Budav{\'a}ri, Tam{\'a}s and Heinis, S{\'e}bastien and Priebe, Carey and Szalay, Alexander S.}, month = mar, year = {2010}, @@ -295,8 +344,10 @@ @article{carrasco_kind_tpz:_2013 title = {{TPZ}: photometric redshift {PDFs} and ancillary information by using prediction trees and random forests}, volume = {432}, issn = {0035-8711}, + url = {http://mnras.oxfordjournals.org/content/432/2/1483}, + doi = {10.1093/mnras/stt574}, number = {2}, - journal = {Monthly Notices of the Royal Astronomical Society}, + journal = {Mon Not R Astron Soc}, author = {Carrasco Kind, M. and Brunner, R. J.}, month = may, year = {2013}, @@ -307,8 +358,10 @@ @article{carrasco_kind_exhausting_2014 title = {Exhausting the information: novel {Bayesian} combination of photometric redshift {PDFs}}, volume = {442}, issn = {0035-8711}, + url = {http://mnras.oxfordjournals.org/content/442/4/3380}, + doi = {10.1093/mnras/stu1098}, number = {4}, - journal = {Monthly Notices of the Royal Astronomical Society}, + journal = {Mon Not R Astron Soc}, author = {Carrasco Kind, M. and Brunner, R. J.}, month = jun, year = {2014}, @@ -319,9 +372,11 @@ @article{dipompeo_quasar_2015 title = {Quasar probabilities and redshifts from {WISE} mid-{IR} through {GALEX} {UV} photometry}, volume = {452}, issn = {0035-8711}, + url = {http://mnras.oxfordjournals.org/content/452/3/3124.full}, + doi = {10.1093/mnras/stv1562}, language = {en}, number = {3}, - journal = {Monthly Notices of the Royal Astronomical Society}, + journal = {Mon Not R Astron Soc}, author = {DiPompeo, M. A. and Bovy, J. and Myers, A. D. and Lang, D.}, month = jul, year = {2015}, @@ -332,8 +387,10 @@ @article{foreman-mackey_emcee_2013 title = {emcee : {The} {MCMC} {Hammer}}, volume = {125}, issn = {00046280}, + url = {http://arxiv.org/abs/1202.3665}, + doi = {10.1086/670067}, number = {925}, - journal = {Publications of the Astronomical Society of the Pacific}, + journal = {PASP}, author = {Foreman-Mackey, Daniel and Hogg, David W. and Lang, Dustin and Goodman, Jonathan}, month = mar, year = {2013}, @@ -344,9 +401,11 @@ @article{foreman-mackey_exoplanet_2014 title = {{EXOPLANET} {POPULATION} {INFERENCE} {AND} {THE} {ABUNDANCE} {OF} {EARTH} {ANALOGS} {FROM} {NOISY}, {INCOMPLETE} {CATALOGS}}, volume = {795}, issn = {1538-4357}, + url = {http://iopscience.iop.org/article/10.1088/0004-637X/795/1/64}, + doi = {10.1088/0004-637X/795/1/64}, language = {en}, number = {1}, - journal = {The Astrophysical Journal}, + journal = {ApJ}, author = {Foreman-Mackey, Daniel and Hogg, David W. and Morton, Timothy D.}, month = oct, year = {2014}, @@ -357,7 +416,9 @@ @article{gorecki_new_2014 title = {A new method to improve photometric redshift reconstruction}, volume = {561}, issn = {0004-6361}, - journal = {Astronomy \& Astrophysics}, + url = {http://arxiv.org/abs/1301.3010}, + doi = {10.1051/0004-6361/201321102}, + journal = {A\&A}, author = {Gorecki, Alexia and Abate, Alexandra and Ansari, R{\'e}za and Barrau, Aur{\'e}lien and Baumont, Sylvain and Moniez, Marc and Ricol, Jean-St{\'e}phane}, month = jan, year = {2014}, @@ -367,6 +428,7 @@ @article{gorecki_new_2014 @article{hogg_data_2012, title = {Data analysis recipes: {Probability} calculus for inference}, volume = {1205.4446}, + url = {http://arxiv.org/abs/1205.4446}, journal = {arXiv}, author = {Hogg, David W.}, month = may, @@ -375,6 +437,7 @@ @article{hogg_data_2012 @incollection{koo_photometric_1999, title = {Photometric {Redshifts}: {A} {Perspective} from an {Old}-{Timer} on {Its} {Past}, {Present}, and {Potential}}, + url = {http://arxiv.org/abs/astro-ph/9907273}, booktitle = {Photom. {Redshifts} {High} {Redshift} {Galaxies}}, publisher = {ASP Conference Series}, author = {Koo, David C.}, @@ -387,8 +450,10 @@ @article{lima_estimating_2008 title = {Estimating the redshift distribution of photometric galaxy samples}, volume = {390}, issn = {00358711}, + url = {http://mnras.oxfordjournals.org/content/390/1/118}, + doi = {10.1111/j.1365-2966.2008.13510.x}, number = {1}, - journal = {Monthly Notices of the Royal Astronomical Society}, + journal = {Mon Not R Astron Soc}, author = {Lima, Marcos and Cunha, Carlos E. and Oyaizu, Hiroaki and Frieman, Joshua and Lin, Huan and Sheldon, Erin S.}, month = oct, year = {2008}, @@ -398,18 +463,21 @@ @article{lima_estimating_2008 @article{menard_clustering-based_2013, title = {Clustering-based redshift estimation: method and application to data}, volume = {1303.4722}, + url = {http://arxiv.org/abs/1303.4722}, journal = {arXiv}, author = {M{\'e}nard, Brice and Scranton, Ryan and Schmidt, Samuel and Morrison, Chris and Jeong, Donghui and Budavari, Tamas and Rahman, Mubdi}, month = mar, - year = {2013}, + year = {2013} } @article{norberg_2df_2002, title = {The {2dF} {Galaxy} {Redshift} {Survey}: the {bJ}-band galaxy luminosity function and survey selection function}, volume = {336}, issn = {0035-8711}, + url = {http://mnras.oxfordjournals.org/content/336/3/907}, + doi = {10.1046/j.1365-8711.2002.05831.x}, number = {3}, - journal = {Monthly Notices of the Royal Astronomical Society}, + journal = {Mon Not R Astron Soc}, author = {Norberg, P. and Cole, S. and Baugh, C. M. and Frenk, C. S. and Baldry, I. and Bland-Hawthorn, J. and Bridges, T. and Cannon, R. and Colless, M. and Collins, C. and Couch, W. and Cross, N. J. G. and Dalton, G. and De Propris, R. and Driver, S. P. and Efstathiou, G. and Ellis, R. S. and Glazebrook, K. and Jackson, C. and Lahav, O. and Lewis, I. and Lumsden, S. and Maddox, S. and Madgwick, D. and Peacock, J. A. and Peterson, B. A. and Sutherland, W. and Taylor, K.}, month = nov, year = {2002}, @@ -420,8 +488,10 @@ @article{sanchez_clustering_2013 title = {The clustering of galaxies in the {SDSS}-{III} {Baryon} {Oscillation} {Spectroscopic} {Survey}: cosmological constraints from the full shape of the clustering wedges}, volume = {433}, issn = {0035-8711}, + url = {http://mnras.oxfordjournals.org/content/433/2/1202}, + doi = {10.1093/mnras/stt799}, number = {2}, - journal = {Mon. Not. R. Astron. Soc.}, + journal = {Mon Not R Astron Soc}, author = {Sanchez, A. G. and Kazin, E. A. and Beutler, F. and Chuang, C.-H. and Cuesta, A. J. and Eisenstein, D. J. and Manera, M. and Montesano, F. and Nichol, R. C. and Padmanabhan, N. and Percival, W. and Prada, F. and Ross, A. J. and Schlegel, D. J. and Tinker, J. and Tojeiro, R. and Weinberg, D. H. and Xu, X. and Brinkmann, J. and Brownstein, J. R. and Schneider, D. P. and Thomas, D.}, month = jun, year = {2013}, @@ -432,8 +502,10 @@ @article{van_breukelen_reliable_2009 title = {A reliable cluster detection technique using photometric redshifts: introducing the {2TecX} algorithm}, volume = {395}, issn = {00358711}, + url = {http://mnras.oxfordjournals.org/content/395/4/1845}, + doi = {10.1111/j.1365-2966.2009.14692.x}, number = {4}, - journal = {Monthly Notices of the Royal Astronomical Society}, + journal = {Mon Not R Astron Soc}, author = {van Breukelen, Caroline and Clewley, Lee}, month = jun, year = {2009}, @@ -444,8 +516,10 @@ @article{viironen_high_2015 title = {High redshift galaxies in the {ALHAMBRA} survey}, volume = {576}, issn = {0004-6361}, + url = {http://www.aanda.org/articles/aa/full_html/2015/04/aa25382-14/aa25382-14.html}, + doi = {10.1051/0004-6361/201425382}, language = {en}, - journal = {Astronomy \& Astrophysics}, + journal = {A\&A}, author = {Viironen, K. and Mar{\'i}n-Franch, A. and L{\'o}pez-Sanjuan, C. and Varela, J. and Chaves-Montero, J. and Crist{\'o}bal-Hornillos, D. and Molino, A. and Fern{\'a}ndez-Soto, A. and Vilella-Rojo, G. and Ascaso, B. and Cenarro, A. J. and Cervi{\~n}o, M. and Cepa, J. and Ederoclite, A. and M{\'a}rquez, I. and Masegosa, J. and Moles, M. and Oteo, I. and Povi{\'c}, M. and Aguerri, J. A. L. and Alfaro, E. and Aparicio-Villegas, T. and Ben{\'i}tez, N. and Broadhurst, T. and Cabrera-Ca{\~n}o, J. and Castander, J. F. and Del Olmo, A. and Gonz{\'a}lez Delgado, R. M. and Husillos, C. and Infante, L. and Mart{\'i}nez, V. J. and Perea, J. and Prada, F. and Quintana, J. M.}, month = mar, year = {2015}, @@ -455,7 +529,8 @@ @article{viironen_high_2015 @article{hildebrandt_phat:_2010, title = {{PHAT}: {PHoto}-z {Accuracy} {Testing}}, volume = {523}, - journal = {Astronomy \& Astrophysics}, + doi = {10.1051/0004-6361/201014885}, + journal = {A\&A}, author = {Hildebrandt, H. and Arnouts, S. and Capak, P. and Moustakas, L. A. and Wolf, C. and Abdalla, F. B. and Assef, R. J. and Banerji, M. and Ben{\'i}tez, N. and Brammer, G. B. and Budav{\'a}ri, T. and Carliles, S. and Coe, D. and Dahlen, T. and Feldmann, R. and Gerdes, D. and Gillis, B. and Ilbert, O. and Kotulla, R. and Lahav, O. and Li, I. H. and Miralles, J.-M. and Purger, N. and Schmidt, S. and Singal, J.}, month = nov, year = {2010}, @@ -465,6 +540,7 @@ @article{hildebrandt_phat:_2010 @article{leistedt_data-driven_2017, title = {Data-driven, {Interpretable} {Photometric} {Redshifts} {Trained} on {Heterogeneous} and {Unrepresentative} {Data}}, volume = {838}, + doi = {10.3847/1538-4357/aa6332}, journal = {The Astrophysical Journal}, author = {Leistedt, B. and Hogg, D. W.}, month = mar, @@ -476,6 +552,8 @@ @article{leistedt_hierarchical_2016 title = {Hierarchical {Bayesian} inference of galaxy redshift distributions from photometric surveys}, volume = {460}, issn = {0035-8711}, + url = {https://academic.oup.com/mnras/article/460/4/4258/2609193}, + doi = {10.1093/mnras/stw1304}, language = {en}, number = {4}, journal = {Mon Not R Astron Soc}, @@ -485,10 +563,22 @@ @article{leistedt_hierarchical_2016 pages = {4258--4267}, } +@misc{malz_probabilistic_2019, + type = {Dance {Video}}, + title = {Probabilistic methods for cosmological analysis with uncertainty-dominated data}, + url = {https://www.youtube.com/watch?v=vKs3PYqZWg8}, + publisher = {YouTube}, + author = {Malz, Alex I.}, + month = jan, + year = {2019}, + note = {https://youtu.be/vKs3PYqZWg8}, +} + @article{jain_whole_2015, title = {The {Whole} is {Greater} than the {Sum} of the {Parts}: {Optimizing} the {Joint} {Science} {Return} from {LSST}, {Euclid} and {WFIRST}}, volume = {1501.07897}, shorttitle = {The {Whole} is {Greater} than the {Sum} of the {Parts}}, + url = {http://arxiv.org/abs/1501.07897}, journal = {arXiv}, author = {Jain, B. and Spergel, D. and Bean, R. and Connolly, A. and Dell'antonio, I. and Frieman, J. and Gawiser, E. and Gehrels, N. and Gladney, L. and Heitmann, K. and Helou, G. and Hirata, C. and Ho, S. and Ivezi{\'c}, {\v Z} and Jarvis, M. and Kahn, S. and Kalirai, J. and Kim, A. and Lupton, R. and Mandelbaum, R. and Marshall, P. and Newman, J. A. and Perlmutter, S. and Postman, M. and Rhodes, J. and Strauss, M. A. and Tyson, J. A. and Walkowicz, L. and Wood-Vasey, W. M.}, month = jan, @@ -498,14 +588,18 @@ @article{jain_whole_2015 @misc{rohatgi_webplotdigitizer_2019, address = {San Francisco, California, USA}, title = {{WebPlotDigitizer}}, + url = {https://automeris.io/WebPlotDigitizer}, author = {Rohatgi, Ankit}, year = {2019}, + note = {https://automeris.io/WebPlotDigitizer}, } @article{yang_calibrating_2018, title = {Calibrating magnification bias for the {EG} statistic to test general relativity}, volume = {481}, issn = {0035-8711}, + url = {https://academic.oup.com/mnras/article/481/2/1441/5092610}, + doi = {10.1093/mnras/sty2353}, language = {en}, number = {2}, journal = {Mon Not R Astron Soc}, @@ -519,6 +613,8 @@ @article{abruzzo_impact_2019 title = {The impact of photometric redshift errors on lensing statistics in ray-tracing simulations}, volume = {486}, issn = {0035-8711}, + url = {https://academic.oup.com/mnras/article/486/2/2730/5475123}, + doi = {10.1093/mnras/stz1016}, language = {en}, number = {2}, journal = {Mon Not R Astron Soc}, @@ -533,6 +629,8 @@ @article{hildebrandt_cfhtlens:_2012 volume = {421}, issn = {0035-8711}, shorttitle = {{CFHTLenS}}, + url = {https://academic.oup.com/mnras/article/421/3/2355/1078193}, + doi = {10.1111/j.1365-2966.2012.20468.x}, language = {en}, number = {3}, journal = {Mon Not R Astron Soc}, @@ -547,6 +645,8 @@ @article{benjamin_cfhtlens_2013 volume = {431}, issn = {0035-8711}, shorttitle = {{CFHTLenS} tomographic weak lensing}, + url = {https://academic.oup.com/mnras/article/431/2/1547/1459040}, + doi = {10.1093/mnras/stt276}, language = {en}, number = {2}, journal = {Mon Not R Astron Soc}, @@ -560,6 +660,8 @@ @article{kelly_weighing_2014 title = {Weighing the {Giants} {\textendash} {II}. {Improved} calibration of photometry from stellar colours and accurate photometric redshifts}, volume = {439}, issn = {0035-8711}, + url = {https://academic.oup.com/mnras/article/439/1/28/963026}, + doi = {10.1093/mnras/stt1946}, language = {en}, number = {1}, journal = {Mon Not R Astron Soc}, @@ -573,6 +675,8 @@ @article{ma_effects_2006 title = {Effects of {Photometric} {Redshift} {Uncertainties} on {Weak}-{Lensing} {Tomography}}, volume = {636}, issn = {0004-637X}, + url = {https://iopscience.iop.org/article/10.1086/497068/meta}, + doi = {10.1086/497068}, language = {en}, number = {1}, journal = {ApJ}, @@ -586,6 +690,8 @@ @article{leistedt_hierarchical_2019 title = {Hierarchical {Modeling} and {Statistical} {Calibration} for {Photometric} {Redshifts}}, volume = {881}, issn = {0004-637X}, + url = {https://doi.org/10.3847%2F1538-4357%2Fab2d29}, + doi = {10.3847/1538-4357/ab2d29}, language = {en}, number = {1}, journal = {ApJ}, @@ -595,20 +701,13 @@ @article{leistedt_hierarchical_2019 pages = {80}, } -@article{awan_angular_2019, - title = {Angular {Correlation} {Function} {Estimators} {Accounting} for {Contamination} from {Probabilistic} {Distance} {Measurements}}, - volume = {1911.07832}, - journal = {arXiv}, - author = {Awan, Humna and Gawiser, Eric}, - month = dec, - year = {2019}, -} - @article{hoyle_dark_2018, title = {Dark {Energy} {Survey} {Year} 1 {Results}: redshift distributions of the weak-lensing source galaxies}, volume = {478}, issn = {0035-8711}, shorttitle = {Dark {Energy} {Survey} {Year} 1 {Results}}, + url = {https://academic.oup.com/mnras/article/478/1/592/4975790}, + doi = {10.1093/mnras/sty957}, language = {en}, number = {1}, journal = {Mon Not R Astron Soc}, @@ -621,32 +720,27 @@ @article{hoyle_dark_2018 @article{dahlen_critical_2013, title = {A {Critical} {Assessment} of {Photometric} {Redshift} {Methods}: {A} {CANDELS} {Investigation}}, volume = {775}, - journal = {The Astrophysical Journal}, + doi = {10.1088/0004-637X/775/2/93}, + journal = {ApJ}, author = {Dahlen, T. and Mobasher, B. and Faber, S. M. and Ferguson, H. C. and Barro, G. and Finkelstein, S. L. and Finlator, K. and Fontana, A. and Gruetzbauch, R. and Johnson, S. and Pforr, J. and Salvato, M. and Wiklind, T. and Wuyts, S. and Acquaviva, V. and Dickinson, M. E. and Guo, Y. and Huang, J. and Huang, K.-H. and Newman, J. A. and Bell, E. F. and Conselice, C. J. and Galametz, A. and Gawiser, E. and Giavalisco, M. and Grogin, N. A. and Hathi, N. and Kocevski, D. and Koekemoer, A. M. and Koo, D. C. and Lee, K.-S. and McGrath, E. J. and Papovich, C. and Peth, M. and Ryan, R. and Somerville, R. and Weiner, B. and Wilson, G.}, month = oct, year = {2013}, pages = {93} } -@article{schmidt_evaluation_2020, - title = {Evaluation of probabilistic photometric redshift estimation approaches for {LSST}}, - volume = {2001.03621}, - journal = {arXiv}, - author = {Schmidt, S. J. and Malz, A. I. and Soo, J. Y. H. and Almosallam, I. A. and Brescia, M. and Cavuoti, S. and Cohen-Tanugi, J. and Connolly, A. J. and DeRose, J. and Freeman, P. E. and Graham, M. L. and Iyer, K. G. and Jarvis, M. J. and Kalmbach, J. B. and Kovacs, E. and Lee, A. B. and Longo, G. and Morrison, C. B. and Newman, J. A. and Nourbakhsh, E. and Nuss, E. and Pospisil, T. and Tranin, H. and Wechsler, R. H. and Zhou, R. and Izbicki, R. and Collaboration, The LSST Dark Energy Science}, - month = jan, - year = {2020}, -} - @article{hogg_inferring_2010, title = {{INFERRING} {THE} {ECCENTRICITY} {DISTRIBUTION}}, volume = {725}, issn = {0004-637X}, + url = {https://doi.org/10.1088%2F0004-637x%2F725%2F2%2F2166}, + doi = {10.1088/0004-637X/725/2/2166}, language = {en}, number = {2}, journal = {ApJ}, author = {Hogg, David W. and Myers, Adam D. and Bovy, Jo}, month = dec, year = {2010}, + note = {Publisher: IOP Publishing}, pages = {2166--2175}, } @@ -654,6 +748,8 @@ @article{gelman_inference_1992 title = {Inference from {Iterative} {Simulation} {Using} {Multiple} {Sequences}}, volume = {7}, issn = {0883-4237, 2168-8745}, + url = {https://projecteuclid.org/euclid.ss/1177011136}, + doi = {10.1214/ss/1177011136}, language = {EN}, number = {4}, journal = {Statist. Sci.}, @@ -661,5 +757,110 @@ @article{gelman_inference_1992 month = nov, year = {1992}, zmnumber = {06853057}, + note = {Publisher: Institute of Mathematical Statistics}, pages = {457--472}, } + +@article{padmanabhan_calibrating_2005, + title = {Calibrating photometric redshifts of luminous red galaxies}, + volume = {359}, + issn = {0035-8711}, + url = {https://academic.oup.com/mnras/article/359/1/237/983557}, + doi = {10.1111/j.1365-2966.2005.08915.x}, + language = {en}, + number = {1}, + journal = {Mon Not R Astron Soc}, + author = {Padmanabhan, Nikhil and Budav{\'a}ri, Tam{\'a}s and Schlegel, David J. and Bridges, Terry and Brinkmann, Jonathan and Cannon, Russell and Connolly, Andrew J. and Croom, Scott M. and Csabai, Istv{\'a}n and Drinkwater, Michael and Eisenstein, Daniel J. and Hewett, Paul C. and Loveday, Jon and Nichol, Robert C. and Pimbblet, Kevin A. and De Propris, Roberto and Schneider, Donald P. and Scranton, Ryan and Seljak, Uro{\v s} and Shanks, Tom and Szapudi, Istv{\'a}n and Szalay, Alexander S. and Wake, David}, + month = may, + year = {2005}, + note = {Publisher: Oxford Academic}, + pages = {237--250}, +} + +@article{awan_angular_2020, + title = {Angular {Correlation} {Function} {Estimators} {Accounting} for {Contamination} from {Probabilistic} {Distance} {Measurements}}, + volume = {890}, + issn = {0004-637X}, + url = {https://doi.org/10.3847%2F1538-4357%2Fab63c8}, + doi = {10.3847/1538-4357/ab63c8}, + language = {en}, + number = {1}, + journal = {ApJ}, + author = {Awan, Humna and Gawiser, Eric}, + month = feb, + year = {2020}, + note = {Publisher: American Astronomical Society}, + pages = {78}, +} + +@article{schmidt_evaluation_2020, + title = {Evaluation of probabilistic photometric redshift estimation approaches for {The} {Rubin} {Observatory} {Legacy} {Survey} of {Space} and {Time} ({LSST})}, + volume = {499}, + issn = {0035-8711}, + url = {https://academic.oup.com/mnras/article/499/2/1587/5905416}, + doi = {10.1093/mnras/staa2799}, + language = {en}, + number = {2}, + journal = {Mon Not R Astron Soc}, + author = {Schmidt, S. J. and Malz, A. I. and Soo, J. Y. H. and Almosallam, I. A. and Brescia, M. and Cavuoti, S. and Cohen-Tanugi, J. and Connolly, A. J. and DeRose, J. and Freeman, P. E. and Graham, M. L. and Iyer, K. G. and Jarvis, M. J. and Kalmbach, J. B. and Kovacs, E. and Lee, A. B. and Longo, G. and Morrison, C. B. and Newman, J. A. and Nourbakhsh, E. and Nuss, E. and Pospisil, T. and Tranin, H. and Wechsler, R. H. and Zhou, R. and Izbicki, R. and LSST~Dark~Energy~Science~Collaboration, The}, + month = oct, + year = {2020}, + note = {Publisher: Oxford Academic}, + pages = {1587--1606}, +} + +@article{malz_how_2021, + title = {How not to obtain the redshift distribution from probabilistic redshift estimates: {Under} what conditions is it not inappropriate to estimate the redshift distribution {N}(z) by stacking photo-z {PDFs}?}, + volume = {103}, + shorttitle = {How not to obtain the redshift distribution from probabilistic redshift estimates}, + url = {http://arxiv.org/abs/2101.04675}, + doi = {10.1103/PhysRevD.103.083502}, + number = {8}, + journal = {PRD}, + author = {Malz, Alex I.}, + month = apr, + year = {2021}, +} + +@article{storey-fisher_two-point_2020, + title = {Two-point statistics without bins: {A} continuous-function generalization of the correlation function estimator for large-scale structure}, + shorttitle = {Two-point statistics without bins}, + url = {http://arxiv.org/abs/2011.01836}, + journal = {arXiv:2011.01836 [astro-ph]}, + author = {Storey-Fisher, Kate and Hogg, David W.}, + month = nov, + year = {2020}, + note = {arXiv: 2011.01836}, +} + +@article{the_lsst_dark_energy_science_collaboration_lsst_2018, + title = {The {LSST} {Dark} {Energy} {Science} {Collaboration} ({DESC}) {Science} {Requirements} {Document}}, + url = {http://arxiv.org/abs/1809.01669}, + journal = {arXiv:1809.01669 [astro-ph]}, + author = {The LSST Dark Energy Science Collaboration and Mandelbaum, Rachel and Eifler, Tim and Hlo{\v z}ek, Ren{\'e}e and Collett, Thomas and Gawiser, Eric and Scolnic, Daniel and Alonso, David and Awan, Humna and Biswas, Rahul and Blazek, Jonathan and Burchat, Patricia and Chisari, Nora Elisa and Dell'Antonio, Ian and Digel, Seth and Frieman, Josh and Goldstein, Daniel A. and Hook, Isobel and Ivezi{\'c}, {\v Z}eljko and Kahn, Steven M. and Kamath, Sowmya and Kirkby, David and Kitching, Thomas and Krause, Elisabeth and Leget, Pierre-Fran{\c c}ois and Marshall, Philip J. and Meyers, Joshua and Miyatake, Hironao and Newman, Jeffrey A. and Nichol, Robert and Rykoff, Eli and Sanchez, F. Javier and Slosar, An{\v z}e and Sullivan, Mark and Troxel, M. A.}, + month = sep, + year = {2018}, +} + +@article{chang_effective_2013, + title = {The effective number density of galaxies for weak lensing measurements in the {LSST} project}, + volume = {434}, + issn = {0035-8711}, + url = {https://doi.org/10.1093/mnras/stt1156}, + doi = {10.1093/mnras/stt1156}, + number = {3}, + journal = {Monthly Notices of the Royal Astronomical Society}, + author = {Chang, C. and Jarvis, M. and Jain, B. and Kahn, S. M. and Kirkby, D. and Connolly, A. and Krughoff, S. and Peng, E.-H. and Peterson, J. R.}, + month = sep, + year = {2013}, + pages = {2121--2135}, +} + +@book{brooks_handbook_2011, + address = {Boca Raton, FL}, + title = {Handbook of {Markov} {Chain} {Monte} {Carlo}}, + isbn = {978-1-4200-7941-8}, + publisher = {Chapman and Hall/CRC}, + author = {Brooks, Steve and Gelman, Andrew and Jones, Galin L. and Meng, Xiao-Li}, + year = {2011}, +} diff --git a/research/paper/figures/chippr/pgm.png b/research/paper/pgm.png similarity index 100% rename from research/paper/figures/chippr/pgm.png rename to research/paper/pgm.png diff --git a/research/paper/figures/chippr/thesis_neghivarbias_log_estimators.png b/research/paper/results_bias.png similarity index 100% rename from research/paper/figures/chippr/thesis_neghivarbias_log_estimators.png rename to research/paper/results_bias.png diff --git a/research/paper/figures/chippr/thesis_eout_log_estimators.png b/research/paper/results_outlier_template.png similarity index 100% rename from research/paper/figures/chippr/thesis_eout_log_estimators.png rename to research/paper/results_outlier_template.png diff --git a/research/paper/figures/chippr/thesis_rout_log_estimators.png b/research/paper/results_outlier_training.png similarity index 100% rename from research/paper/figures/chippr/thesis_rout_log_estimators.png rename to research/paper/results_outlier_training.png diff --git a/research/paper/figures/chippr/single_uout_log_estimators.png b/research/paper/results_outlier_uniform.png similarity index 100% rename from research/paper/figures/chippr/single_uout_log_estimators.png rename to research/paper/results_outlier_uniform.png diff --git a/research/paper/figures/chippr/single_lsst_tmpr_log_estimators.png b/research/paper/results_prior_template.png similarity index 100% rename from research/paper/figures/chippr/single_lsst_tmpr_log_estimators.png rename to research/paper/results_prior_template.png diff --git a/research/paper/figures/chippr/single_lsst_trpr_log_estimators.png b/research/paper/results_prior_training.png similarity index 100% rename from research/paper/figures/chippr/single_lsst_trpr_log_estimators.png rename to research/paper/results_prior_training.png diff --git a/research/paper/figures/chippr/thesis_hivarsig_log_estimators.png b/research/paper/results_scatter_hi.png similarity index 100% rename from research/paper/figures/chippr/thesis_hivarsig_log_estimators.png rename to research/paper/results_scatter_hi.png diff --git a/research/paper/figures/chippr/single_varsigmas_log_estimators.png b/research/paper/results_scatter_lo.png similarity index 100% rename from research/paper/figures/chippr/single_varsigmas_log_estimators.png rename to research/paper/results_scatter_lo.png diff --git a/research/paper/figures/chippr/single_lsst_tmpr_wrong_log_estimators.png b/research/paper/results_wrongprior_template.png similarity index 100% rename from research/paper/figures/chippr/single_lsst_tmpr_wrong_log_estimators.png rename to research/paper/results_wrongprior_template.png diff --git a/research/paper/figures/chippr/single_lsst_trpr_wrong_log_estimators.png b/research/paper/results_wrongprior_training.png similarity index 100% rename from research/paper/figures/chippr/single_lsst_trpr_wrong_log_estimators.png rename to research/paper/results_wrongprior_training.png diff --git a/research/scripts/basic_replot_script.py b/research/scripts/basic_replot_script.py index 6b18665..f9d7c93 100644 --- a/research/scripts/basic_replot_script.py +++ b/research/scripts/basic_replot_script.py @@ -101,6 +101,7 @@ def just_plot(given_key): result_dir = os.path.join('..', 'results') test_name = 'thesis_neghibias' + watermark = None all_tests = {} test_info = {} diff --git a/research/scripts/cosmolike_sandbox.ipynb b/research/scripts/cosmolike_sandbox.ipynb index 7b659f1..c48eecb 100644 --- a/research/scripts/cosmolike_sandbox.ipynb +++ b/research/scripts/cosmolike_sandbox.ipynb @@ -16,6 +16,7 @@ "mpl.use('PS')\n", "import matplotlib.pyplot as plt\n", "from matplotlib.patches import Ellipse\n", + "from mpl_toolkits.axes_grid1 import make_axes_locatable\n", "\n", "mpl.rcParams['mathtext.rm'] = 'serif'\n", "plt.rcParams[\"mathtext.fontset\"] = \"dejavuserif\"\n", @@ -43,10 +44,22 @@ "\n", "import chippr\n", "import chippr.plot_utils as pu\n", + "import chippr.defaults as d\n", "\n", "%matplotlib inline" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import pylab\n", + "from scipy.stats import kde\n", + "import multiprocessing as mp" + ] + }, { "cell_type": "code", "execution_count": null, @@ -106,7 +119,10 @@ "outputs": [], "source": [ "chippr_bins = 35\n", - "colors = ['k', pu.colors[2], pu.colors[-3], pu.colors[0]]#'#332288', '#117733', '#999933', '#FF6700']" + "# bin_colors = [pu.colors[1], pu.colors[2], pu.colors[5], pu.colors[0]]\n", + "bin_colors = [pu.colors[5], pu.colors[2], pu.colors[1], pu.colors[6]]\n", + "colors = bin_colors\n", + "# colors = ['k', pu.colors[2], pu.colors[-3], pu.colors[0]]#'#332288', '#117733', '#999933', '#FF6700']" ] }, { @@ -117,9 +133,7 @@ "source": [ "cl_input = np.genfromtxt('../results/CosmoLike/nz_histo.txt')\n", "n_tomobins = np.shape(cl_input)[1] - 1\n", - "fine_dif = np.mean(cl_input.T[0])\n", - "cl_data = np.empty((chippr_bins, 5))\n", - "factor = len(cl_input.T[0]) / chippr_bins" + "fine_dif = np.mean(cl_input.T[0])" ] }, { @@ -128,10 +142,12 @@ "metadata": {}, "outputs": [], "source": [ - "for i in range(chippr_bins):\n", - " cl_data[i][0] = cl_input[i * factor][0]\n", - " cl_data[i][1:] = np.sum(cl_input[i * factor:(i + 1) * factor, 1:], axis=0) / (factor)\n", - "np.savetxt('../results/CosmoLike/nz_format_test.txt', cl_data)" + "# cl_data = np.empty((chippr_bins, 5))\n", + "# factor = len(cl_input.T[0]) / chippr_bins\n", + "# for i in range(chippr_bins):\n", + "# cl_data[i][0] = cl_input[i * factor][0]\n", + "# cl_data[i][1:] = np.sum(cl_input[i * factor:(i + 1) * factor, 1:], axis=0) / (factor)\n", + "# np.savetxt('../results/CosmoLike/nz_format_test.txt', cl_data)\n" ] }, { @@ -140,15 +156,24 @@ "metadata": {}, "outputs": [], "source": [ - "# this is the truth!\n", - "cl_data = np.empty((chippr_bins, 5))\n", - "for i in range(chippr_bins):\n", - " cl_data[i][1:] = np.sum(cl_input[i * 10:(i + 1) * 10, 1:], axis=0)\n", - " cl_data[i][0] = cl_input[i * 10][0]\n", - "for i in range(n_tomobins):\n", - " plt.plot(cl_data.T[0], cl_data.T[i+1])\n", + "# # this is the truth!\n", + "# cl_data = np.empty((chippr_bins, 5))\n", + "# for i in range(chippr_bins):\n", + "# cl_data[i][1:] = np.sum(cl_input[i * 10:(i + 1) * 10, 1:], axis=0)\n", + "# cl_data[i][0] = cl_input[i * 10][0]\n", + "# # for i in range(n_tomobins):\n", + "# # plt.plot(cl_data.T[0], cl_data.T[i+1])\n", " \n", - "# np.savetxt('../results/CosmoLike/nz_format_test.txt', cl_data)" + "# # np.savetxt('../results/CosmoLike/nz_format_test.txt', cl_data)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cl_data = np.loadtxt('../results/CosmoLike/nz_format_test.txt')" ] }, { @@ -184,7 +209,15 @@ "metadata": {}, "outputs": [], "source": [ - "print(bin_ends)" + "# all_true_nzs = [bin_mids]\n", + "# for i in range(n_tomobins):\n", + "# true_zs = np.genfromtxt('../results/CosmoLike/' + str(i) + 'single_lsst-in-thesis/data/true_vals.txt').T[0]\n", + "# true_nz = np.histogram(true_zs, bins=bin_ends)[0] / float(len(true_zs)) / bin_difs\n", + "# all_true_nzs.append(true_nz)\n", + "# all_true_nzs = np.array(all_true_nzs)\n", + "# np.savetxt('../results/CosmoLike/true_nz.txt', all_true_nzs.T)\n", + "\n", + "all_true_nzs = np.loadtxt('../results/CosmoLike/true_nz.txt').T" ] }, { @@ -193,13 +226,15 @@ "metadata": {}, "outputs": [], "source": [ - "all_true_nzs = [bin_mids]\n", - "for i in range(n_tomobins):\n", - " true_zs = np.genfromtxt('../results/CosmoLike/' + str(i) + 'single_lsst-in-thesis/data/true_vals.txt').T[0]\n", - " true_nz = np.histogram(true_zs, bins=bin_ends)[0] / float(len(true_zs)) / bin_difs\n", - " all_true_nzs.append(true_nz)\n", - "all_true_nzs = np.array(all_true_nzs)\n", - "np.savetxt('../results/CosmoLike/true_nz.txt', all_true_nzs.T)" + "# for i in range(n_tomobins):\n", + "# plt.plot(all_true_nzs[0], all_true_nzs[i+1])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## plotting the drawn values" ] }, { @@ -208,8 +243,51 @@ "metadata": {}, "outputs": [], "source": [ + "f = plt.figure(figsize=(7.5, 5))\n", + "ax = f.add_subplot(1, 1, 1)\n", + "ax.set_title(r'')\n", + "ax.set_xlabel(r'$z$', fontsize=18)\n", + "ax.set_ylabel(r'$n(z)$', fontsize=18)\n", + "\n", + "# ax.plot([-1., -2], [0., 1.], color='k', lw=1., alpha=1., label=r'underlying tomographic $n(z)$')\n", + "# ax.plot([-1., -2], [0., 1.], color='k', lw=2., alpha=0.5, linestyle='-', dashes=[3, 2], label=r'binned tomographic $n(z)$')\n", + "# ax.plot([-1., -2], [0., 1.], color='k', lw=2., alpha=1., linestyle='-', dashes=[1, 1], label=r'drawn tomographic $n(z)$')\n", + "# histy.plot([-1.], [-1.], color='k', lw=1.5, alpha=0.8, label=r'$n(z_{true})$')\n", + "ax.plot([-1.], [-1.], color='k', lw=5., alpha=0.4, label=r'$n(z_{est})$')\n", + "ax.plot([-1., -2], [0., 1.], color='k', lw=0.75, alpha=1., label=r'underlying tomographic $n(z)$')\n", + "ax.plot([-1., -2], [0., 1.], color='k', lw=1., alpha=0.6, linestyle='-', dashes=[3, 2], label=r'binned tomographic $n(z)$')\n", + "ax.plot([-1., -2], [0., 1.], color='k', lw=1.5, alpha=0.8, label=r'drawn tomographic $n(z_{true})$')\n", + "\n", "for i in range(n_tomobins):\n", - " plt.plot(all_true_nzs[0], all_true_nzs[i+1])" + " zs = np.genfromtxt('../results/CosmoLike/'+str(i)+'single_lsst-in-thesis/data/true_vals.txt')\n", + " zs = zs.T\n", + " true_zs = zs[0]\n", + " obs_zs = zs[1]\n", + "# ax.plot(cl_input.T[0], cl_input.T[i+1], color=colors[i], lw=1., alpha=1.)\n", + "# norm_factor = np.sum(bin_difs * cl_data.T[i+1])\n", + "# # print(norm_factor)\n", + "# pu.plot_step(ax, cl_data.T[0], cl_data.T[i+1] / norm_factor, c=colors[i], w=2., s='--', a=0.5, d=[(0,(3, 2))])\n", + "# pu.plot_step(ax, bin_ends, all_true_nzs[i+1], c=colors[i], w=2., s='--', a=0.5, d=[(0,(1, 1))])\n", + " ax.hist(true_zs, bins=bin_ends, histtype='step', linewidth=1.5, alpha=0.8, color=colors[i], stacked=False, density=True)\n", + "# histx.plot([-1., -2], [0., 1.], color='k', lw=1., alpha=1., label=r'underlying tomographic $n(z)$')\n", + " ax.plot(cl_input.T[0], cl_input.T[i+1], color=colors[i], lw=0.75, alpha=1.)\n", + " norm_factor = np.sum(bin_difs * cl_data.T[i+1])\n", + "# histx.plot([-1., -2], [0., 1.], color='k', lw=1., alpha=0.5, linestyle='-', dashes=[3, 2], label=r'binned tomographic $n(z)$')\n", + " pu.plot_step(ax, cl_data.T[0], cl_data.T[i+1] / norm_factor, c=colors[i], w=1., s='-', a=0.6, d=[(0,(3, 2))])\n", + "# histx.plot([-1., -2], [0., 1.], color='k', lw=1., alpha=1., linestyle='--', dashes=[1, 1], label=r'drawn tomographic $n(z)$')\n", + " pu.plot_step(ax, bin_ends, all_true_nzs[i+1], c=colors[i], w=1.5, a=0.8)#s='--', a=1., d=[(0,(1, 1))])\n", + " ax.hist(obs_zs, bins=bin_ends, alpha=0.4, color=colors[i], stacked=False, density=True)\n", + "ax.set_xlim(0., 3.5)\n", + "ax.set_ylim(0., 3.5)\n", + "ax.legend(loc='upper right', fontsize=16)\n", + "f.savefig('../results/CosmoLike/cosmolike_inputs.png', bbox_inches='tight', pad_inches = 0, dpi=250)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "warning: super slow!" ] }, { @@ -218,7 +296,18 @@ "metadata": {}, "outputs": [], "source": [ - "cl_input.T[0]" + "# def help_kde(i):\n", + "# zs = np.genfromtxt('../results/CosmoLike/'+str(i)+'single_lsst-in-thesis/data/true_vals.txt').T\n", + "# k = kde.gaussian_kde(zs)\n", + "# with open('../results/CosmoLike/'+str(i)+'single_lsst-in-thesis/data/scatter_kdes_small.pkl', 'wb') as dumpfile:\n", + "# cpkl.dump(k, dumpfile)\n", + "# return None\n", + "\n", + "# nps = mp.cpu_count()\n", + "# pool = mp.Pool(nps - 4)\n", + "# pool.map(help_kde, range(n_tomobins))\n", + "# # with open('../results/CosmoLike/scatter_kdes.pkl', 'wb') as dumpfile:\n", + "# # cpkl.dump(kdelist, dumpfile)" ] }, { @@ -227,26 +316,147 @@ "metadata": {}, "outputs": [], "source": [ - "f = plt.figure(figsize=(10, 5))\n", - "ax = f.add_subplot(1, 1, 1)\n", - "ax.set_title(r'')\n", - "ax.set_xlabel(r'$z$', fontsize=18)\n", - "ax.set_ylabel(r'$n(z)$', fontsize=18)\n", + "# def help_eval(i):\n", + "# z_grid = cl_input.T[0]\n", + "# nbins = len(z_grid)\n", + "# mingrid, maxgrid = z_grid[0], z_grid[-1]\n", + "# xi, yi = np.mgrid[mingrid:maxgrid:nbins*1j, mingrid:maxgrid:nbins*1j]\n", + "# # k = kde.gaussian_kde(zs)\n", + "# # k = kdelist[i]\n", + "# with open('../results/CosmoLike/'+str(i)+'single_lsst-in-thesis/data/scatter_kdes_small.pkl', 'rb') as dumpfile:\n", + "# k = cpkl.load(dumpfile)\n", + "# zi = k(np.vstack([xi.flatten(), yi.flatten()]))\n", "\n", - "ax.plot([-1., -2], [0., 1.], color='k', lw=1., alpha=1., label=r'underlying tomographic $n(z)$')\n", - "ax.plot([-1., -2], [0., 1.], color='k', lw=2., alpha=0.5, linestyle='-', dashes=[3, 2], label=r'binned tomographic $n(z)$')\n", - "ax.plot([-1., -2], [0., 1.], color='k', lw=2., alpha=1., linestyle='-', dashes=[1, 1], label=r'drawn tomographic $n(z)$')\n", + "# with open('../results/CosmoLike/'+str(i)+'single_lsst-in-thesis/data/scatter_contours.pkl', 'wb') as dumpfile:\n", + "# cpkl.dump(zi, dumpfile)\n", + "# return None\n", "\n", - "for i in range(n_tomobins):\n", - " ax.plot(cl_input.T[0], cl_input.T[i+1], color=colors[i], lw=1., alpha=1.)\n", - " norm_factor = np.sum(bin_difs * cl_data.T[i+1])\n", - " print(norm_factor)\n", - " pu.plot_step(ax, cl_data.T[0], cl_data.T[i+1] / norm_factor, c=colors[i], w=2., s='--', a=0.5, d=[(0,(3, 2))])\n", - " pu.plot_step(ax, bin_ends, all_true_nzs[i+1], c=colors[i], w=2., s='--', a=0.5, d=[(0,(1, 1))])\n", + "# nps = mp.cpu_count()\n", + "# pool = mp.Pool(nps-4)\n", + "# pool.map(help_eval, range(n_tomobins))\n", + "# # with open('../results/CosmoLike/scatter_kdes.pkl', 'wb') as dumpfile:\n", + "# # cpkl.dump(kdelist, dumpfile)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# bincolors = ['Purples', 'Blues', 'Greens', 'Oranges']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# with open('../results/CosmoLike/'+str(0)+'single_lsst-in-thesis/data/scatter_contours.pkl', 'rb') as dumpfile:\n", + "# # # k = cpkl.load(dumpfile)\n", + "# # # zi = k(np.vstack([xi.flatten(), yi.flatten()]))\n", + "# zi = np.log(cpkl.load(dumpfile))\n", + "# plt.hist(zi)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def plot_all_scatter(plot_loc='', prepend='', plot_name='combo_mega_scatter.png'):\n", + " z_grid = cl_input.T[0]\n", + " nbins = len(z_grid)\n", + " grid_ends, grid_mids = bin_ends, bin_mids\n", + " pu.set_up_plot()\n", + " f, scatplot = plt.subplots(1, 1, figsize=(7.5, 7.5))\n", + " divider = make_axes_locatable(scatplot)\n", + " histx = divider.append_axes('top', 1.2, pad=0., sharex=scatplot)\n", + " histy = divider.append_axes('right', 1.2, pad=0., sharey=scatplot)\n", " \n", - "ax.set_xlim(-0.1, 3.6)\n", - "ax.legend(loc='upper right', fontsize=20)\n", - "f.savefig('cosmolike_inputs.png', bbox_inches='tight', pad_inches = 0, dpi=250)" + " scatplot.plot(z_grid, z_grid, color='k', alpha=0.5, linewidth=1.5)\n", + " histy.plot([-1.], [-1.], color='k', lw=1.25, alpha=0.8, label=r'$n(z_{spec})$')\n", + " histy.plot([-1.], [-1.], color='k', lw=5., alpha=0.4, label=r'$n(z_{phot})$')\n", + " histx.plot([-1., -2], [0., 1.], color='k', lw=0.75, alpha=1., \n", + " label=r'underlying $n(z_{spec})$')\n", + " histx.plot([-1., -2], [0., 1.], color='k', lw=1., alpha=0.6, linestyle='-', dashes=[3, 2], \n", + " label=r'binned $n(z_{spec})$')\n", + " histx.plot([-1., -2], [0., 1.], color='k', lw=1.25, alpha=0.8, label=r'sampled $n(z_{spec})$')\n", + "\n", + " for i in range(n_tomobins):\n", + " scatplot.scatter([-1], [-1], color=colors[i], label='bin '+str(i))\n", + "\n", + " zs = np.genfromtxt(prepend+str(i)+'single_lsst-in-thesis/data/true_vals.txt')\n", + "\n", + " n = len(zs)\n", + " zs = zs.T\n", + " true_zs = zs[0]\n", + " obs_zs = zs[1]\n", + "\n", + " scatplot.scatter(true_zs, obs_zs, marker='.', s=2, alpha=0.15, color=colors[i])\n", + "\n", + " histy.hist(true_zs, bins=grid_ends, histtype='step', \n", + " linewidth=1.25, alpha=0.8, color=colors[i], orientation='horizontal', \n", + " stacked=False, density=True)\n", + " histy.hist(obs_zs, bins=grid_ends, orientation='horizontal', \n", + " alpha=0.4, color=colors[i], stacked=False, density=True)\n", + "\n", + "# histx.plot([-1., -2], [0., 1.], color='k', lw=1., alpha=1., label=r'underlying tomographic $n(z)$')\n", + " histx.plot(cl_input.T[0], cl_input.T[i+1], color=colors[i], lw=0.75, alpha=1.)\n", + " \n", + " norm_factor = np.sum(bin_difs * cl_data.T[i+1])\n", + " # print(norm_factor)\n", + " \n", + "# histx.plot([-1., -2], [0., 1.], color='k', lw=1., alpha=0.5, linestyle='-', dashes=[3, 2], label=r'binned tomographic $n(z)$')\n", + " pu.plot_step(histx, cl_data.T[0], cl_data.T[i+1] / norm_factor, \n", + " c=colors[i], w=1., s='-', a=0.6, d=[(0,(3, 2))])\n", + "# histx.plot([-1., -2], [0., 1.], color='k', lw=1., alpha=1., linestyle='--', dashes=[1, 1], label=r'drawn tomographic $n(z)$')\n", + " pu.plot_step(histx, bin_ends, all_true_nzs[i+1], \n", + " c=colors[i], w=1.25, a=0.8)#s='--', a=1., d=[(0,(1, 1))])\n", + " scatplot.legend(loc='lower right')\n", + " scatplot.set_xlim([min(grid_ends), max(grid_ends)])\n", + " scatplot.set_ylim([min(grid_ends), max(grid_ends)])\n", + " scatplot.set_xticks(np.linspace(np.floor(min(grid_ends)), np.floor(max(grid_ends)), 4))\n", + " scatplot.set_yticks(np.linspace(np.floor(min(grid_ends)), np.floor(max(grid_ends)), 4))\n", + " scatplot.set_xlabel(r'$z_{spec}$', fontsize=20)\n", + " scatplot.set_ylabel(r'$z_{phot}$', fontsize=20)\n", + " scatplot.text(0.25, 3., r'LSST-like mock $p(z\\mid \\mathrm{``data\"})$', rotation=0, size=20)\n", + "\n", + " # scatplot.set_aspect(1.)\n", + " histx.xaxis.set_tick_params(labelbottom=False)\n", + " histy.yaxis.set_tick_params(labelleft=False)\n", + " histx.set_ylabel(r'$n(z_{spec})$')\n", + " histy.set_xlabel(r'$n(z_{phot})$')\n", + " histx.set_yticks([])\n", + " histy.set_xticks([])\n", + " histx.set_xlim([min(grid_ends), max(grid_ends)])\n", + " histy.set_ylim([min(grid_ends), max(grid_ends)])\n", + " histx.set_ylim([0., 3.5])\n", + " histy.set_xlim([0., 3.5])\n", + " histx.legend(loc='upper right')\n", + " histy.legend(loc='upper right')\n", + " \n", + " f.subplots_adjust(hspace=0.01, wspace=0.01)\n", + " pylab.savefig(os.path.join(plot_loc, plot_name), bbox_inches='tight', pad_inches=0, dpi=d.dpi)\n", + " return" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_all_scatter(plot_loc='../results/CosmoLike/', prepend='../results/CosmoLike/')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## now run chippr separately" ] }, { @@ -266,19 +476,9 @@ "metadata": {}, "outputs": [], "source": [ - "all_stk_nzs, all_mle_nzs, all_map_nzs = [bin_mids], [bin_mids], [bin_mids]\n", - "for i in range(n_tomobins):\n", - " with open('../results/' + str(i) + 'single_lsst-in-thesis/results/nz.p', 'rb') as test_file:\n", - " test_info = cpkl.load(test_file)\n", - " all_stk_nzs.append(np.exp(test_info['estimators']['log_stacked_nz']))\n", - " all_mle_nzs.append(np.exp(test_info['estimators']['log_mmle_nz']))\n", - " all_map_nzs.append(np.exp(test_info['estimators']['log_mmap_nz']))\n", - "all_stk_nzs = np.array(all_stk_nzs)\n", - "all_mle_nzs = np.array(all_mle_nzs)\n", - "all_map_nzs = np.array(all_map_nzs)\n", - "np.savetxt('../results/CosmoLike/stack_nz.txt', all_stk_nzs.T)\n", - "np.savetxt('../results/CosmoLike/mmle_nz.txt', all_mle_nzs.T)\n", - "np.savetxt('../results/CosmoLike/mmap_nz.txt', all_map_nzs.T)" + "# with open('../results/' + str(0) + 'single_lsst/results/nz.p', 'rb') as test_file:\n", + "# test_info = cpkl.load(test_file)\n", + "# print(test_info)" ] }, { @@ -287,19 +487,86 @@ "metadata": {}, "outputs": [], "source": [ - "for i in range(n_tomobins):\n", - " plt.plot(all_true_nzs[0], all_true_nzs[i+1], c='k')\n", - " plt.plot(all_true_nzs[0], all_stk_nzs[i+1], c=pu.colors[0])\n", - " plt.plot(all_true_nzs[0], all_mle_nzs[i+1], c='b')\n", - " plt.plot(all_true_nzs[0], all_map_nzs[i+1], c=pu.colors[-3])\n", - " " + "# all_stk_nzs, all_mle_nzs, all_map_nzs = [bin_mids], [bin_mids], [bin_mids]\n", + "# for i in range(n_tomobins):\n", + "# with open('../results/' + str(i) + 'single_lsst/results/nz.p', 'rb') as test_file:\n", + "# test_info = cpkl.load(test_file)\n", + "# all_stk_nzs.append(np.exp(test_info['estimators']['log_stacked_nz']))\n", + "# all_mle_nzs.append(np.exp(test_info['estimators']['log_mmle_nz']))\n", + "# all_map_nzs.append(np.exp(test_info['estimators']['log_mmap_nz']))\n", + "# all_stk_nzs = np.array(all_stk_nzs)\n", + "# all_mle_nzs = np.array(all_mle_nzs)\n", + "# all_map_nzs = np.array(all_map_nzs)\n", + "# np.savetxt('../results/CosmoLike/stack_nz.txt', all_stk_nzs.T)\n", + "# np.savetxt('../results/CosmoLike/mmle_nz.txt', all_mle_nzs.T)\n", + "# np.savetxt('../results/CosmoLike/mmap_nz.txt', all_map_nzs.T)\n", + "\n", + "all_stk_nzs = np.loadtxt('../results/CosmoLike/stack_nz.txt').T\n", + "all_mle_nzs = np.loadtxt('../results/CosmoLike/mmle_nz.txt').T\n", + "all_map_nzs = np.loadtxt('../results/CosmoLike/mmap_nz.txt').T" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "# what to do with CosmoLike output" + "testnames = ['true_nz', 'stack_nz', 'mmle_nz', 'mmap_nz']\n", + "meth_colors = ['k', pu.colors[0], pu.colors[3], pu.colors[4]]\n", + "colors = meth_colors\n", + "styles = [[(0, (1, 0.001))], [(0, (2, 2))], [(0, (1, 2))], [(0, (2, 1))]]\n", + "names = ['true', 'stacked', 'CHIPPR', 'modes']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# f = plt.figure(figsize=(10, 5))\n", + "f, axs = plt.subplots(4, 1, figsize=(10, 20), sharex=True, sharey=True)\n", + "# ax = f.add_subplot(1, 1, 1)\n", + "# ax.set_title(r'')\n", + "\n", + "# ax.set_xlabel(r'$z$', fontsize=18)\n", + "ax.set_ylabel(r'$n(z)$', fontsize=18)\n", + "\n", + "# ax.plot([-1., -2], [0., 1.], color=colors[0], lw=1., alpha=1., label=r'true sample $n(z)$')\n", + "# ax.plot([-1., -2], [0., 1.], color=colors[1], lw=1.5, alpha=1., linestyle='--', dashes=[2, 2], label=r'stacked $n(z)$')\n", + "# ax.plot([-1., -2], [0., 1.], color=colors[2], lw=2., alpha=1., linestyle='-', dashes=[2, 1], label=r'CHIPPR $n(z)$')\n", + "# ax.plot([-1., -2], [0., 1.], color=colors[3], lw=1.5, alpha=1., linestyle='--', dashes=[2, 1], label=r'modes $n(z)$')\n", + "\n", + "# ax.plot([-1., -2], [0., 1.], color='k', lw=1., alpha=1., label=r'underlying tomographic $n(z)$')\n", + "# ax.plot([-1., -2], [0., 1.], color='k', lw=2., alpha=0.5, linestyle='-', dashes=[3, 2], label=r'binned tomographic $n(z)$')\n", + "# ax.plot([-1., -2], [0., 1.], color='k', lw=2., alpha=1., linestyle='-', dashes=[1, 1], label=r'drawn tomographic $n(z)$')\n", + "\n", + "for i in range(n_tomobins):\n", + " ax = axs[i]\n", + "# ax.plot(cl_input.T[0], cl_input.T[i+1], color=colors[i], lw=0.75, alpha=1.)\n", + "# norm_factor = np.sum(bin_difs * cl_data.T[i+1])\n", + "# print(norm_factor)\n", + "# pu.plot_step(ax, cl_data.T[0], cl_data.T[i+1] / norm_factor, c=colors[i], w=1., s='-', a=0.6, d=[(0,(3, 2))])\n", + " pu.plot_step(ax, bin_ends, all_true_nzs[i+1], c=colors[0], w=1., a=1., l=r'true sample $n(z)$')\n", + " pu.plot_step(ax, bin_ends, all_stk_nzs[i+1], c=colors[1], w=1.5, s='--', a=1., d=[(0,(2, 2))], l=r'stacked $n(z)$')\n", + " pu.plot_step(ax, bin_ends, all_mle_nzs[i+1], c=colors[2], w=2., s='-', a=1, d=[(0, (2, 1))], l=r'CHIPPR $n(z)$')\n", + " pu.plot_step(ax, bin_ends, all_map_nzs[i+1], c=colors[3], w=1.4, s='--', a=1., d=[(0,(2, 1))], l=r'modes $n(z)$')\n", + " \n", + "# ax.set_xlim(-0.1, 3.6)\n", + " ax.semilogy()\n", + " ax.set_ylabel(r'$\\ln[n(z)]$', fontsize=18)\n", + " ax.legend(loc='lower left', fontsize=14)\n", + " if i == 3:\n", + " ax.set_xlabel(r'$z$', fontsize=18)\n", + " ax.set_xlim(0., 3.5)\n", + " ax.set_xticks([0, 1, 2, 3])\n", + " ax.set_ylim(1e-4, 1e1)\n", + " ax.set_yticks([1e-3, 1e-2, 1e-1, 1e0])\n", + " ax.set_yticklabels([-3, -2, -1, 0])\n", + " ax.text(3., 1., r'bin '+str(i), rotation=0, size=20, backgroundcolor=bin_colors[i])\n", + "f.subplots_adjust(hspace=0., wspace=0.)\n", + "f.savefig('../results/CosmoLike/binned_nz_recovered.png', bbox_inches='tight', pad_inches = 0, dpi=250)" ] }, { @@ -308,7 +575,42 @@ "metadata": {}, "outputs": [], "source": [ - "testnames = ['true_nz', 'stack_nz', 'mmle_nz', 'mmap_nz']" + "# f = plt.figure(figsize=(10, 5))\n", + "# ax = f.add_subplot(1, 1, 1)\n", + "# ax.set_title(r'')\n", + "# ax.set_xlabel(r'$z$', fontsize=18)\n", + "# ax.set_ylabel(r'$n(z)$', fontsize=18)\n", + "\n", + "# ax.plot([-1., -2], [0., 1.], color=colors[0], lw=1., alpha=1., label=r'true sample $n(z)$')\n", + "# ax.plot([-1., -2], [0., 1.], color=colors[1], lw=1.5, alpha=1., linestyle='--', dashes=[2, 2], label=r'stacked $n(z)$')\n", + "# ax.plot([-1., -2], [0., 1.], color=colors[2], lw=2., alpha=1., linestyle='-', dashes=[2, 1], label=r'CHIPPR $n(z)$')\n", + "# ax.plot([-1., -2], [0., 1.], color=colors[3], lw=1.5, alpha=1., linestyle='--', dashes=[2, 1], label=r'modes $n(z)$')\n", + "\n", + "# # ax.plot([-1., -2], [0., 1.], color='k', lw=1., alpha=1., label=r'underlying tomographic $n(z)$')\n", + "# # ax.plot([-1., -2], [0., 1.], color='k', lw=2., alpha=0.5, linestyle='-', dashes=[3, 2], label=r'binned tomographic $n(z)$')\n", + "# # ax.plot([-1., -2], [0., 1.], color='k', lw=2., alpha=1., linestyle='-', dashes=[1, 1], label=r'drawn tomographic $n(z)$')\n", + "\n", + "# for i in range(n_tomobins):\n", + "# # ax.plot(cl_input.T[0], cl_input.T[i+1], color=colors[i], lw=0.75, alpha=1.)\n", + "# # norm_factor = np.sum(bin_difs * cl_data.T[i+1])\n", + "# # print(norm_factor)\n", + "# # pu.plot_step(ax, cl_data.T[0], cl_data.T[i+1] / norm_factor, c=colors[i], w=1., s='-', a=0.6, d=[(0,(3, 2))])\n", + "# pu.plot_step(ax, bin_ends, all_true_nzs[i+1], c=colors[0], w=1., a=1.)\n", + "# pu.plot_step(ax, bin_ends, all_stk_nzs[i+1], c=colors[1], w=1.5, s='--', a=1., d=[(0,(2, 2))])\n", + "# pu.plot_step(ax, bin_ends, all_mle_nzs[i+1], c=colors[2], w=2., s='-', a=1, d=[(0, (2, 1))])\n", + "# pu.plot_step(ax, bin_ends, all_map_nzs[i+1], c=colors[3], w=1.5, s='--', a=1., d=[(0,(2, 1))])\n", + " \n", + "# ax.semilogy()\n", + "# ax.set_xlim(-0.1, 3.6)\n", + "# ax.legend(loc='upper right', fontsize=20)\n", + "# f.savefig('../results/CosmoLike/binned_nz_recovered.png', bbox_inches='tight', pad_inches = 0, dpi=250)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# what to do with CosmoLike output" ] }, { @@ -531,9 +833,6 @@ "metadata": {}, "outputs": [], "source": [ - "colors = ['k', pu.colors[0], pu.colors[2], pu.colors[4]] \n", - "styles = [[(0, (1, 0.001))], [(0, (2, 2))], [(0, (1, 2))], [(0, (2, 1))]]\n", - "names = ['true', 'stacked', 'CHIPPR', 'modes']\n", "hatches = [None, '--', None, '||']\n", "fills = [False, False, True, False]\n", "alphas = [1., 1., 0.5, 1.]\n", @@ -557,6 +856,7 @@ " labels = keys\n", " \n", " ncol = len(keys)\n", + " plt.ticklabel_format(style='sci')\n", " fig = plt.figure(figsize=(ncol*(5), ncol*5))\n", " ax = [[fig.add_subplot(ncol, ncol, ncol * i + j + 1) for j in range(i+1)] for i in range(ncol)]\n", " to_keep = range(ncol)#[0, 1, 2, 4]\n", @@ -580,6 +880,8 @@ " lim_val_i = np.sqrt(all_invfisher[k][i][i])\n", " for j in range(i + 1):#to_keep[:k]:#range(i+1):\n", " lim_val_j = np.sqrt(all_invfisher[k][j][j])\n", + " ax[i][j].set_yticks([])\n", + " ax[i][j].set_yticks([])\n", "# ax[i][j].set_xlim(-1. * extrema[i], extrema[i])\n", " if i == j:\n", " lim_val = np.sqrt(all_invfisher[k][i][i])\n", @@ -588,20 +890,26 @@ " ax[i][j].plot(x_grid, func.pdf(x_grid), \n", " color=colors[k], label=names[k], alpha=1., linestyle=styles[k][0], linewidth=3.)\n", " ax[i][j].set_xlabel(labels[i], fontsize=24)\n", - " ax[i][j].set_yticks([])\n", " else:\n", " ellipse = Ellipse(xy=(0., 0.), height=2.*a[i][j], width=2.*b[i][j], angle=(90. - t[i][j]*180./np.pi), \n", " alpha=alphas[k], color=colors[k], linestyle=styles[k][0], lw=2., hatch=hatches[k], fill=fills[k])\n", " ax[i][j].add_artist(ellipse)\n", - " if j == 0:\n", - " ax[i][j].set_ylabel(labels[i], fontsize=24)\n", - " if i == ncol-1:\n", - " ax[i][j].set_xlabel(labels[j], fontsize=24)\n", + " \n", " # note x and y swapping here. . .\n", " extreme_y = a[i][j]*abs(np.cos(t[i][j]))+b[i][j]*abs(np.sin(t[i][j]))\n", " extreme_x = a[i][j]*abs(np.sin(t[i][j]))+b[i][j]*abs(np.cos(t[i][j]))\n", " ax[i][j].set_xlim(-1.5*extreme_y, 1.5*extreme_y)\n", " ax[i][j].set_ylim(-1.5*extreme_x, 1.5*extreme_x)\n", + " if j == 0:\n", + " ax[i][j].set_ylabel(labels[i], fontsize=24)\n", + " xsca = round(extreme_x, 4)\n", + " ax[i][j].set_yticks([-1.*extreme_x, 0., extreme_x])\n", + " ax[i][j].set_yticklabels([-1.*xsca, 0, xsca])\n", + " if i == ncol-1:\n", + " ax[i][j].set_xlabel(labels[j], fontsize=24)\n", + " ysca = round(extreme_y, 4)\n", + " ax[i][j].set_xticks([-1.*extreme_y, 0., extreme_y])\n", + " ax[i][j].set_xticklabels([-1.*ysca, 0, ysca])\n", " plt.subplots_adjust(hspace=0., wspace=0.)\n", " plt.savefig('../results/CosmoLike/final_plot'+extra+'.png', dpi=100, bbox_inches='tight', pad_inches = 0)\n", " return" @@ -619,13 +927,6 @@ "mycorner(fish_params, keys, names, hatches, fills, alphas, fancy)" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, { "cell_type": "markdown", "metadata": {}, diff --git a/research/scripts/multi_inf_script.py b/research/scripts/multi_inf_script.py index b7b718b..7993cb0 100644 --- a/research/scripts/multi_inf_script.py +++ b/research/scripts/multi_inf_script.py @@ -160,6 +160,7 @@ def do_inference(given_key): import pickle import os import multiprocessing as mp + nps = mp.cpu_count() import chippr from chippr import * @@ -174,6 +175,5 @@ def do_inference(given_key): test_info['name'] = test_name all_tests[test_name] = test_info - nps = mp.cpu_count() - pool = mp.Pool(nps) + pool = mp.Pool(nps-1) pool.map(do_inference, all_tests.keys()) diff --git a/research/scripts/multi_replot_script.py b/research/scripts/multi_replot_script.py index 6eec47a..54854b6 100644 --- a/research/scripts/multi_replot_script.py +++ b/research/scripts/multi_replot_script.py @@ -104,6 +104,10 @@ def just_plot(given_key): result_dir = os.path.join('..', 'results') name_file = 'which_plt_tests.txt' + watermarks = { + + } + with open(name_file) as tests_to_run: all_tests = {} for test_name in tests_to_run: diff --git a/research/scripts/multi_sim_script.py b/research/scripts/multi_sim_script.py index 1d25c85..75fce61 100644 --- a/research/scripts/multi_sim_script.py +++ b/research/scripts/multi_sim_script.py @@ -65,26 +65,26 @@ def make_true(given_key): if test_info['params']['smooth_truth'] == 1: # true_shape = test_info['params']['true_shape'] true_loc = test_info['params']['true_loc'] - true_scale = test_info['params']['true_scale'] + true_scale = test_info['params1']['true_scale'] true_nz = gamma(true_loc, true_scale**2, bounds=(min(test_info['bin_ends']), max(test_info['bin_ends'])))#sps.erlang(true_shape, true_loc, true_scale) true_dict = {'amps': [1.], 'means': [true_loc], 'sigmas': [true_scale]} else: bin_range = max(test_info['bin_ends']) - min(test_info['bin_ends']) - true_amps = np.array([0.20, 0.35, 0.55]) + true_amps = np.array([0.35, 0.55, 0.20]) true_means = np.array([0.5, 0.2, 0.75]) * bin_range + min(test_info['bin_ends']) - true_sigmas = np.array([0.4, 0.2, 0.1]) * bin_range + true_sigmas = np.array([0.2, 0.1, 0.4]) * bin_range n_mix_comps = len(true_amps) true_funcs = [] for c in range(n_mix_comps): true_funcs.append(chippr.gauss(true_means[c], true_sigmas[c]**2)) - true_nz = chippr.gmix(true_amps, true_funcs, limits=(min(test_info['bin_ends']), max(test_info['bin_ends']))) + true_nz = chippr.gmix(true_amps, true_funcs, limits=(min(test_info['bin_ends']), max(test_info['bin_ends']))) - true_dict = {'amps': true_amps, 'means': true_means, 'sigmas': true_sigmas} + true_dict = {'amps': true_amps, 'means': true_means, 'sigmas': true_sigmas} true_dict['bins'] = test_info['bin_ends'] # true_zs = true_nz.sample(test_info['params']['n_galaxies']) - # true_dict['zs'] = true_zs + # true_dict['zs'] = true_zsnps = mp.cpu_count() test_info['truth'] = true_dict return(test_info, true_nz) @@ -231,6 +231,7 @@ def make_catalog(given_key): import os import shutil import multiprocessing as mp + nps = mp.cpu_count() import scipy.interpolate as spi import chippr @@ -246,6 +247,8 @@ def make_catalog(given_key): test_info['name'] = test_name[:-1] all_tests[test_name] = test_info - nps = mp.cpu_count() - pool = mp.Pool(nps) - pool.map(make_catalog, all_tests.keys()) + if len(all_tests.keys()) > 1: + pool = mp.Pool(nps-1) + pool.map(make_catalog, all_tests.keys()) + else: + make_catalog(all_tests.keys()[0]) diff --git a/research/scripts/pedagogical-final.ipynb b/research/scripts/pedagogical-final.ipynb new file mode 100644 index 0000000..95cc5cb --- /dev/null +++ b/research/scripts/pedagogical-final.ipynb @@ -0,0 +1,1057 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import cPickle as cpkl\n", + "import csv\n", + "import numpy as np\n", + "import os\n", + "import pandas as pd\n", + "import pickle\n", + "from scipy.stats import kde\n", + "import scipy.stats as sps\n", + "from sklearn.neighbors import KernelDensity\n", + "import sys\n", + "\n", + "import chippr\n", + "import chippr.plot_utils as pu" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib as mpl\n", + "mpl.use('PS')\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.patches import Ellipse, Polygon\n", + "import matplotlib.cm as cm\n", + "from mpl_toolkits.axes_grid1 import make_axes_locatable\n", + "import matplotlib.hatch\n", + "from matplotlib.colors import ListedColormap, LinearSegmentedColormap\n", + "\n", + "\n", + "# from cosmolike_sandbox\n", + "mpl.rcParams['mathtext.rm'] = 'serif'\n", + "plt.rcParams[\"mathtext.fontset\"] = \"dejavuserif\"\n", + "mpl.rcParams['font.family'] = 'serif'\n", + "mpl.rcParams['font.serif'] = 'Times New Roman'\n", + "# mpl.rcParams['text.usetex'] = False\n", + "\n", + "# from pedagogical\n", + "mpl.rcParams['text.usetex'] = True\n", + "# mpl.rcParams['mathtext.rm'] = 'serif'\n", + "# mpl.rcParams['font.family'] = 'serif'\n", + "# mpl.rcParams['font.serif'] = 'DejaVu Serif'\n", + "\n", + "title = 18\n", + "label = 16\n", + "mpl.rcParams['axes.titlesize'] = title\n", + "mpl.rcParams['axes.labelsize'] = label\n", + "mpl.pyplot.rcParams['xtick.labelsize'] = label\n", + "mpl.pyplot.rcParams['ytick.labelsize'] = label\n", + "mpl.rcParams['figure.subplot.left'] = 0.2\n", + "mpl.rcParams['figure.subplot.right'] = 0.9\n", + "mpl.rcParams['figure.subplot.bottom'] = 0.2\n", + "mpl.rcParams['figure.subplot.top'] = 0.9\n", + "mpl.rcParams['figure.subplot.wspace'] = 0.5\n", + "mpl.rcParams['figure.subplot.hspace'] = 0.5\n", + "mpl.rcParams['savefig.dpi'] = 250\n", + "mpl.rcParams['savefig.format'] = 'pdf'\n", + "mpl.rcParams['savefig.bbox'] = 'tight'\n", + "\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# pedagogical plot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "points = pd.read_csv('jain05.csv', delimiter=',', names=['redshift', '\"data\"'])\n", + "\n", + "# remove outliers\n", + "to_use = points[np.logical_and(points['redshift'] > 0., points['\"data\"'] > 0)].values.T\n", + "\n", + "X, Y = np.linspace(0., 4.25, 100), np.linspace(0., 4.25, 100)\n", + "X, Y = np.meshgrid(X, Y)\n", + "xy = np.vstack([Y.ravel(), X.ravel()]).T\n", + "kde = KernelDensity(kernel='gaussian', bandwidth=0.1).fit(to_use.T)\n", + "\n", + "# slow!\n", + "Z = np.exp(kde.score_samples(xy))\n", + "\n", + "Z = Z.reshape(X.shape)\n", + "\n", + "for_posterior = points[np.logical_and(points['\"data\"'] > 2.9, points['\"data\"'] < 3.1)]['redshift'].values[:, np.newaxis]\n", + "for_likelihood = points[np.logical_and(points['redshift'] > 1.9, points['redshift'] < 2.1)]['\"data\"'].values[:, np.newaxis]\n", + "\n", + "plotgrid = np.linspace(0., 4.25, 100)\n", + "\n", + "kde_p = KernelDensity(kernel='gaussian', bandwidth=0.05).fit(for_posterior)\n", + "posterior = np.exp(kde_p.score_samples(plotgrid[:, np.newaxis]))\n", + "kde_l = KernelDensity(kernel='gaussian', bandwidth=0.05).fit(for_likelihood)\n", + "likelihood = np.exp(kde_l.score_samples(plotgrid[:, np.newaxis]))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "discrete_colors = np.array([(86, 180, 233), (0, 114, 178), (213, 94, 0), (230, 159, 0)])/256.#, axis=0)\n", + "# ['Sky blue', 'Blue', 'Vermilion', 'Orange']\n", + "# discrete_colors = np.array([(86, 180, 233), (0, 158, 115), (240, 228, 66)])/256.#, axis=0)\n", + "col_hor = np.array((213, 94, 0))/256.\n", + "col_ver = np.array((0, 114, 178))/256.\n", + "for i in range(len(discrete_colors)):\n", + " plt.scatter(i, i, color=discrete_colors[i])\n", + "\n", + " \n", + "mycmap = matplotlib.colors.LinearSegmentedColormap.from_list('mycmap', discrete_colors)#['#1B3B52', '#45C0CE',#'#258599', \n", + "# # '#5E8253', \n", + "# '#F7DB3F', \n", + "# '#F26811', \n", + "# # '#8F3311', \n", + "# # '#B43809', \n", + "# '#821308'])\n", + "# mycmap = ListedColormap([discrete_colors[-2], discrete_colors[1]])\n", + "plt.scatter(np.linspace(0., 6., 100), np.linspace(0., 6, 100), color=mycmap(np.linspace(0., 6., 100)/6.))\n", + "plt.hlines(3., 0., 6., color=col_hor)\n", + "plt.vlines(3., 0., 6., color=col_ver)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 2d KDE\n", + "f, scatplot = plt.subplots(figsize=(8.5, 7.5))\n", + "f.subplots_adjust(hspace=0)\n", + "plt.scatter(points['redshift'], points['\"data\"'], marker='.', s=2, color='k')\n", + "Z = Z.reshape(X.shape)\n", + "scatplot.contourf(Y, X, Z, cmap=mycmap, alpha=0.75, linestyles=None)# scatplot.colorbar()\n", + "plt.plot([0, 4.25], [0, 4.25], color=discrete_colors[-1], alpha=1.)\n", + "scatplot.vlines(3., 0., 4.25, linewidth=10., alpha=0.5, color=col_hor)\n", + "scatplot.hlines(2., 0., 4.25, linewidth=10., alpha=0.5, color=col_ver)\n", + "scatplot.set_xlabel(r'true redshift $z$', fontsize=20)\n", + "scatplot.set_ylabel('nonlinear projection of photometric data', fontsize=20)\n", + "scatplot.set_xlim(0., 4.25)\n", + "scatplot.set_ylim(0., 4.25)\n", + "scatplot.set_yticks([])#0., 1., 2., 3., 4])\n", + "scatplot.set_xticks([0., 1., 2., 3., 4])\n", + "divider = make_axes_locatable(scatplot)\n", + "histx = divider.append_axes('top', 1.2, pad=0., sharex=scatplot)\n", + "histy = divider.append_axes('right', 1.2, pad=0., sharey=scatplot)\n", + "histx.xaxis.set_tick_params(labelbottom=False)\n", + "histy.yaxis.set_tick_params(labelleft=False)\n", + "histx.plot(plotgrid, likelihood, color=col_ver)\n", + "histy.plot(posterior, plotgrid, color=col_hor)\n", + "histx.set_yticks([])\n", + "histx.set_xlim(0., 4.25)\n", + "histy.set_xticks([])\n", + "histy.set_ylim(0., 4.25)\n", + "histx.text(0.5, 1.25, r'posterior\\\\$p(z\\mid \\mathrm{``data\"})$', rotation=0, size=16)\n", + "histy.text(0.85, 2., r'likelihood\\\\$p(\\mathrm{``data\"}\\mid z)$', rotation=-90, size=16)\n", + "f.savefig('jain15.png', bbox_inches='tight', pad_inches=0, dpi=250)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# scratch after here, actual cosmolike work/plots in other notebook" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Examining CosmoLike input" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# # print(np.shape(cl_input))\n", + "# # n_tomobins = np.shape(cl_input)[1] - 1\n", + "# # 1 column of redshifts, 4 tomographic bins, 350 redshift bins\n", + "# # too many redshift bins, I'll set the truth to that binned down" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# chippr_bins = 35\n", + "# colors = ['k', pu.colors[2], pu.colors[-3], pu.colors[0]]#'#332288', '#117733', '#999933', '#FF6700']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# cl_input = np.genfromtxt('../results/CosmoLike/nz_histo.txt')\n", + "# n_tomobins = np.shape(cl_input)[1] - 1\n", + "# fine_dif = np.mean(cl_input.T[0])\n", + "# cl_data = np.empty((chippr_bins, 5))\n", + "# factor = len(cl_input.T[0]) / chippr_bins" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# for i in range(chippr_bins):\n", + "# cl_data[i][0] = cl_input[i * factor][0]\n", + "# cl_data[i][1:] = np.sum(cl_input[i * factor:(i + 1) * factor, 1:], axis=0) / (factor)\n", + "# np.savetxt('../results/CosmoLike/nz_format_test.txt', cl_data)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# # this is the truth!\n", + "# cl_data = np.empty((chippr_bins, 5))\n", + "# for i in range(chippr_bins):\n", + "# cl_data[i][1:] = np.sum(cl_input[i * 10:(i + 1) * 10, 1:], axis=0)\n", + "# cl_data[i][0] = cl_input[i * 10][0]\n", + "# for i in range(n_tomobins):\n", + "# plt.plot(cl_data.T[0], cl_data.T[i+1])\n", + " \n", + "# # np.savetxt('../results/CosmoLike/nz_format_test.txt', cl_data)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# make CHIPPR output into CosmoLike format" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## truth\n", + "\n", + "in units of counts" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# bin_ends = np.genfromtxt('../results/CosmoLike/0single_lsst-in-thesis/data/metadata.txt')\n", + "# bin_difs = bin_ends[1:] - bin_ends[:-1]\n", + "# bin_mids = (bin_ends[1:] + bin_ends[:-1]) / 2." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# all_true_nzs = [bin_mids]\n", + "# for i in range(n_tomobins):\n", + "# true_zs = np.genfromtxt('../results/CosmoLike/' + str(i) + 'single_lsst-in-thesis/data/true_vals.txt').T[0]\n", + "# true_nz = np.histogram(true_zs, bins=bin_ends)[0] / float(len(true_zs)) / bin_difs\n", + "# all_true_nzs.append(true_nz)\n", + "# all_true_nzs = np.array(all_true_nzs)\n", + "# np.savetxt('../results/CosmoLike/true_nz.txt', all_true_nzs.T)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# for i in range(n_tomobins):\n", + "# plt.plot(all_true_nzs[0], all_true_nzs[i+1])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# cl_input.T[0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# f = plt.figure(figsize=(10, 5))\n", + "# ax = f.add_subplot(1, 1, 1)\n", + "# ax.set_xlabel(r'$z$', fontsize=18)\n", + "# ax.set_ylabel(r'$n(z)$', fontsize=18)\n", + "\n", + "# ax.plot([-1., -2], [0., 1.], color='k', lw=1., alpha=1., label=r'underlying tomographic $n(z)$')\n", + "# ax.plot([-1., -2], [0., 1.], color='k', lw=2., alpha=0.5, linestyle='-', dashes=[3, 2], label=r'binned tomographic $n(z)$')\n", + "# ax.plot([-1., -2], [0., 1.], color='k', lw=2., alpha=1., linestyle='-', dashes=[1, 1], label=r'drawn tomographic $n(z)$')\n", + "\n", + "# for i in range(n_tomobins):\n", + "# ax.plot(cl_input.T[0], cl_input.T[i+1], color=colors[i], lw=1., alpha=1.)\n", + "# norm_factor = np.sum(bin_difs * cl_data.T[i+1])\n", + "# print(norm_factor)\n", + "# pu.plot_step(ax, cl_data.T[0], cl_data.T[i+1] / norm_factor, c=colors[i], w=2., s='--', a=0.5, d=[(0,(3, 2))])\n", + "# pu.plot_step(ax, bin_ends, all_true_nzs[i+1], c=colors[i], w=2., s='--', a=0.5, d=[(0,(1, 1))])\n", + " \n", + "# ax.set_xlim(-0.1, 3.6)\n", + "# ax.legend(loc='upper right', fontsize=20)\n", + "# f.savefig('cosmolike_inputs.png', bbox_inches='tight', pad_inches = 0, dpi=250)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## estimators\n", + "\n", + "log n(z) piecewise constant and separate files per tomobin each with all formats\n", + "\n", + "need single file with z values, point evaluations in each tomobin" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# all_stk_nzs, all_mle_nzs, all_map_nzs = [bin_mids], [bin_mids], [bin_mids]\n", + "# for i in range(n_tomobins):\n", + "# with open('../results/' + str(i) + 'single_lsst-in-thesis/results/nz.p', 'rb') as test_file:\n", + "# test_info = cpkl.load(test_file)\n", + "# all_stk_nzs.append(np.exp(test_info['estimators']['log_stacked_nz']))\n", + "# all_mle_nzs.append(np.exp(test_info['estimators']['log_mmle_nz']))\n", + "# all_map_nzs.append(np.exp(test_info['estimators']['log_mmap_nz']))\n", + "# all_stk_nzs = np.array(all_stk_nzs)\n", + "# all_mle_nzs = np.array(all_mle_nzs)\n", + "# all_map_nzs = np.array(all_map_nzs)\n", + "# np.savetxt('../results/CosmoLike/stack_nz.txt', all_stk_nzs.T)\n", + "# np.savetxt('../results/CosmoLike/mmle_nz.txt', all_mle_nzs.T)\n", + "# np.savetxt('../results/CosmoLike/mmap_nz.txt', all_map_nzs.T)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# for i in range(n_tomobins):\n", + "# plt.plot(all_true_nzs[0], all_true_nzs[i+1], c='k')\n", + "# plt.plot(all_true_nzs[0], all_stk_nzs[i+1], c=pu.colors[0])\n", + "# plt.plot(all_true_nzs[0], all_mle_nzs[i+1], c='b')\n", + "# plt.plot(all_true_nzs[0], all_map_nzs[i+1], c=pu.colors[-3])\n", + " " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# what to do with CosmoLike output" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# testnames = ['true_nz', 'stack_nz', 'mmle_nz', 'mmap_nz']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# all_covariances = []\n", + "# all_invcovariances = []\n", + "# magfactors = []\n", + "# for testname in testnames:\n", + "# covariances = np.zeros((200, 200)) + sys.float_info.min\n", + "# # with open(os.path.join('../results/CosmoLike', 'Cl_cov.nz'+testname+'.txt'), 'rb') as cosmolike_file:\n", + "# # cosmolike_reader = csv.reader(cosmolike_file, delimiter=' ')\n", + "# # cosmolike_reader.next()\n", + "# # for row in cosmolike_reader:\n", + "# # # covariance matrix is filled with positive values << floating point precision.\n", + "# # # I'm going to inflate them and deflate before reporting results\n", + "# # covariances[int(row[0])][int(row[1])] = float(row[2])\n", + "# # covariances[int(row[1])][int(row[0])] = float(row[2])\n", + "# covariance_table = np.genfromtxt('../results/CosmoLike/Cl_cov.'+testname+'.txt')\n", + "# magfactor = np.asarray(np.min(covariance_table.T[-1]))\n", + "# print(magfactor)\n", + "# magfactors.append(magfactor)\n", + "# for row in covariance_table:\n", + "# # print(row)\n", + "# covariances[int(row[0])][int(row[1])] = row[2]# / magfactor\n", + "# covariances[int(row[1])][int(row[0])] = row[2]# / magfactor\n", + "# # covariances = covariances - sys.float_info.min\n", + "# # assert(np.all(covariances >= 0.))\n", + "# # covariances += sys.float_info.min\n", + "# # covariances = covariances / magfactor\n", + "# assert(np.all(covariances > 0.))\n", + "# all_covariances.append(covariances)\n", + "# invcovariances = np.linalg.inv(covariances)\n", + "# mask = np.where((invcovariances < sys.float_info.min))\n", + "# invcovariances[mask] = sys.float_info.min\n", + "# print(np.shape(invcovariances))\n", + "# assert(np.all(invcovariances > 0.))\n", + "# all_invcovariances.append(invcovariances)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# # # plt.imshow(np.log(covariances))\n", + "# # # print(covariances)\n", + "\n", + "# # invcovariances = np.linalg.pinv(covariances)# / sys.float_info.epsilon) * sys.float_info.epsilon\n", + "# # # print(invcovariances)\n", + "\n", + "# # # # invcovariances = np.linalg.pinv(1.e15 * covariances) * 1.e15\n", + "# # # # print(invcovariances)\n", + "\n", + "# # # plt.imshow(np.log(invcovariances * sys.float_info.epsilon))\n", + "# # # print(invcovariances * sys.float_info.epsilon)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# # Cl, dCldOm, dClds8, dCldns, dCldw0, dCldwa, dCldOb, dCldH0 = np.zeros(200), np.zeros(200), np.zeros(200), np.zeros(200), np.zeros(200), np.zeros(200), np.zeros(200)\n", + "# all_deriv_info = []\n", + "# for testname in testnames:\n", + "# # deriv_info = np.zeros((200, 7)) + sys.float_info.min\n", + "# # with open(os.path.join('../results/CosmoLike', 'Cl_derivs.nz'+testname+'.txt'), 'rb') as cosmolike_file:\n", + "# # cosmolike_reader = csv.reader(cosmolike_file, delimiter=' ')\n", + "# # cosmolike_reader.next()\n", + "# # i = 0\n", + "# # while i < 200:\n", + "# # for row in cosmolike_reader:\n", + "# # for j in range(7):\n", + "# # deriv_info[i][j] = float(row[j])\n", + "# # i += 1\n", + "# deriv_info = np.genfromtxt('../results/CosmoLike/Cl_derivs.'+testname+'.txt')#, skip_header=0)\n", + "# # print(np.shape(deriv_info))\n", + "# deriv_info = np.vstack((np.zeros(7), deriv_info))\n", + "# deriv_info = deriv_info.T\n", + "# all_deriv_info.append(deriv_info)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# # # deriv_info = deriv_info.T\n", + "# # dCldOm, dClds8, dCldns, dCldw0, dCldwa, dCldOb, dCldH0 = deriv_info[0], deriv_info[1], deriv_info[2], deriv_info[3], deriv_info[4], deriv_info[5], deriv_info[6]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# all_fisher = []\n", + "# all_invfisher = []\n", + "# for k in range(len(testnames)):\n", + "# ls = np.arange(200)\n", + "# fisher = np.eye(7)\n", + "# for i in range(7):\n", + "# for j in range(i+1):\n", + "# fisher[i][j] = np.sum((2. * ls[1:] + 1.) * all_deriv_info[k][i][1:] * \n", + "# (all_invcovariances[k][1:, 1:]) * all_deriv_info[k][j][1:])\n", + "# fisher[j][i] = fisher[i][j]\n", + "# all_fisher.append(fisher)\n", + "# invfisher = np.linalg.pinv(fisher)#np.linalg.inv(fisher/sys.float_info.epsilon) * sys.float_info.epsilon\n", + "# all_invfisher.append(invfisher)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# # # print(fisher/sys.float_info.epsilon)\n", + "# # # plt.imshow(fisher/sys.float_info.epsilon)\n", + "\n", + "# # invfisher = np.linalg.inv(fisher/sys.float_info.epsilon) * sys.float_info.epsilon\n", + "# # plt.imshow(np.log(all_fisher[0]))\n", + "# # plt.colorbar()\n", + "# # print(invfisher)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# def make_ellipse_params(fisher):\n", + "# diag_elems = np.diag(fisher)\n", + "# term1 = (diag_elems[:, np.newaxis] + diag_elems[np.newaxis, :]) / 2.\n", + "# term2 = np.sqrt((diag_elems[:, np.newaxis] - diag_elems[np.newaxis, :])**2 / 4. + fisher * fisher.T)\n", + "# a = np.sqrt(term1 + term2)\n", + "# b = np.sqrt(term1 - term2)\n", + "# t = np.arctan((fisher + fisher.T) / np.abs(diag_elems[:, np.newaxis] - diag_elems[np.newaxis, :])) / 2.\n", + "# try:\n", + "# assert(np.all(a >= 0.))\n", + "# except:\n", + "# print(a)\n", + "# try:\n", + "# assert(np.all(b >= 0.))\n", + "# except:\n", + "# print(b)\n", + "# return(a, b, t)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# fish_params = []\n", + "# for k in range(len(testnames)):\n", + "# fisher = all_invfisher[k]\n", + "# (a, b, t) = make_ellipse_params(fisher)\n", + "# fish_params.append((a, b, t))\n", + "# fish_params = np.array(fish_params)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# def kld(p, q):\n", + "# term1 = np.trace(np.matmul(np.linalg.pinv(q), p))\n", + "# # print(term1)\n", + "# term2 = float(len(keys))\n", + "# term3 = np.log(np.linalg.det(q) / np.linalg.det(p))\n", + "# # print(term3)\n", + "# return (term1 - term2 + term3) / 2." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# def area(a, b):\n", + "# return np.pi * a * b" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# for k in range(1, len(testnames)):\n", + "# print(testnames[k])\n", + "# the_kld = kld(all_invfisher[0], all_invfisher[k])\n", + "# print(the_kld)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# colors = ['k', pu.colors[0], pu.colors[2], pu.colors[4]] \n", + "# styles = [[(0, (1, 0.001))], [(0, (2, 2))], [(0, (1, 2))], [(0, (2, 1))]]\n", + "# names = ['true', 'stack', 'CHIPPR', 'MAP']\n", + "# hatches = [None, '--', '..', '||']\n", + "# keys = ['dCl/dOm', 'dCl/ds8', 'dCl/dns', 'dCl/dw0', 'dCl/dwa', 'dCl/dOb', 'dCl/dH0']\n", + "# fancy = [r'$dC_{l}/d\\Omega_{m}$', r'$dC_{l}/ds_{8}$', r'$dC_{l}/dn_{s}$', r'$dC_{l}/dw_{0}$', r'$dC_{l}/dw_{a}$', r'$dC_{l}/d\\Omega_{b}$',r'$dC_{l}/dH_{0}$']\n", + "# mini_inds = [0, 3, 4, -1]\n", + "# mini_keys = [keys[mini_ind] for mini_ind in mini_inds]\n", + "# mini_fancy = [fancy[mini_ind] for mini_ind in mini_inds]\n", + "# # keys = [r'$\\Omega_{m}$', r'$S_{8}$', r'$n_{s}$', r'$w_{0}$', r'$w_{a}$', r'$\\Omega_{b}$', r'$H_{0}$']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# def mycorner(fish_params, keys, names, hatches, labels=None, extra=''):\n", + " \n", + "# if labels==None:\n", + "# labels = keys\n", + " \n", + "# ncol = len(keys)\n", + "# fig = plt.figure(figsize=(ncol*(5), ncol*5))\n", + "# ax = [[fig.add_subplot(ncol, ncol, ncol * i + j + 1) for j in range(i+1)] for i in range(ncol)]\n", + "# to_keep = range(ncol)#[0, 1, 2, 4]\n", + " \n", + "# extrema = np.zeros(ncol)\n", + "# for k in range(len(testnames)):\n", + "# if k != 0:\n", + "# the_kld = '\\n KLD='+\"{0:.3e}\".format(kld(all_invfisher[0], all_invfisher[k]))\n", + "# else:\n", + "# the_kld = ''\n", + "# axl = fig.add_subplot(ncol, ncol, ncol+ncol-1)\n", + "# axl.text(0.1, 0.9 - 0.25 * k, names[k]+r' $n(z)$'+the_kld, fontsize=30, color=colors[k])\n", + "# axl.set_xticks([])\n", + "# axl.set_yticks([])\n", + "# # fisher = fishers_info[k]\n", + " \n", + "# a, b, t = fish_params[k]\n", + "# extrema = np.max(np.sqrt(a), axis=0) / 2.\n", + " \n", + "# for i in range(ncol):\n", + "# lim_val_i = np.sqrt(all_invfisher[k][i][i])\n", + "# for j in range(i + 1):#to_keep[:k]:#range(i+1):\n", + "# lim_val_j = np.sqrt(all_invfisher[k][j][j])\n", + "# # ax[i][j].set_xlim(-1. * extrema[i], extrema[i])\n", + "# if i == j:\n", + "# lim_val = np.sqrt(all_invfisher[k][i][i])\n", + "# x_grid = np.linspace(-3. * lim_val_i, 3. * lim_val_i, 100)\n", + "# func = sps.norm(0., lim_val)\n", + "# ax[i][j].plot(x_grid, func.pdf(x_grid), \n", + "# color=colors[k], label=names[k], alpha=0.75, linestyle=styles[k][0], linewidth=2.)\n", + "# ax[i][j].set_xlabel(labels[i], fontsize=24)\n", + "# ax[i][j].set_yticks([])\n", + "# else:\n", + "# ellipse = Ellipse(xy=(0., 0.), height=2.*a[i][j], width=2.*b[i][j], angle=(90. - t[i][j]*180./np.pi), \n", + "# alpha=1., color=colors[k], linestyle=styles[k][0], lw=2., hatch=hatches[k], fill=False)\n", + "# ax[i][j].add_artist(ellipse)\n", + "# if j == 0:\n", + "# ax[i][j].set_ylabel(labels[i], fontsize=24)\n", + "# if i == ncol-1:\n", + "# ax[i][j].set_xlabel(labels[j], fontsize=24)\n", + "# # note x and y swapping here. . .\n", + "# extreme_y = a[i][j]*abs(np.cos(t[i][j]))+b[i][j]*abs(np.sin(t[i][j]))\n", + "# extreme_x = a[i][j]*abs(np.sin(t[i][j]))+b[i][j]*abs(np.cos(t[i][j]))\n", + "# ax[i][j].set_xlim(-2.*extreme_y, 2.*extreme_y)\n", + "# ax[i][j].set_ylim(-2.*extreme_x, 2.*extreme_x)\n", + "# plt.subplots_adjust(hspace=0., wspace=0.)\n", + "# plt.savefig('../results/CosmoLike/final_plot'+extra+'.png', dpi=100, bbox_inches='tight', pad_inches = 0)\n", + "# return" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# # mycorner(fish_params, mini_keys, names, mini_fancy, 'mini')\n", + "# mycorner(fish_params, keys, names, hatches, fancy)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Calculate the KLD between ellipses" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# making CosmoLike input" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# test_dir = '../results/single'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# simulated_posteriors = chippr.catalog(params='single.txt', loc=test_dir)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# data = simulated_posteriors.read(loc='data', style='.txt')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# params = chippr.utils.ingest('single.txt')\n", + "# def check_prob_params(params):\n", + "# \"\"\"\n", + "# Sets parameter values pertaining to components of probability\n", + "\n", + "# Parameters\n", + "# ----------\n", + "# params: dict\n", + "# dictionary containing key/value pairs for probability\n", + "\n", + "# Returns\n", + "# -------\n", + "# params: dict\n", + "# dictionary containing key/value pairs for probability\n", + "# \"\"\"\n", + "# if 'prior_mean' not in params:\n", + "# params['prior_mean'] = 'interim'\n", + "# else:\n", + "# params['prior_mean'] = params['prior_mean'][0]\n", + "# if 'no_prior' not in params:\n", + "# params['no_prior'] = 0\n", + "# else:\n", + "# params['no_prior'] = int(params['no_prior'][0])\n", + "# if 'no_data' not in params:\n", + "# params['no_data'] = 0\n", + "# else:\n", + "# params['no_data'] = int(params['no_data'][0])\n", + "# return params\n", + "# params = check_prob_params(params)\n", + "# def set_up_prior(data, params):\n", + "# \"\"\"\n", + "# Function to create prior distribution from data\n", + "\n", + "# Parameters\n", + "# ----------\n", + "# data: dict\n", + "# catalog dictionary containing bin endpoints, log interim prior, and log\n", + "# interim posteriors\n", + "# params: dict\n", + "# dictionary of parameter values for creation of prior\n", + "\n", + "# Returns\n", + "# -------\n", + "# prior: chippr.mvn object\n", + "# prior distribution as multivariate normal\n", + "# \"\"\"\n", + "# zs = data['bin_ends']\n", + "# log_nz_intp = data['log_interim_prior']\n", + "# log_z_posts = data['log_interim_posteriors']\n", + "\n", + "# z_difs = zs[1:]-zs[:-1]\n", + "# z_mids = (zs[1:]+zs[:-1])/2.\n", + "# n_bins = len(z_mids)\n", + "\n", + "# n_pdfs = len(log_z_posts)\n", + "\n", + "# a = 1.# / n_bins\n", + "# b = 20.#1. / z_difs ** 2\n", + "# c = a / n_pdfs\n", + "# prior_var = np.eye(n_bins)\n", + "# for k in range(n_bins):\n", + "# prior_var[k] = a * np.exp(-0.5 * b * (z_mids[k] - z_mids) ** 2)\n", + "# prior_var += c * np.identity(n_bins)\n", + "\n", + "# prior_mean = log_nz_intp\n", + "# prior = chippr.mvn(prior_mean, prior_var)\n", + "# if params['prior_mean'] == 'sample':\n", + "# new_mean = prior.sample_one()\n", + "# prior = chippr.mvn(new_mean, prior_var)\n", + "# print(params['prior_mean'], prior_mean, new_mean)\n", + "# else:\n", + "# print(params['prior_mean'], prior_mean)\n", + "\n", + "# return (prior, prior_var)\n", + "# (prior, cov) = set_up_prior(data, params)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# with open(os.path.join('../results/single/data', 'true_params.p'), 'r') as true_file:\n", + "# true_nz_info = pickle.load(true_file)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# true_nz_info['amps']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# true_funcs = [chippr.discrete(np.array([true_nz_info['bins'][0], true_nz_info['bins'][-1]]), true_nz_info['amps'])]\n", + "# true_amps = true_nz_info['amps']\n", + "# # # print(true_amps)\n", + "# true_nz = true_funcs[0]#chippr.gmix(true_amps, true_funcs, limits=(min(true_nz_info['bins']), max(true_nz_info['bins'])))\n", + "# grid_mids = (true_nz_info['bins'][1:] + true_nz_info['bins'][:-1])/2.\n", + "# plt.plot(grid_mids, true_nz.evaluate(grid_mids))\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# nz = chippr.log_z_dens(data, prior, truth=true_nz, loc='../results/single', vb=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# nz.read('nz.p')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# bins_to_write = np.linspace(0.0101, 3.5001, 350)\n", + "# empty_bins = np.random.random((350, 4))/350.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# bin_mids = (nz.info['bin_ends'][1:] + nz.info['bin_ends'][:-1])/2.\n", + "# with open(os.path.join('../results/single/results', 'nz_mmle_test.txt'), 'wb') as cosmolike_file:\n", + "# # cosmolike_file.write(zip(bin_mids, np.exp(nz.info['estimators']['log_mmle_nz']))\n", + "# cosmolike_writer = csv.writer(cosmolike_file, delimiter=' ')\n", + "# cosmolike_writer.writerows(zip(bin_mids, np.exp(nz.info['estimators']['log_mmle_nz'])))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Placeholder" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Revisiting!" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "How to run cosmolike. . ." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# plot samples from prior" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# z_0 = 0.3\n", + "# def smooth_func(z):\n", + "# return 1/(2 * z_0) * (z/z_0)**2 * np.exp(-z/z_0)\n", + "# zs = np.linspace(0., 1., 100)\n", + "\n", + "# nz = smooth_func(zs[:-1])\n", + "# nz /= np.dot(nz, zs[1:]-zs[:-1])\n", + "# plt.plot(zs[:-1], nz)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# print(np.dot(true_funcs[0].evaluate(grid_mids), true_nz_info['grid'][1:]-true_nz_info['grid'][:-1]))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# print(np.dot(true_funcs[0].distweights, true_funcs[0].dbins))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# func = sps.norm(0.25, 0.05)\n", + "# print(func.std())\n", + "# x = np.linspace(0., 1., 100)\n", + "# plt.plot(x, func.pdf(x))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# n_bins = 100\n", + "# z_mids = np.linspace(0., 1., n_bins)\n", + "\n", + "# a = 1.\n", + "# b = 20.#mid-scale variations, larger means more peaks\n", + "# c = 1.e-6#longest-scale variation, lower increases amplitude relative to small-scale\n", + "# prior_var = np.eye(n_bins)\n", + "# for k in range(n_bins):\n", + "# prior_var[k] = a * np.exp(-0.5 * b * (z_mids[k] - z_mids) ** 2)\n", + "# prior_var += c * np.identity(n_bins)\n", + "\n", + "# prior_mean = np.zeros(n_bins)\n", + "# prior = chippr.mvn(prior_mean, prior_var)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# samples = prior.sample(7)\n", + "# for each in samples:\n", + "# plt.plot(z_mids, each)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "chippr (Python 2)", + "language": "python", + "name": "chippr_2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.15rc1" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/research/scripts/single_lsst.txt b/research/scripts/single_lsst.txt index ba7fdf4..14d6cd8 100644 --- a/research/scripts/single_lsst.txt +++ b/research/scripts/single_lsst.txt @@ -1,16 +1,14 @@ n_gals 4 -n_bins 15 +n_bins 35 bin_min 0.0101 bin_max 3.4101 smooth_truth 1 true_loc 3 true_scale 0.3 variable_sigmas 1 -constant_sigma 0.05 +constant_sigma 0.02 catastrophic_outliers uniform outlier_fraction 0.1 -outlier_mean 1.0 -outlier_sigma 0.01 ez_bias 1 variable_bias 1 ez_bias_val -0.003 diff --git a/research/scripts/single_lsst_tmpr_wrong.txt b/research/scripts/single_lsst_tmpr_wrong.txt new file mode 100644 index 0000000..1e76df9 --- /dev/null +++ b/research/scripts/single_lsst_tmpr_wrong.txt @@ -0,0 +1,18 @@ +n_gals 4 +n_bins 35 +bin_min 0.0101 +bin_max 3.4101 +smooth_truth 1 +true_loc 3 +true_scale 0.3 +variable_sigmas 1 +constant_sigma 0.02 +catastrophic_outliers uniform +outlier_fraction 0.1 +ez_bias 1 +variable_bias 1 +ez_bias_val -0.003 +interim_prior template +wrong_prior 1 +n_walkers 100 +gr_threshold 1.5 diff --git a/research/scripts/single_uout.txt b/research/scripts/single_uout.txt new file mode 100644 index 0000000..5b79e37 --- /dev/null +++ b/research/scripts/single_uout.txt @@ -0,0 +1,15 @@ +n_gals 4 +n_bins 35 +bin_min 0.0 +bin_max 3.5 +smooth_truth 1 +true_loc 3 +true_scale 0.3 +variable_sigmas 0 +constant_sigma 0.02 +catastrophic_outliers uniform +outlier_fraction 0.1 +ez_bias 0 +variable_bias 0 +n_walkers 100 +gr_threshold 1.5 diff --git a/research/scripts/single_varsigmas.txt b/research/scripts/single_varsigmas.txt new file mode 100644 index 0000000..c650a77 --- /dev/null +++ b/research/scripts/single_varsigmas.txt @@ -0,0 +1,13 @@ +n_gals 4 +n_bins 35 +bin_min 0.0 +bin_max 3.5 +smooth_truth 1 +true_loc 3 +true_scale 0.3 +variable_sigmas 1 +constant_sigma 0.02 +catastrophic_outliers 0 +ez_bias 0 +n_walkers 100 +gr_threshold 1.5 diff --git a/research/scripts/thesis_eout.txt b/research/scripts/thesis_eout.txt new file mode 100644 index 0000000..87358b4 --- /dev/null +++ b/research/scripts/thesis_eout.txt @@ -0,0 +1,18 @@ +n_gals 4 +n_bins 35 +bin_min 0.0 +bin_max 3.5 +smooth_truth 1 +true_loc 3 +true_scale 0.3 +variable_sigmas 0 +constant_sigma 0.02 +catastrophic_outliers template +outlier_fraction 0.1 +outlier_mean 1.75 +outlier_sigma 0.02 +ez_bias 0 +variable_bias 0 +ez_bias_val 0. +n_walkers 100 +gr_threshold 1.2 diff --git a/research/scripts/thesis_hivarsig.txt b/research/scripts/thesis_hivarsig.txt new file mode 100644 index 0000000..50f9fcf --- /dev/null +++ b/research/scripts/thesis_hivarsig.txt @@ -0,0 +1,13 @@ +n_gals 4 +n_bins 35 +bin_min 0.0 +bin_max 3.5 +smooth_truth 1 +true_loc 3 +true_scale 0.3 +variable_sigmas 1 +constant_sigma 0.04 +catastrophic_outliers 0 +ez_bias 0 +n_walkers 100 +gr_threshold 1.5 diff --git a/research/scripts/thesis_neghivarbias.txt b/research/scripts/thesis_neghivarbias.txt new file mode 100644 index 0000000..7da3a9d --- /dev/null +++ b/research/scripts/thesis_neghivarbias.txt @@ -0,0 +1,15 @@ +n_gals 4 +n_bins 35 +bin_min 0.0 +bin_max 3.5 +smooth_truth 1 +true_loc 3 +true_scale 0.3 +variable_sigmas 0 +constant_sigma 0.02 +catastrophic_outliers 0 +ez_bias 1 +variable_bias 1 +ez_bias_val -0.03 +n_walkers 100 +gr_threshold 1.5 diff --git a/research/scripts/thesis_negvarbias.txt b/research/scripts/thesis_negvarbias.txt new file mode 100644 index 0000000..2950ff7 --- /dev/null +++ b/research/scripts/thesis_negvarbias.txt @@ -0,0 +1,15 @@ +n_gals 4 +n_bins 35 +bin_min 0.0 +bin_max 3.5 +smooth_truth 1 +true_loc 3 +true_scale 0.3 +variable_sigmas 0 +constant_sigma 0.02 +catastrophic_outliers 0 +ez_bias 1 +variable_bias 1 +ez_bias_val -0.003 +n_walkers 100 +gr_threshold 1.5 diff --git a/research/scripts/thesis_rout.txt b/research/scripts/thesis_rout.txt new file mode 100644 index 0000000..628404c --- /dev/null +++ b/research/scripts/thesis_rout.txt @@ -0,0 +1,18 @@ +n_gals 4 +n_bins 35 +bin_min 0.0 +bin_max 3.5 +smooth_truth 1 +true_loc 3 +true_scale 0.3 +variable_sigmas 0 +constant_sigma 0.02 +catastrophic_outliers training +outlier_fraction 0.1 +outlier_mean 1.75 +outlier_sigma 0.02 +ez_bias 0 +variable_bias 0 +ez_bias_val 0. +n_walkers 100 +gr_threshold 1.2