<commit_before>\documentclass[iop]{emulateapj}

\usepackage{tikz}
\usepackage{natbib}
\usepackage{amsmath}

\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]

\newcommand{\myemail}{<EMAIL>}
\newcommand{\textul}{\underline}
\newcommand{\chippr}{\texttt{CHIPPR} }

\shorttitle{How to obtain the redshift distribution from probabilistic redshift 
estimates}
\shortauthors{Malz, et al.}

\begin{document}

\title{How to obtain the redshift distribution from probabilistic redshift 
estimates}

\author{<NAME>\altaffilmark{1}, <NAME>\altaffilmark{1,2,3,4}, Phil 
Marshall\altaffilmark{5}, \& Others}
\email{<EMAIL>}

\altaffiltext{1}{Center for Cosmology and Particle Physics, Department of 
Physics,
  New York University, 726 Broadway, 9th floor, New York, NY 10003, USA}
\altaffiltext{2}{Simons Center for Computational Astrophysics, 162 Fifth 
Avenue, 7th floor, New York, NY 10010, USA}
\altaffiltext{3}{Center for Data Science, New York University, 60 Fifth Avenue, 
7th floor, New York, NY 10003, USA}
\altaffiltext{4}{Max-Planck-Institut f\"ur Astronomie, K\"onigstuhl 17, D-69117 
Heidelberg, Germany}
\altaffiltext{5}{[SLAC]}

\begin{abstract}
The redshift distribution $n(z)$ is a crucial ingredient for weak lensing 
cosmology.  Spectroscopically confirmed redshifts for the dim and numerous 
galaxies observed by weak lensing surveys are expected to be inaccessible, 
making photometric redshifts (photo-$z$s) the next best alternative.  Because 
of the nontrivial inference involved in their determination, photo-$z$ point 
estimates are being superseded by photo-$z$ probability distribution functions 
(PDFs).  However, analytic methods for utilizing these new data products in 
cosmological inference are still evolving.  This paper presents a novel 
approach to estimating the posterior distribution over $n(z)$ from a survey of 
galaxy photo-$z$ PDFs based upon a probabilistic graphical model of 
hierarchical inference.  We present the Cosmological Hierarchical Inference 
with Probabilistic Photometric Redshifts (\chippr) code implementing this 
technique, as well as its validation on mock data and testing on the 
\texttt{Buzzard} simulations.  \chippr yields a more accurate characterization 
of $n(z)$ containing information beyond the best-fit estimator produced by 
traditional procedures.  The publicly available code is easily extensible to 
other one-point statistics that depend on redshift.

\end{abstract}

\keywords{catalogs --- cosmology: cosmological parameters --- galaxies: 
statistics --- gravitational lensing: weak --- methods: analytical --- methods: 
data analysis --- methods: statistical --- techniques: photometric}

\section{Introduction}
\label{sec:introduction}

After a brief literature review addressing how photo-$z$ PDFs are currently 
used in cosmology, this paper aims to answer the following questions:

\begin{itemize}
	\item Why should we question existing methods?
	\item How can we improve the effectiveness of using photo-$z$ PDFs in 
inference?
	\item How does the result of \chippr compare to established estimators 
in terms of the accuracy of $n(z)$?
	\item How significant is the effect of the discrepancy between $n(z)$ 
estimators on cosmological constraints?
\end{itemize}

\section{Method}
\label{sec:method}

This paper presents a mathematically consistent method for obtaining the 
posterior probability distribution over the redshift distribution $n(z)$ using 
a catalog of photo-$z$ PDFs.  We start by introducing some nomenclature, 
definitions, and symbols that will be used throughout Secs. \ref{sec:model}, 
\ref{sec:others}, and \ref{sec:implementation}.

We shall say the parameters comprising $\vec{\phi}$ define the redshift 
distribution under some functional form about which this method is agnostic, so 
long as the functional form evaluated with some true values of the parameters 
accurately describes the true $n(z)$.  Since the redshift distribution is 
itself a probability distribution, it may be written as $p(z | \vec{\phi})$.

Before jumping into the details of the model, it is crucial to settle on an 
interpretation of what a photo-$z$ PDF actually is.  When a photo-$z$ PDF is 
reported, it is an object containing information about the redshift $z_{i}$ of 
galaxy $i$ based on its photometric data $\vec{d}_{i}$, making a photo-$z$ PDF 
a \textit{posterior} probability distribution.  The data $\vec{d}_{i}$ may be 
fluxes, magnitudes, colors, or any combination thereof.  We assume that the 
redshifts of the galaxies are not related to the redshifts or photometry of 
other galaxies in the survey, i.e. they are \textit{independent}.

In addition to the data $\vec{d}_{i}$, the photo-$z$ PDF also contains 
information imparted by assumptions that went into the process by which the 
photo-$z$ was made.  Those assumptions are only relevant insomuch as they 
impose a preference on the redshift distribution. Thus, they can be reduced to 
a sort of prior on the distribution of redshifts, which can be parametrized 
under the chosen functional form by parameters $\vec{\phi}^{*}$.  The redshift 
distribution $p(z | \vec{\phi}^{*})$ associated with the parameters in 
$\vec{\phi}^{*}$ shall be called the interim prior because in many cases it was 
our best guess as to the true redshift distribution before observing any data.  
For example, in template-based photo-$z$ PDF methods, the interim prior is a 
linear combination of the redshift distributions of the classes of templates 
based on previous observations.  (In training-based photo-$z$ PDF methods, the 
interim prior is related to the redshift distribution of the training set and 
thus may not be shared among all galaxies for which photo-$z$ PDFs are 
produced; however, in this work, we will assume that it is shared among all 
galaxies in a survey.)  Because the choice of parameters $\vec{\phi}^{*}$ is 
not causally connected to the data or the true redshift distribution but 
nonetheless contributes to the photo-$z$ PDF, we call it an \textit{interim 
prior}.  

A photo-$z$ PDF is then an \textit{interim posterior} probability distribution 
of the redshift of a galaxy given its observed photometric data and the interim 
prior: $p(z_{i} | \vec{d}_{i}, \vec{\phi}^{*})$.  
We present the two attributes of a probabilistic graphical model (PGM) for the 
redshift distribution $n(z)$: the directed acyclic graph (DAG) and its 
mathematical interpretation in terms of Bayesian hierarchical inference.  

The fundamental assumption underlying the concept of photo-$z$ estimation is 
that each galaxy $i$ has some observed photometric data $\vec{d}_{i}$ (fluxes, 
magnitudes, or colors) that is drawn from a function of its redshift $z_{i}$, 
which is a parameter in this model.  This function constitutes a forward model 
for the observations.  The redshifts $\{z_{i}\}$ for all galaxies in a survey 
are random draws from the redshift distribution, whose parameters  $\vec{\phi}$ 
under a chosen functional form are called \textit{hyperparameters} because they 
are shared among all $N$ galaxies in the survey.  Fig. \ref{fig:pgm} 
illustrates these relationships.  PGMs of this structure are 
\textit{hierarchical} in that, though the data are only directly influenced by 
their redshift parameters, we can still use them to infer something about the 
global hyperparameters by way of a higher level in the graph.

\begin{figure}
	\begin{center}We would like to learn about $p(z | \vec{\phi})$ from a 
catalog of photo-$z$ interim posteriors.
		
		\subsection{Model}
		\label{sec:model}
		
		\includegraphics[width=0.5\textwidth]{fig/pgm.png}
		\caption{This directed acyclic graph corresponds to a 
probabilistic graphical model for a Bayesian hierarchical inference of 
$p(\vec{\phi}|\{\vec{d}_{i}\})$.  In this graph, all random variables are shown 
in circles, with observed random variables shown in shaded circles.  
Relationships between variables are indicated by arrows from parameters to the 
variables distributed according to functions of them.  The box indicates that 
there are a number of copies of the relationships between boxed parameters, 
each independent of all others.  The hyperparameters comprising $\vec{\phi}$, 
which define $n(z)$, are at the top.  Independently drawn from a function of 
the hyperparameters $\vec{\phi}$ are galaxy redshifts $\{z_{i}\}$ below.  The 
observed galaxy photometry $\{\vec{d}_{i}\}$, shown in shaded circles, is 
determined by the redshifts above.}
	\label{fig:pgm}
	\end{center}
\end{figure}

The DAG of Fig. \ref{fig:pgm} translates directly into a mathematical 
expression for the posterior distribution $p(\vec{\phi} | \{\vec{d}_{i}\})$ of 
the hyperparameters given the entire set of interim posteriors $\{p(z_{i} | 
\vec{d}_{i}, \vec{\phi}^{*})\}$ by way of the following derivation.  From this 
point forward, we will use log probabilities.

According to Bayes' rule, the posterior distribution of interest is
\begin{align}
\label{eq:hyper_bayes}
\log[p(\vec{\phi} | \{\vec{d}_{i}\})] &\propto \log[p(\{\vec{d}_{i}\} | 
\vec{\phi})] + \log[p(\vec{\phi})] ,
\end{align}
where $p(\vec{\phi})$ is the \textit{hyperprior} probability distribution over 
possible values of the hyperparameters.  The hyperprior is a choice made by 
those performing the inference and is not to be confused with the interim prior 
parameters that influence the photo-$z$ PDF data product.

Focusing on the likelihood term $p(\{\vec{d}_{i}\} | \vec{\phi})$, we employ 
the power of the hierarchical model through \textit{marginalization} of the 
unobserved redshift parameters by integrating over them to obtain 
\begin{align}
\label{eq:marginalize}
\begin{split}
\log[p(\{\vec{d}_{i}\} | \vec{\phi})] = & 
\log\left[\int\exp\left[\log[p(\{\vec{d}_{i}\} | \{z_{i}\})]\right.\right. \\
& \left.\left.+ \log[p(\{z_{i}\} | \vec{\phi})]\right]\ d\{z_{i}\}\right] .
\end{split}
\end{align}
The redshifts $\{z_{i}\}$ are statistically independent of each other, and each 
galaxy's photometry depends solely on its own redshift, so we may invoke 
statistical independence to get
\begin{align}
\label{eq:independence_hyper}
\log[p(\{z_{i}\} | \vec{\phi})] &= \sum_{i}^{N} \log[p(z_{i} | \vec{\phi})]
\end{align}
and
\begin{align}
\label{eq:independence}
\log[p(\{\vec{d}_{i}\} | \{z_{i}\})] &= \sum_{i}^{N} \log[p(\vec{d}_{i} | 
z_{i})] .
\end{align}
Having chosen the functional form mapping $\vec{\phi}$ to $n(z)$, the terms on 
the righthand side of Eq. \ref{eq:independence_hyper} are known.  However, the 
terms on the righthand side of Eq. \ref{eq:independence} are likelihoods, and 
the photo-$z$ PDFs we have are interim posteriors.

To transform the interim posteriors we have into the likelihoods we need, we 
start with the vacuously true statement
\begin{align}
\begin{split}
\label{eq:unity}
\log[p(\vec{d}_{i} | z_{i})] =& \log[p(\vec{d}_{i} | z_{i})] + \log[p(z_{i} | 
\vec{d}_{i}, \vec{\phi}^{*})]\\
& - \log[p(z_{i} | \vec{d}_{i}, \vec{\phi}^{*})].
\end{split}
\end{align}
We expand the last term using Bayes' Rule to get
\begin{align}
\begin{split}
\label{eq:bayes}
\log[p(\vec{d}_{i} | z_{i})] =& \log[p(\vec{d}_{i} | z_{i})] + \log[p(z_{i} | 
\vec{d}_{i}, \vec{\phi}^{*})]\\
& - \log[p(\vec{d}_{i} | z_{i}, \vec{\phi}^{*})]\\
& + \log[p(\vec{d}_{i} | \vec{\phi}^{*})] - \log[p(z_{i} | \vec{\phi}^{*})].
\end{split}
\end{align}
Fig. \ref{fig:pgm} has no direct connection between $\vec{d}_{i}$ and 
$\vec{\phi}$, so we can break up $\log[p(\vec{d}_{i} | z_{i}, \vec{\phi}^{*})]$ 
as follows:
\begin{align}
\begin{split}
\label{eq:split}
\log[p(\vec{d}_{i} | z_{i})] =& \log[p(\vec{d}_{i} | z_{i})] + \log[p(z_{i} | 
\vec{d}_{i}, \vec{\phi}^{*})]\\
& - \log[p(\vec{d}_{i} | z_{i})] - \log[p(\vec{d}_{i} | \vec{\phi}^{*})]\\
& + \log[p(\vec{d}_{i} | \vec{\phi}^{*})] - \log[p(z_{i} | \vec{\phi}^{*})].
\end{split}
\end{align}
Canceling terms, we find
\begin{align}
\label{eq:cancel}
\log[p(\vec{d}_{i} | z_{i})] =& \log[p(z_{i} | \vec{d}_{i}, \vec{\phi}^{*})] - 
\log[p(z_{i} | \vec{\phi}^{*})].
\end{align}

Now we combine Eqs. \ref{eq:hyper_bayes}, \ref{eq:marginalize}, 
\ref{eq:independence_hyper}, \ref{eq:independence}, and \ref{eq:cancel} to 
arrive at 
\begin{align}
\begin{split}
\label{eq:final}
\log[p(\vec{\phi} | \{\vec{d}_{i}\})] \propto& \log[p(\vec{\phi})] + 
\log\left[\int \exp\left[\sum_{i}^{N} \left(\right.\right.\right.\\
& \left.\left.\left. + \log[p(z_{i} | \vec{\phi})] - \log[p(z_{i} | 
\vec{\phi}^{*})]\right.\right.\right.\\
& \left.\left.\left. + \log[p(z_{i} | \vec{d}_{i}, \vec{\phi}^{*})]
\right)\right] d\{z_{i}\}\right].
\end{split}
\end{align}
Thus, if we have a catalog of photo-$z$ interim posteriors $\{p(z_{i} | 
\vec{d}_{i}, \vec{\phi}^{*})\}$, the interim prior parameters $\vec{\phi}^{*}$, 
and a prior distribution $p(\vec{\phi})$, we may find the posterior 
$\log[p(\vec{\phi} | \{\vec{d}_{i}\})]$ on the redshift distribution $n(z)$.

This framework entails a number of choices and assumptions that must be 
addressed explicitly:
\begin{enumerate}
	\item The chosen functional form of $n(z)$ with parameters $\vec{\phi}$ 
must be capable of describing the true redshift distribution.
	\item The photometric data and redshift for each galaxy must be 
independent of the photometric data and redshifts of all other galaxies.
	\item There is one interim prior shared among all galaxies in the 
survey, and its parameters with known parameters $\vec{\phi}^{*}$ are known.
	\item We must choose a hyperprior probability distribution 
$p(\vec{\phi})$ over the hyperparameters.
	\item While we advocate for the approach of hierarchical inference, the 
probabilistic graphical model presented here is not the only one that could be 
proposed.
\end{enumerate}

\subsection{Alternative Approaches}
\label{sec:others}

It is useful to translate some popular existing methods for deriving $n(z)$ 
from photo-$z$ PDFs into the mathematical framework of Sec. \ref{sec:model}.  
We briefly discuss the conditions under which they are equivalent to Eq. 
\ref{eq:final}, which are elaborated upon in the Appendix.

\subsubsection{Stacking}
\label{sec:stacking}

The most common way to combine photo-$z$ interim posteriors into an estimator 
$\hat{n}_{stack}(z)$ for $n(z)$ is to "stack" them, which corresponds to
\begin{align}
\label{eq:stack}
\hat{n}^{stack}(z) &= \frac{1}{N}\ \sum_{i}^{N} \exp\left[\log[p(z_{i} | 
\vec{d}_{i}, \vec{\phi}^{*})]\right].
\end{align}

\subsubsection{Point Estimation}
\label{sec:point_estimation}

In some cases, estimators $\hat{n}^{point}(z)$ are obtained by assuming the 
photo-$z$ PDFs are effectively delta functions by reducing them to point 
estimates according to
\begin{align}
\label{eq:delta}
p(z_{i} | \vec{d}_{i}, \vec{\phi}^{*}) &\approx \delta(z, \hat{z}_{i}^{point})
\end{align}
before applying Eq. \ref{eq:stack}.  The most popular redshift point estimators 
are the mean $\hat{z}^{mean}_{i}$, median $\hat{z}^{median}_{i}$, and mode 
$\hat{z}^{mode}_{i}$ of the original photo-$z$ interim posterior.

\subsection{Implementation}
\label{sec:implementation}

The publicly available \chippr code implements the probabilistic graphical 
model presented in Sec. \ref{sec:model}.  

In addition to the choices and assumptions underlying the probabilistic 
graphical model, the implementation of \chippr makes choices and assumptions of 
its own.

\begin{enumerate}
	\item \chippr is only applicable if the redshift interim posteriors and 
interim prior are accurate.
	\item \chippr currently only accepts photo-$z$ PDFs and produces $n(z)$ 
samples of a single format, that of the piecewise constant parametrization, 
also referred to as a binned histogram parametrization and a sum of top hat 
functions.
\end{enumerate}

\section{Validation on Simple Mock Data}
\label{sec:validation}

We demonstrate the superiority of \chippr over alternative approaches in a 
number of compelling test cases on mock data.  Each experiment is characterized 
by a single change to a fiducial case in order to isolate the influence of 
systematic effects known to be relevant to photo-$z$ estimation and propagation 
in analysis.  

\subsection{Mock Data \& Metrics}
\label{sec:validintro}

The mock data in these tests consists of photo-$z$ interim posteriors rather 
than photometric data because the various existing methods for deriving 
photo-$z$ interim posteriors do not in general yield results that are 
consistent with one another, indicating that their systematics are not 
well-understood.  Because the mock data is independent of any choice of 
photo-$z$ PDF production method, we not only ensure that our photo-$z$ interim 
priors are perfectly understood but also deter readers from assuming that 
\chippr has any preference over the method by which photo-$z$ interim 
posteriors are derived from photometric data.

\subsubsection{Mock Data}
\label{sec:mockdata}

The mock data used here are produced by the following steps.

\begin{enumerate}
	\item \label{itm:true} Choose a true $n(z)$ that is a continuous 
function with known parameters $\vec{\phi}'$.
	\item \label{itm:sample} Sample the true $n(z)$ to create a catalog of 
$N$ true redshifts $z'_{i}$.
	\item \label{itm:scatter} Choose a true intrinsic scatter parameter 
$\sigma$, and sample Gaussians $p(z''_{i} | z'_{i}, \sigma) = 
\mathcal{N}(z'_{i}, \sigma)$ to get "observed" redshifts $z'_{i}$ defining 
Gaussian likelihoods $p(\vec{d}_{i} | z''_{i}, \sigma) = \mathcal{N}(z''_{i}, 
\sigma)$.
	\item \label{itm:interim} Choose a parametrization and the parameters 
$\vec{\phi}^{*}$ of the interim prior $p(z_{i} | \vec{\phi}^{*})$.
	\item \label{itm:posterior} Create interim posteriors $\log[p(z_{i} | 
\vec{d}_{i}, \vec{\phi}^{*})] = \log[p(\vec{d}_{i} | z'_{i})]\ \log[p(z_{i} | 
\vec{\phi}^{*})]$ in this parametrization.
\end{enumerate}

In all of the following validation tests, we use a piecewise constant 
parametrization in log-space with $10$ bins.  $N = 10,000$ galaxies.  In Secs. 
\ref{sec:truth} and \ref{sec:likelihoods}, we use a flat interim prior.  This 
method for deriving mock data is referred to as the fiducial case, and 
variations on it will refer directly to the steps that are altered.

\subsubsection{$n(z)$ Accuracy Metric}
\label{sec:accuracy}

The Kullback-Leibler divergence is our primary measure of the accuracy of 
estimators of $n(z)$ in cases of mock data with known true redshifts.

\textbf{Review precision and bias from \texttt{kld.ipynb} and interpret in 
terms of a \% difference.}

\subsection{Underlying $n(z)$ Effects}
\label{sec:truth}

Existing $n(z)$ estimators are systematically smoother than the true $n(z)$.  
Here we show that the traditional estimators perform better when the true 
$n(z)$ is weakly featured than when the true $n(z)$ is strongly featured by 
experimenting with Step \ref{itm:true} of Sec. \ref{sec:mockdata}.  The 
implication of this issue is quite serious; the consistently smooth, unimodal 
$n(z)$ estimates appearing in the literature could result from much more 
featured true redshift distributions, and there would be no way to catch this 
error without using a fully probabilistic method.

\subsubsection{Featureless $n(z)$}
\label{sec:smooth}

In this test, we choose a weakly featured true $n(z)$ of Fig. \ref{fig:smooth} 
with a smooth, unimodal shape, based on the interim prior used for the SDSS DR8 
photo-$z$ PDFs.  [check and include citation!]

\begin{figure}
	\begin{center}
		
\includegraphics[width=0.5\textwidth]{fig/smooth_truth/scatter.png}\\
		
\includegraphics[width=0.5\textwidth]{fig/smooth_truth/estimators.png}	
		\caption{All estimators perform well when the true $n(z)$ is 
well behaved, exhibiting significant deviation only when $n(z)$ is very small, 
as the sample size of true redshifts in that range will be small.  Top panel: 
The traditional MAP reduction of a photo-$z$ PDF against the true redshifts 
with a few rescaled photo-$z$ interim posteriors are overplotted in solid 
lines, with a dotted line indicating zero probability.  Bottom panel: Various 
estimators of $\ln[n(z)]$, the interim prior, and the true $\ln[n(z)]$ as a 
continuous function and under a binned parametrization.}
		\label{fig:featured}
	\end{center}
\end{figure}

\subsubsection{Featured $n(z)$}
\label{sec:featured}

In this test, we choose the true $n(z)$ of Fig. \ref{fig:featured} with 
nontrivial structure.

\begin{figure}
	\begin{center}
		
\includegraphics[width=0.5\textwidth]{fig/fiducial/scatter.png}\\
		
\includegraphics[width=0.5\textwidth]{fig/fiducial/estimators.png}		
		\caption{[This plot doesn't really show what I want because the 
intrinsic scatter is too low!  I used $\sigma=0.03$ because that's what has 
been quoted as what LSST, etc. needs, but that's not what comes out of 
photo-$z$ estimation methods before they impose aggressive cuts\dots]  Top 
panel: The traditional MAP reduction of a photo-$z$ PDF against the true 
redshifts with a few rescaled photo-$z$ interim posteriors are overplotted in 
solid lines, with a dotted line indicating zero probability.  Bottom panel: 
Various estimators of $\ln[n(z)]$, the interim prior, and the true $\ln[n(z)]$ 
as a continuous function and under a binned parametrization.}
		\label{fig:featured}
	\end{center}
\end{figure}

This featured true $n(z)$ will henceforth be referred to as the fiducial $n(z)$

\subsection{Emulated Data Quality Effects}
\label{sec:likelihoods}

In the following test cases, we vary the properties of the mock photo-$z$ 
likelihoods in an effort to emulate known systematics in photo-$z$ estimation.  
These tests vary Step \ref{itm:scatter} of Sec. \ref{sec:mockdata}.

\subsubsection{Intrinsic Scatter}
\label{sec:intscat}

One major concern about photo-$z$s is the intrinsic scatter of point 
estimators, including those derived from photo-$z$ PDFs, that is observed to 
varying extents with every existing photo-$z$ algorithm and illustrated in Fig. 
\ref{fig:intscat}.  To emulate intrinsic scatter, we modify the fiducial case 
to simply broaden the single Gaussian component of the likelihood.  To enforce 
self-consistency, the mean is a random variable drawn from a Gaussian 
distribution with the newly increased variance.

\begin{figure}
	\begin{center}
		
\includegraphics[width=0.5\textwidth]{fig/high_scatter/scatter.png}\\
		
\includegraphics[width=0.5\textwidth]{fig/high_scatter/estimators.png}		
		\caption{As the intrinsic scatter increases, the discrepancy 
between estimators increases.  In particular, the stacked estimator and 
marginalized point estimators predict $\ln[n(z)]$ to be smoother than the 
truth, while the  [This would be a lot more compelling with more galaxies.  
Also, the weird edge effects in the top panel are real, because the point 
estimator is the MAP, not the center of the Gaussian likelihood, and there's no 
requirement that the mean of the likelihood be within the true redshift range.] 
 Top panel: The traditional MAP reduction of a photo-$z$ PDF against the true 
redshifts with a few rescaled photo-$z$ interim posteriors are overplotted in 
solid lines, with a dotted line indicating zero probability.  Bottom panel: 
Various estimators of $\ln[n(z)]$, the interim prior, and the true $\ln[n(z)]$ 
as a continuous function and under a binned parametrization.}
		\label{fig:intscat}
	\end{center}
\end{figure}

To emulate intrinsic scatter, we modify the fiducial case to simply broaden the 
single Gaussian component of the likelihood.  To enforce self-consistency, the 
mean is a random variable drawn from a Gaussian distribution with the newly 
increased variance.

\subsubsection{Template-like Catastrophic Outliers}
\label{sec:tempcatout}

In addition to intrinsic scatter, photo-$z$ methods employing template fitting 
tend to produce catastrophic outliers that are distributed to be broad in 
$z_{spec}$ and narrow in $z_{phot}$, as in Fig. \ref{fig:tempcatout}.  The 
systematic behind these catastrophic outliers may be described as an attractor 
in the space of $z_{phot}$; some galaxies at a range of $z_{spec}$ map onto a 
single $z_{phot}$ (with some scatter) if their true SED does not have 
sufficiently strong features (as is the case for blue galaxies), leading 
galaxies of that SED type at many $z_{spec}$ to have similar colors.

\begin{figure}
	\begin{center}
		
\includegraphics[width=0.5\textwidth]{fig/template_outliers/scatter.png}\\
		
\includegraphics[width=0.5\textwidth]{fig/template_outliers/estimators.png}	
	
		\caption{Top panel: The traditional MAP reduction of a 
photo-$z$ PDF against the true redshifts with a few rescaled photo-$z$ interim 
posteriors are overplotted in solid lines, with a dotted line indicating zero 
probability.  Bottom panel: Various estimators of $\ln[n(z)]$, the interim 
prior, and the true $\ln[n(z)]$ as a continuous function and under a binned 
parametrization.}
		\label{fig:tempcatout}
	\end{center}
\end{figure}

\subsubsection{Training-like Catastrophic Outliers}
\label{sec:traincatout}

Data driven photo-$z$ methods tend to suffer from a different form of 
catastrophic outliers that are distributed to be narrow in $z_{spec}$ and broad 
in $z_{phot}$, as in Fig. \ref{fig:traincatout}.  The systematic behind these 
catastrophic outliers may be described as an attractor in the space of 
$z_{spec}$; some galaxies near a particular $z_{spec}$ map to a range of 
$z_{phot}$ if the training set galaxies at that $z_{spec}$ have inconsistent 
$z_{phot}$, as might occur if their SED's features fall between photometric 
filters.

\begin{figure}
	\begin{center}
		
\includegraphics[width=0.5\textwidth]{fig/training_outliers/scatter.png}\\
		
\includegraphics[width=0.5\textwidth]{fig/training_outliers/estimators.png}	
		\caption{[This has proven to be the most challenging test case 
to implement, and there's clearly still a bug in the function that makes the 
likelihoods.]}
		\label{fig:traincatout}
	\end{center}
\end{figure}

\subsection{Emulated Interim Prior Effects}
\label{sec:priors}

The interim prior encapsulates the the relationship between observed photometry 
and redshift information upon which a photo-$z$ estimate is based.  Interim 
priors are in general not identical to the true $n(z)$ we wish to estimate; if 
they were, we would not need any data!  For template fitting photo-$z$ methods, 
the interim prior is usually an input chosen by the researcher.  However, for 
machine learning methods, the interim prior is some function of the training 
set data that in many cases may be influenced by random numbers and is rarely 
output with the redshift estimates.  Interim priors for template fitting 
methods tend to have incomplete coverage in the space of true photometry, 
because they are limited by the choice of the library of SEDs.  Interim priors 
for machine learning methods tend to have incomplete coverage in the space of 
redshifts, because there are fewer galaxies with spectroscopically confirmed 
redshifts at high redshifts than low redshifts.

Existing $n(z)$ estimation routines will always produce a biased estimator when 
the interim prior is not equal to the true $n(z)$.  We demonstrate here that 
regardless of the appropriateness of the interim prior as an approximation to 
the true $n(z)$, \chippr is not affected by the choice of the interim prior so 
long as it has nontrivial coverage in the space of redshift.  These tests 
modify Step \ref{itm:interim} of Sec. \ref{sec:mockdata}.

\subsubsection{Template-like Interim Prior}
\label{sec:tempintpr}

An interim prior based on a template library may be a sum of smooth functions 
representing $n(z)$ for each SED type in the library.  Template libraries do 
not include every possible galaxy SED, and the $n(z)$ used for each SED type 
may not be accurate.  The interim prior shown in Fig. \ref{fig:tempintpr} is an 
emulation of an interim prior corresponding to a template library of this type.

\begin{figure}
	\begin{center}
		
\includegraphics[width=0.5\textwidth]{fig/template_prior/scatter.png}\\
		
\includegraphics[width=0.5\textwidth]{fig/template_prior/estimators.png}	
		\caption{[The case of a multimodal interim prior was a very 
compelling test in the previous version but somehow isn't anymore.]}
		\label{fig:tempintpr}
	\end{center}
\end{figure}

\subsubsection{Training-like Interim Prior}
\label{sec:trainintpr}

An interim prior based on a training set may be biased toward low redshifts due 
to the dearth of distant galaxies with spectroscopic redshifts.  The interim 
prior shown in Fig. \ref{fig:trainintpr} is an emulation of an interim prior 
corresponding to a training set biased in this way.  We chose an interpolation 
of the interim prior used for the SDSS DR8 photo-$z$ PDFs.  

\begin{figure}
	\begin{center}
		
\includegraphics[width=0.5\textwidth]{fig/training_prior/scatter.png}\\
		
\includegraphics[width=0.5\textwidth]{fig/training_prior/estimators.png}	
		\caption{[The case of a low-$z$ favoring interim prior was a 
very compelling test in the previous version but somehow isn't anymore.]}
		\label{fig:trainintpr}
	\end{center}
\end{figure}

\section{Application to Realistic Mock Data}
\label{sec:application}

To show how the choice of $n(z)$ estimator propagates to cosmological 
constraints, we apply \chippr to a data from a realistic cosmological 
simulation (probably \texttt{Buzzard}) with photo-$z$ PDFs produced by a 
popular method (probably \texttt{BPZ}).

\subsection{Mock Data \& Metrics}
\label{sec:appintro}

As in Sec. \ref{sec:validintro}, the mock data takes the form of photo-$z$ 
interim posteriors, but that is where the similarity ends.  These photo-$z$ 
interim posteriors are derived from the photometry resulting from the 
\texttt{Buzzard} simulation by way of \texttt{BPZ}.  Because the simulation 
begins with setting true values of the cosmological parameters, we can 
propagate the different estimators of $n(z)$ through a forecasting code (which 
one?) to generate error ellipses on the cosmological parameters.

\subsubsection{Mock Data}
\label{sec:buzzard}

We summarize the details of the \texttt{Buzzard} simulation here.

\subsubsection{Cosmological Constraint Metric}
\label{sec:cosmo}

Because the realistic mock data of Sec. \ref{sec:buzzard} is associated with 
true values of the cosmological parameters, we may compare the quality of 
cosmological constraints under different estimators of $n(z)$, creating a 
figure like Fig. \ref{fig:scgsr}.  We perform this analysis using a forecasting 
code that takes $n(z)$ and produces projected error ellipses in the space of 
cosmological parameters.  

\begin{figure}
	\begin{center}
		Not my figure! (<NAME> via <NAME>)
		\includegraphics[width=0.5\textwidth]{cosmo_constraints.png}
		\caption{[I'd like to be able to make a plot like this with the 
different $n(z)$ estimators produced by \texttt{chippr} from the same set of 
photo-$z$ PDFs.  Then I can calculate some metrics of the accuracy and 
precision of the error distributions relative to the true values of the 
cosmological parameters that produced the mock data.]}
		\label{fig:scgsr}
	\end{center}
\end{figure}


\subsection{Results}
\label{sec:results}

\begin{figure}
	\begin{center}
		\caption{[moneyplot of error ellipses resulting from different 
$n(z)$ estimators]}
		\label{fig:money}
	\end{center}
\end{figure}

\section{Discussion}
\label{sec:discussion}

The results of Fig \ref{fig:money} have significant implications for the 
developing data analysis pipelines of next-generation telescope surveys.  
However, the method presented here has its own limitations, which are 
reiterated to discourage the community from applying this work inappropriately.

We intend to pursue a number of extensions of the work presented in this paper 
in future work.

\section{Conclusion}
\label{sec:conclusion}

We now summarize answers to the questions posed in the introduction:

\begin{itemize}
	\item Existing $n(z)$ estimation methods produce biased estimators that 
propagate to inaccuracies in characterizing the cosmological parameters.
	\item Photo-$z$ PDFs are probabilistic data products so must be handled 
in a mathematically consistent manner such as the probabilistic graphical model 
outlined in this paper.
	\item In addition to coming with its own error distribution, the $n(z)$ 
estimator produced by \chippr is quantifiably more accurate than established 
estimators.
	\item Propagation of the \chippr result leads to a quantifiable 
improvement in the constraints on cosmological parameters.
\end{itemize}

In conclusion, we discourage the community from continuing to use the stacked 
estimator and reductions of photo-$z$ PDFs to redshift point estimates in 
obtaining estimators of $n(z)$.  \chippr is freely available to the community 
for incorporation into evolving data analysis pipelines.  


\appendix{Catalog Production}
\label{app:catalog}

Sec. \ref{sec:mockdata} outlines one way to make mock data, but it's not 
necessarily a logically consistent way to do it, a problem identified in trying 
to use it for the test cases in Sec. \ref{sec:likelihoods}.  If the goal is to 
emulate the traditional $z_{spec}$ vs. $z_{phot}$ scatterplots, we need to 
start from probability distributions in that space.  This means there are two 
parameters $z_{spec}$ and $z_{phot}$ drawn from a two-dimensional distribution 
$p(z_{spec}, z_{phot} | \textul{\sigma}, \vec{\phi})$, where $\textul{\sigma}$ 
contains parameters concerning the relationship between $z_{spec}$ and 
$z_{phot}$, and $\vec{\phi}$ contains hyperparameters concerning $n(z)$.  We 
know from basic probability that
\begin{align}
\label{eq:trivial}
p(z_{spec}, z_{phot} | \textul{\sigma}, \vec{\phi}) =& p(z_{phot} | z_{spec}, 
\textul{\sigma}, \vec{\phi})\ p(z_{sped} | \textul{\sigma}, \vec{\phi}).
\end{align}
Because $z_{spec}$ depends only on $\vec{\phi}$ and $z_{phot}$ is independent 
of $\vec{\phi}$, this becomes
\begin{align}
\label{eq:simpler}
p(z_{spec}, z_{phot} | \textul{\sigma}, \vec{\phi}) =& p(z_{phot} | z_{spec}, 
\textul{\sigma})\ p(z_{spec} | \vec{\phi}).
\end{align}
We know how to make reasonable choices for both terms on the righthand side of 
Eq. \ref{eq:simpler}.

In the $z_{spec}$ vs. $z_{phot}$ scatterplot, $p(z_{phot} | z_{spec}, 
\textul{\sigma}, \vec{\phi})$ would be represented by vertical slices, but we 
are interested in the horizontal slices that represent $p(z_{spec} | z_{phot}, 
\textul{\sigma}, \vec{\phi})$.  Now we must distinguish the true 
hyperparameters $\vec{\phi}_{true}$ from the interim hyperparameters 
$\vec{\phi}_{int}$.  If the horizontal slices carry information about 
$\vec{\phi}$, how can we get a catalog of $p(z_{spec} | z_{phot}, 
\textul{\sigma}, \vec{\phi}_{int})$ that's consistent with $p(z_{spec} | 
z_{phot}, \textul{\sigma}, \vec{\phi}_{true})$?  Somehow that information needs 
to be propagated through in order to recover it.

We achieve this goal by dividing out the true $n(z)$ and multiplying by the 
interim $n(z)$ as in
\begin{align}
\label{eq:transformation}
p(z_{spec} | z_{phot}, \textul{\sigma}, \vec{\phi}_{int}) =& p(z_{spec} | 
z_{phot}, \textul{\sigma}, \vec{\phi}_{true})\ \frac{p(z_{spec} | 
\vec{\phi}_{int})}{p(z_{spec} | \vec{\phi}_{true})},
\end{align}
since $p(z_{spec} | z_{phot}, \textul{\sigma}, \vec{\phi})$ is separable and 
the term $p(z_{spec} | z_{phot}, \textul{\sigma})$ changes independently of the 
$p(z_{spec} | \vec{\phi})$ term.  The information about $\vec{\phi}_{true}$ is 
retained by the elements of the catalog: only $(z_{spec}, z_{phot})$ pairs 
drawn from the true distribution will be represented in the catalog!

It appears as if this is the only fully self-consistent way to ensure that a 
draw from the posterior distribution $p(z_{spec} | z_{phot}, \textul{\sigma}, 
\vec{\phi})$ does in fact follow that distribution.  In the cases with fiducial 
posteriors $p(z_{spec} | z_{phot}, \sigma, \vec{\phi})$, the method of Sec. 
\ref{sec:mockdata} was equivalent to this, but Sec. \ref{sec:likelihoods} 
requires a more sophisticated approach.

\begin{acknowledgements}
AIM thanks <NAME> for insightful input on statistics, Geoffrey 
Ryan for assistance in debugging, and <NAME> for helpful comments 
provided in the preparation of this paper.
\end{acknowledgements}


\end{document}<commit_msg>updated abstract<commit_after>\documentclass[iop]{emulateapj}

\usepackage{tikz}
\usepackage{natbib}
\usepackage{amsmath}

\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]

\newcommand{\myemail}{<EMAIL>}
\newcommand{\textul}{\underline}
\newcommand{\chippr}{\texttt{CHIPPR} }

\shorttitle{How to obtain the redshift distribution from probabilistic redshift 
estimates}
\shortauthors{Malz, et al.}

\begin{document}

\title{How to obtain the redshift distribution from probabilistic redshift 
estimates}

\author{<NAME>\altaffilmark{1}, <NAME>\altaffilmark{1,2,3,4}, Phil 
Marshall\altaffilmark{5}}
\email{<EMAIL>}

\altaffiltext{1}{Center for Cosmology and Particle Physics, Department of 
Physics,
  New York University, 726 Broadway, 9th floor, New York, NY 10003, USA}
\altaffiltext{2}{Simons Center for Computational Astrophysics, 162 Fifth 
Avenue, 7th floor, New York, NY 10010, USA}
\altaffiltext{3}{Center for Data Science, New York University, 60 Fifth Avenue, 
7th floor, New York, NY 10003, USA}
\altaffiltext{4}{Max-Planck-Institut f\"ur Astronomie, K\"onigstuhl 17, D-69117 
Heidelberg, Germany}
\altaffiltext{5}{[SLAC]}

\begin{abstract}
A trustworthy estimate of the redshift distribution $n(z)$ is crucial for weak 
lensing cosmology.  Spectroscopically confirmed redshifts for the dim and 
numerous galaxies of weak lensing surveys are expected to be inaccessible, 
making photometric redshifts (photo-$z$s) the next best alternative.  The 
nontrivial systematics affecting photo-$z$ estimation have motivated the weak 
lensing community to favor photo-$z$ probability distribution functions (PDFs) 
as a more comprehensive alternative to photo-$z$ point estimates.  However, 
analytic methods for utilizing these new data products in cosmological 
inference are still evolving.  This paper presents the Cosmological 
Hierarchical Inference with Probabilistic Photometric Redshifts (CHIPPR) model, 
a novel approach to estimating the $n(z)$ from a catalog of galaxy photo-$z$ 
PDFs based upon a probabilistic graphical model of hierarchical inference.  We 
present the publicly available \chippr code implementing this technique, as 
well as its validation on mock data.  The CHIPPR model yields a more accurate 
characterization of $n(z)$ containing information beyond the best-fit estimator 
produced by traditional procedures.

\end{abstract}

\keywords{cosmology: cosmological parameters --- galaxies: statistics --- 
gravitational lensing: weak --- methods: data analysis --- methods: statistical}

\section{Introduction}
\label{sec:introduction}

After a brief literature review addressing how photo-$z$ PDFs are currently 
used in cosmology, this paper aims to answer the following questions:

\begin{itemize}
	\item Why should we question existing methods?
	\item How can we improve the effectiveness of using photo-$z$ PDFs in 
inference?
	\item How does the result of \chippr compare to established estimators 
in terms of the accuracy of $n(z)$?
	\item How significant is the effect of the discrepancy between $n(z)$ 
estimators on cosmological constraints?
\end{itemize}

\section{Method}
\label{sec:method}

This paper presents a mathematically consistent method for obtaining the 
posterior probability distribution over the redshift distribution $n(z)$ using 
a catalog of photo-$z$ PDFs.  We start by introducing some nomenclature, 
definitions, and symbols that will be used throughout Secs. \ref{sec:model}, 
\ref{sec:others}, and \ref{sec:implementation}.

We shall say the parameters comprising $\vec{\phi}$ define the redshift 
distribution under some functional form about which this method is agnostic, so 
long as the functional form evaluated with some true values of the parameters 
accurately describes the true $n(z)$.  Since the redshift distribution is 
itself a probability distribution, it may be written as $p(z | \vec{\phi})$.

Before jumping into the details of the model, it is crucial to settle on an 
interpretation of what a photo-$z$ PDF actually is.  When a photo-$z$ PDF is 
reported, it is an object containing information about the redshift $z_{i}$ of 
galaxy $i$ based on its photometric data $\vec{d}_{i}$, making a photo-$z$ PDF 
a \textit{posterior} probability distribution.  The data $\vec{d}_{i}$ may be 
fluxes, magnitudes, colors, or any combination thereof.  We assume that the 
redshifts of the galaxies are not related to the redshifts or photometry of 
other galaxies in the survey, i.e. they are \textit{independent}.

In addition to the data $\vec{d}_{i}$, the photo-$z$ PDF also contains 
information imparted by assumptions that went into the process by which the 
photo-$z$ was made.  Those assumptions are only relevant insomuch as they 
impose a preference on the redshift distribution. Thus, they can be reduced to 
a sort of prior on the distribution of redshifts, which can be parametrized 
under the chosen functional form by parameters $\vec{\phi}^{*}$.  The redshift 
distribution $p(z | \vec{\phi}^{*})$ associated with the parameters in 
$\vec{\phi}^{*}$ shall be called the interim prior because in many cases it was 
our best guess as to the true redshift distribution before observing any data.  
For example, in template-based photo-$z$ PDF methods, the interim prior is a 
linear combination of the redshift distributions of the classes of templates 
based on previous observations.  (In training-based photo-$z$ PDF methods, the 
interim prior is related to the redshift distribution of the training set and 
thus may not be shared among all galaxies for which photo-$z$ PDFs are 
produced; however, in this work, we will assume that it is shared among all 
galaxies in a survey.)  Because the choice of parameters $\vec{\phi}^{*}$ is 
not causally connected to the data or the true redshift distribution but 
nonetheless contributes to the photo-$z$ PDF, we call it an \textit{interim 
prior}.  

A photo-$z$ PDF is then an \textit{interim posterior} probability distribution 
of the redshift of a galaxy given its observed photometric data and the interim 
prior: $p(z_{i} | \vec{d}_{i}, \vec{\phi}^{*})$.  
We present the two attributes of a probabilistic graphical model (PGM) for the 
redshift distribution $n(z)$: the directed acyclic graph (DAG) and its 
mathematical interpretation in terms of Bayesian hierarchical inference.  

The fundamental assumption underlying the concept of photo-$z$ estimation is 
that each galaxy $i$ has some observed photometric data $\vec{d}_{i}$ (fluxes, 
magnitudes, or colors) that is drawn from a function of its redshift $z_{i}$, 
which is a parameter in this model.  This function constitutes a forward model 
for the observations.  The redshifts $\{z_{i}\}$ for all galaxies in a survey 
are random draws from the redshift distribution, whose parameters  $\vec{\phi}$ 
under a chosen functional form are called \textit{hyperparameters} because they 
are shared among all $N$ galaxies in the survey.  Fig. \ref{fig:pgm} 
illustrates these relationships.  PGMs of this structure are 
\textit{hierarchical} in that, though the data are only directly influenced by 
their redshift parameters, we can still use them to infer something about the 
global hyperparameters by way of a higher level in the graph.

\begin{figure}
	\begin{center}We would like to learn about $p(z | \vec{\phi})$ from a 
catalog of photo-$z$ interim posteriors.
		
		\subsection{Model}
		\label{sec:model}
		
		\includegraphics[width=0.5\textwidth]{fig/pgm.png}
		\caption{This directed acyclic graph corresponds to a 
probabilistic graphical model for a Bayesian hierarchical inference of 
$p(\vec{\phi}|\{\vec{d}_{i}\})$.  In this graph, all random variables are shown 
in circles, with observed random variables shown in shaded circles.  
Relationships between variables are indicated by arrows from parameters to the 
variables distributed according to functions of them.  The box indicates that 
there are a number of copies of the relationships between boxed parameters, 
each independent of all others.  The hyperparameters comprising $\vec{\phi}$, 
which define $n(z)$, are at the top.  Independently drawn from a function of 
the hyperparameters $\vec{\phi}$ are galaxy redshifts $\{z_{i}\}$ below.  The 
observed galaxy photometry $\{\vec{d}_{i}\}$, shown in shaded circles, is 
determined by the redshifts above.}
	\label{fig:pgm}
	\end{center}
\end{figure}

The DAG of Fig. \ref{fig:pgm} translates directly into a mathematical 
expression for the posterior distribution $p(\vec{\phi} | \{\vec{d}_{i}\})$ of 
the hyperparameters given the entire set of interim posteriors $\{p(z_{i} | 
\vec{d}_{i}, \vec{\phi}^{*})\}$ by way of the following derivation.  From this 
point forward, we will use log probabilities.

According to Bayes' rule, the posterior distribution of interest is
\begin{align}
\label{eq:hyper_bayes}
\log[p(\vec{\phi} | \{\vec{d}_{i}\})] &\propto \log[p(\{\vec{d}_{i}\} | 
\vec{\phi})] + \log[p(\vec{\phi})] ,
\end{align}
where $p(\vec{\phi})$ is the \textit{hyperprior} probability distribution over 
possible values of the hyperparameters.  The hyperprior is a choice made by 
those performing the inference and is not to be confused with the interim prior 
parameters that influence the photo-$z$ PDF data product.

Focusing on the likelihood term $p(\{\vec{d}_{i}\} | \vec{\phi})$, we employ 
the power of the hierarchical model through \textit{marginalization} of the 
unobserved redshift parameters by integrating over them to obtain 
\begin{align}
\label{eq:marginalize}
\begin{split}
\log[p(\{\vec{d}_{i}\} | \vec{\phi})] = & 
\log\left[\int\exp\left[\log[p(\{\vec{d}_{i}\} | \{z_{i}\})]\right.\right. \\
& \left.\left.+ \log[p(\{z_{i}\} | \vec{\phi})]\right]\ d\{z_{i}\}\right] .
\end{split}
\end{align}
The redshifts $\{z_{i}\}$ are statistically independent of each other, and each 
galaxy's photometry depends solely on its own redshift, so we may invoke 
statistical independence to get
\begin{align}
\label{eq:independence_hyper}
\log[p(\{z_{i}\} | \vec{\phi})] &= \sum_{i}^{N} \log[p(z_{i} | \vec{\phi})]
\end{align}
and
\begin{align}
\label{eq:independence}
\log[p(\{\vec{d}_{i}\} | \{z_{i}\})] &= \sum_{i}^{N} \log[p(\vec{d}_{i} | 
z_{i})] .
\end{align}
Having chosen the functional form mapping $\vec{\phi}$ to $n(z)$, the terms on 
the righthand side of Eq. \ref{eq:independence_hyper} are known.  However, the 
terms on the righthand side of Eq. \ref{eq:independence} are likelihoods, and 
the photo-$z$ PDFs we have are interim posteriors.

To transform the interim posteriors we have into the likelihoods we need, we 
start with the vacuously true statement
\begin{align}
\begin{split}
\label{eq:unity}
\log[p(\vec{d}_{i} | z_{i})] =& \log[p(\vec{d}_{i} | z_{i})] + \log[p(z_{i} | 
\vec{d}_{i}, \vec{\phi}^{*})]\\
& - \log[p(z_{i} | \vec{d}_{i}, \vec{\phi}^{*})].
\end{split}
\end{align}
We expand the last term using Bayes' Rule to get
\begin{align}
\begin{split}
\label{eq:bayes}
\log[p(\vec{d}_{i} | z_{i})] =& \log[p(\vec{d}_{i} | z_{i})] + \log[p(z_{i} | 
\vec{d}_{i}, \vec{\phi}^{*})]\\
& - \log[p(\vec{d}_{i} | z_{i}, \vec{\phi}^{*})]\\
& + \log[p(\vec{d}_{i} | \vec{\phi}^{*})] - \log[p(z_{i} | \vec{\phi}^{*})].
\end{split}
\end{align}
Fig. \ref{fig:pgm} has no direct connection between $\vec{d}_{i}$ and 
$\vec{\phi}$, so we can break up $\log[p(\vec{d}_{i} | z_{i}, \vec{\phi}^{*})]$ 
as follows:
\begin{align}
\begin{split}
\label{eq:split}
\log[p(\vec{d}_{i} | z_{i})] =& \log[p(\vec{d}_{i} | z_{i})] + \log[p(z_{i} | 
\vec{d}_{i}, \vec{\phi}^{*})]\\
& - \log[p(\vec{d}_{i} | z_{i})] - \log[p(\vec{d}_{i} | \vec{\phi}^{*})]\\
& + \log[p(\vec{d}_{i} | \vec{\phi}^{*})] - \log[p(z_{i} | \vec{\phi}^{*})].
\end{split}
\end{align}
Canceling terms, we find
\begin{align}
\label{eq:cancel}
\log[p(\vec{d}_{i} | z_{i})] =& \log[p(z_{i} | \vec{d}_{i}, \vec{\phi}^{*})] - 
\log[p(z_{i} | \vec{\phi}^{*})].
\end{align}

Now we combine Eqs. \ref{eq:hyper_bayes}, \ref{eq:marginalize}, 
\ref{eq:independence_hyper}, \ref{eq:independence}, and \ref{eq:cancel} to 
arrive at 
\begin{align}
\begin{split}
\label{eq:final}
\log[p(\vec{\phi} | \{\vec{d}_{i}\})] \propto& \log[p(\vec{\phi})] + 
\log\left[\int \exp\left[\sum_{i}^{N} \left(\right.\right.\right.\\
& \left.\left.\left. + \log[p(z_{i} | \vec{\phi})] - \log[p(z_{i} | 
\vec{\phi}^{*})]\right.\right.\right.\\
& \left.\left.\left. + \log[p(z_{i} | \vec{d}_{i}, \vec{\phi}^{*})]
\right)\right] d\{z_{i}\}\right].
\end{split}
\end{align}
Thus, if we have a catalog of photo-$z$ interim posteriors $\{p(z_{i} | 
\vec{d}_{i}, \vec{\phi}^{*})\}$, the interim prior parameters $\vec{\phi}^{*}$, 
and a prior distribution $p(\vec{\phi})$, we may find the posterior 
$\log[p(\vec{\phi} | \{\vec{d}_{i}\})]$ on the redshift distribution $n(z)$.

This framework entails a number of choices and assumptions that must be 
addressed explicitly:
\begin{enumerate}
	\item The chosen functional form of $n(z)$ with parameters $\vec{\phi}$ 
must be capable of describing the true redshift distribution.
	\item The photometric data and redshift for each galaxy must be 
independent of the photometric data and redshifts of all other galaxies.
	\item There is one interim prior shared among all galaxies in the 
survey, and its parameters with known parameters $\vec{\phi}^{*}$ are known.
	\item We must choose a hyperprior probability distribution 
$p(\vec{\phi})$ over the hyperparameters.
	\item While we advocate for the approach of hierarchical inference, the 
probabilistic graphical model presented here is not the only one that could be 
proposed.
\end{enumerate}

\subsection{Alternative Approaches}
\label{sec:others}

It is useful to translate some popular existing methods for deriving $n(z)$ 
from photo-$z$ PDFs into the mathematical framework of Sec. \ref{sec:model}.  
We briefly discuss the conditions under which they are equivalent to Eq. 
\ref{eq:final}, which are elaborated upon in the Appendix.

\subsubsection{Stacking}
\label{sec:stacking}

The most common way to combine photo-$z$ interim posteriors into an estimator 
$\hat{n}_{stack}(z)$ for $n(z)$ is to "stack" them, which corresponds to
\begin{align}
\label{eq:stack}
\hat{n}^{stack}(z) &= \frac{1}{N}\ \sum_{i}^{N} \exp\left[\log[p(z_{i} | 
\vec{d}_{i}, \vec{\phi}^{*})]\right].
\end{align}

\subsubsection{Point Estimation}
\label{sec:point_estimation}

In some cases, estimators $\hat{n}^{point}(z)$ are obtained by assuming the 
photo-$z$ PDFs are effectively delta functions by reducing them to point 
estimates according to
\begin{align}
\label{eq:delta}
p(z_{i} | \vec{d}_{i}, \vec{\phi}^{*}) &\approx \delta(z, \hat{z}_{i}^{point})
\end{align}
before applying Eq. \ref{eq:stack}.  The most popular redshift point estimators 
are the mean $\hat{z}^{mean}_{i}$, median $\hat{z}^{median}_{i}$, and mode 
$\hat{z}^{mode}_{i}$ of the original photo-$z$ interim posterior.

\subsection{Implementation}
\label{sec:implementation}

The publicly available \chippr code implements the probabilistic graphical 
model presented in Sec. \ref{sec:model}.  

In addition to the choices and assumptions underlying the probabilistic 
graphical model, the implementation of \chippr makes choices and assumptions of 
its own.

\begin{enumerate}
	\item \chippr is only applicable if the redshift interim posteriors and 
interim prior are accurate.
	\item \chippr currently only accepts photo-$z$ PDFs and produces $n(z)$ 
samples of a single format, that of the piecewise constant parametrization, 
also referred to as a binned histogram parametrization and a sum of top hat 
functions.
\end{enumerate}

\section{Validation on Simple Mock Data}
\label{sec:validation}

We demonstrate the superiority of \chippr over alternative approaches in a 
number of compelling test cases on mock data.  Each experiment is characterized 
by a single change to a fiducial case in order to isolate the influence of 
systematic effects known to be relevant to photo-$z$ estimation and propagation 
in analysis.  

\subsection{Mock Data \& Metrics}
\label{sec:validintro}

The mock data in these tests consists of photo-$z$ interim posteriors rather 
than photometric data because the various existing methods for deriving 
photo-$z$ interim posteriors do not in general yield results that are 
consistent with one another, indicating that their systematics are not 
well-understood.  Because the mock data is independent of any choice of 
photo-$z$ PDF production method, we not only ensure that our photo-$z$ interim 
priors are perfectly understood but also deter readers from assuming that 
\chippr has any preference over the method by which photo-$z$ interim 
posteriors are derived from photometric data.

\subsubsection{Mock Data}
\label{sec:mockdata}

The mock data used here are produced by the following steps.

\begin{enumerate}
	\item \label{itm:true} Choose a true $n(z)$ that is a continuous 
function with known parameters $\vec{\phi}'$.
	\item \label{itm:sample} Sample the true $n(z)$ to create a catalog of 
$N$ true redshifts $z'_{i}$.
	\item \label{itm:scatter} Choose a true intrinsic scatter parameter 
$\sigma$, and sample Gaussians $p(z''_{i} | z'_{i}, \sigma) = 
\mathcal{N}(z'_{i}, \sigma)$ to get "observed" redshifts $z'_{i}$ defining 
Gaussian likelihoods $p(\vec{d}_{i} | z''_{i}, \sigma) = \mathcal{N}(z''_{i}, 
\sigma)$.
	\item \label{itm:interim} Choose a parametrization and the parameters 
$\vec{\phi}^{*}$ of the interim prior $p(z_{i} | \vec{\phi}^{*})$.
	\item \label{itm:posterior} Create interim posteriors $\log[p(z_{i} | 
\vec{d}_{i}, \vec{\phi}^{*})] = \log[p(\vec{d}_{i} | z'_{i})]\ \log[p(z_{i} | 
\vec{\phi}^{*})]$ in this parametrization.
\end{enumerate}

In all of the following validation tests, we use a piecewise constant 
parametrization in log-space with $10$ bins.  $N = 10,000$ galaxies.  In Secs. 
\ref{sec:truth} and \ref{sec:likelihoods}, we use a flat interim prior.  This 
method for deriving mock data is referred to as the fiducial case, and 
variations on it will refer directly to the steps that are altered.

\subsubsection{$n(z)$ Accuracy Metric}
\label{sec:accuracy}

The Kullback-Leibler divergence is our primary measure of the accuracy of 
estimators of $n(z)$ in cases of mock data with known true redshifts.

\textbf{Review precision and bias from \texttt{kld.ipynb} and interpret in 
terms of a \% difference.}

\subsection{Underlying $n(z)$ Effects}
\label{sec:truth}

Existing $n(z)$ estimators are systematically smoother than the true $n(z)$.  
Here we show that the traditional estimators perform better when the true 
$n(z)$ is weakly featured than when the true $n(z)$ is strongly featured by 
experimenting with Step \ref{itm:true} of Sec. \ref{sec:mockdata}.  The 
implication of this issue is quite serious; the consistently smooth, unimodal 
$n(z)$ estimates appearing in the literature could result from much more 
featured true redshift distributions, and there would be no way to catch this 
error without using a fully probabilistic method.

\subsubsection{Featureless $n(z)$}
\label{sec:smooth}

In this test, we choose a weakly featured true $n(z)$ of Fig. \ref{fig:smooth} 
with a smooth, unimodal shape, based on the interim prior used for the SDSS DR8 
photo-$z$ PDFs.  [check and include citation!]

\begin{figure}
	\begin{center}
		
\includegraphics[width=0.5\textwidth]{fig/smooth_truth/scatter.png}\\
		
\includegraphics[width=0.5\textwidth]{fig/smooth_truth/estimators.png}	
		\caption{All estimators perform well when the true $n(z)$ is 
well behaved, exhibiting significant deviation only when $n(z)$ is very small, 
as the sample size of true redshifts in that range will be small.  Top panel: 
The traditional MAP reduction of a photo-$z$ PDF against the true redshifts 
with a few rescaled photo-$z$ interim posteriors are overplotted in solid 
lines, with a dotted line indicating zero probability.  Bottom panel: Various 
estimators of $\ln[n(z)]$, the interim prior, and the true $\ln[n(z)]$ as a 
continuous function and under a binned parametrization.}
		\label{fig:featured}
	\end{center}
\end{figure}

\subsubsection{Featured $n(z)$}
\label{sec:featured}

In this test, we choose the true $n(z)$ of Fig. \ref{fig:featured} with 
nontrivial structure.

\begin{figure}
	\begin{center}
		
\includegraphics[width=0.5\textwidth]{fig/fiducial/scatter.png}\\
		
\includegraphics[width=0.5\textwidth]{fig/fiducial/estimators.png}		
		\caption{[This plot doesn't really show what I want because the 
intrinsic scatter is too low!  I used $\sigma=0.03$ because that's what has 
been quoted as what LSST, etc. needs, but that's not what comes out of 
photo-$z$ estimation methods before they impose aggressive cuts\dots]  Top 
panel: The traditional MAP reduction of a photo-$z$ PDF against the true 
redshifts with a few rescaled photo-$z$ interim posteriors are overplotted in 
solid lines, with a dotted line indicating zero probability.  Bottom panel: 
Various estimators of $\ln[n(z)]$, the interim prior, and the true $\ln[n(z)]$ 
as a continuous function and under a binned parametrization.}
		\label{fig:featured}
	\end{center}
\end{figure}

This featured true $n(z)$ will henceforth be referred to as the fiducial $n(z)$

\subsection{Emulated Data Quality Effects}
\label{sec:likelihoods}

In the following test cases, we vary the properties of the mock photo-$z$ 
likelihoods in an effort to emulate known systematics in photo-$z$ estimation.  
These tests vary Step \ref{itm:scatter} of Sec. \ref{sec:mockdata}.

\subsubsection{Intrinsic Scatter}
\label{sec:intscat}

One major concern about photo-$z$s is the intrinsic scatter of point 
estimators, including those derived from photo-$z$ PDFs, that is observed to 
varying extents with every existing photo-$z$ algorithm and illustrated in Fig. 
\ref{fig:intscat}.  To emulate intrinsic scatter, we modify the fiducial case 
to simply broaden the single Gaussian component of the likelihood.  To enforce 
self-consistency, the mean is a random variable drawn from a Gaussian 
distribution with the newly increased variance.

\begin{figure}
	\begin{center}
		
\includegraphics[width=0.5\textwidth]{fig/high_scatter/scatter.png}\\
		
\includegraphics[width=0.5\textwidth]{fig/high_scatter/estimators.png}		
		\caption{As the intrinsic scatter increases, the discrepancy 
between estimators increases.  In particular, the stacked estimator and 
marginalized point estimators predict $\ln[n(z)]$ to be smoother than the 
truth, while the  [This would be a lot more compelling with more galaxies.  
Also, the weird edge effects in the top panel are real, because the point 
estimator is the MAP, not the center of the Gaussian likelihood, and there's no 
requirement that the mean of the likelihood be within the true redshift range.] 
 Top panel: The traditional MAP reduction of a photo-$z$ PDF against the true 
redshifts with a few rescaled photo-$z$ interim posteriors are overplotted in 
solid lines, with a dotted line indicating zero probability.  Bottom panel: 
Various estimators of $\ln[n(z)]$, the interim prior, and the true $\ln[n(z)]$ 
as a continuous function and under a binned parametrization.}
		\label{fig:intscat}
	\end{center}
\end{figure}

To emulate intrinsic scatter, we modify the fiducial case to simply broaden the 
single Gaussian component of the likelihood.  To enforce self-consistency, the 
mean is a random variable drawn from a Gaussian distribution with the newly 
increased variance.

\subsubsection{Template-like Catastrophic Outliers}
\label{sec:tempcatout}

In addition to intrinsic scatter, photo-$z$ methods employing template fitting 
tend to produce catastrophic outliers that are distributed to be broad in 
$z_{spec}$ and narrow in $z_{phot}$, as in Fig. \ref{fig:tempcatout}.  The 
systematic behind these catastrophic outliers may be described as an attractor 
in the space of $z_{phot}$; some galaxies at a range of $z_{spec}$ map onto a 
single $z_{phot}$ (with some scatter) if their true SED does not have 
sufficiently strong features (as is the case for blue galaxies), leading 
galaxies of that SED type at many $z_{spec}$ to have similar colors.

\begin{figure}
	\begin{center}
		
\includegraphics[width=0.5\textwidth]{fig/template_outliers/scatter.png}\\
		
\includegraphics[width=0.5\textwidth]{fig/template_outliers/estimators.png}	
	
		\caption{Top panel: The traditional MAP reduction of a 
photo-$z$ PDF against the true redshifts with a few rescaled photo-$z$ interim 
posteriors are overplotted in solid lines, with a dotted line indicating zero 
probability.  Bottom panel: Various estimators of $\ln[n(z)]$, the interim 
prior, and the true $\ln[n(z)]$ as a continuous function and under a binned 
parametrization.}
		\label{fig:tempcatout}
	\end{center}
\end{figure}

\subsubsection{Training-like Catastrophic Outliers}
\label{sec:traincatout}

Data driven photo-$z$ methods tend to suffer from a different form of 
catastrophic outliers that are distributed to be narrow in $z_{spec}$ and broad 
in $z_{phot}$, as in Fig. \ref{fig:traincatout}.  The systematic behind these 
catastrophic outliers may be described as an attractor in the space of 
$z_{spec}$; some galaxies near a particular $z_{spec}$ map to a range of 
$z_{phot}$ if the training set galaxies at that $z_{spec}$ have inconsistent 
$z_{phot}$, as might occur if their SED's features fall between photometric 
filters.

\begin{figure}
	\begin{center}
		
\includegraphics[width=0.5\textwidth]{fig/training_outliers/scatter.png}\\
		
\includegraphics[width=0.5\textwidth]{fig/training_outliers/estimators.png}	
		\caption{[This has proven to be the most challenging test case 
to implement, and there's clearly still a bug in the function that makes the 
likelihoods.]}
		\label{fig:traincatout}
	\end{center}
\end{figure}

\subsection{Emulated Interim Prior Effects}
\label{sec:priors}

The interim prior encapsulates the the relationship between observed photometry 
and redshift information upon which a photo-$z$ estimate is based.  Interim 
priors are in general not identical to the true $n(z)$ we wish to estimate; if 
they were, we would not need any data!  For template fitting photo-$z$ methods, 
the interim prior is usually an input chosen by the researcher.  However, for 
machine learning methods, the interim prior is some function of the training 
set data that in many cases may be influenced by random numbers and is rarely 
output with the redshift estimates.  Interim priors for template fitting 
methods tend to have incomplete coverage in the space of true photometry, 
because they are limited by the choice of the library of SEDs.  Interim priors 
for machine learning methods tend to have incomplete coverage in the space of 
redshifts, because there are fewer galaxies with spectroscopically confirmed 
redshifts at high redshifts than low redshifts.

Existing $n(z)$ estimation routines will always produce a biased estimator when 
the interim prior is not equal to the true $n(z)$.  We demonstrate here that 
regardless of the appropriateness of the interim prior as an approximation to 
the true $n(z)$, \chippr is not affected by the choice of the interim prior so 
long as it has nontrivial coverage in the space of redshift.  These tests 
modify Step \ref{itm:interim} of Sec. \ref{sec:mockdata}.

\subsubsection{Template-like Interim Prior}
\label{sec:tempintpr}

An interim prior based on a template library may be a sum of smooth functions 
representing $n(z)$ for each SED type in the library.  Template libraries do 
not include every possible galaxy SED, and the $n(z)$ used for each SED type 
may not be accurate.  The interim prior shown in Fig. \ref{fig:tempintpr} is an 
emulation of an interim prior corresponding to a template library of this type.

\begin{figure}
	\begin{center}
		
\includegraphics[width=0.5\textwidth]{fig/template_prior/scatter.png}\\
		
\includegraphics[width=0.5\textwidth]{fig/template_prior/estimators.png}	
		\caption{[The case of a multimodal interim prior was a very 
compelling test in the previous version but somehow isn't anymore.]}
		\label{fig:tempintpr}
	\end{center}
\end{figure}

\subsubsection{Training-like Interim Prior}
\label{sec:trainintpr}

An interim prior based on a training set may be biased toward low redshifts due 
to the dearth of distant galaxies with spectroscopic redshifts.  The interim 
prior shown in Fig. \ref{fig:trainintpr} is an emulation of an interim prior 
corresponding to a training set biased in this way.  We chose an interpolation 
of the interim prior used for the SDSS DR8 photo-$z$ PDFs.  

\begin{figure}
	\begin{center}
		
\includegraphics[width=0.5\textwidth]{fig/training_prior/scatter.png}\\
		
\includegraphics[width=0.5\textwidth]{fig/training_prior/estimators.png}	
		\caption{[The case of a low-$z$ favoring interim prior was a 
very compelling test in the previous version but somehow isn't anymore.]}
		\label{fig:trainintpr}
	\end{center}
\end{figure}

\section{Application to Realistic Mock Data}
\label{sec:application}

To show how the choice of $n(z)$ estimator propagates to cosmological 
constraints, we apply \chippr to a data from a realistic cosmological 
simulation (probably \texttt{Buzzard}) with photo-$z$ PDFs produced by a 
popular method (probably \texttt{BPZ}).

\subsection{Mock Data \& Metrics}
\label{sec:appintro}

As in Sec. \ref{sec:validintro}, the mock data takes the form of photo-$z$ 
interim posteriors, but that is where the similarity ends.  These photo-$z$ 
interim posteriors are derived from the photometry resulting from the 
\texttt{Buzzard} simulation by way of \texttt{BPZ}.  Because the simulation 
begins with setting true values of the cosmological parameters, we can 
propagate the different estimators of $n(z)$ through a forecasting code (which 
one?) to generate error ellipses on the cosmological parameters.

\subsubsection{Mock Data}
\label{sec:buzzard}

We summarize the details of the \texttt{Buzzard} simulation here.

\subsubsection{Cosmological Constraint Metric}
\label{sec:cosmo}

Because the realistic mock data of Sec. \ref{sec:buzzard} is associated with 
true values of the cosmological parameters, we may compare the quality of 
cosmological constraints under different estimators of $n(z)$, creating a 
figure like Fig. \ref{fig:scgsr}.  We perform this analysis using a forecasting 
code that takes $n(z)$ and produces projected error ellipses in the space of 
cosmological parameters.  

\begin{figure}
	\begin{center}
		Not my figure! (<NAME> via <NAME>)
		\includegraphics[width=0.5\textwidth]{cosmo_constraints.png}
		\caption{[I'd like to be able to make a plot like this with the 
different $n(z)$ estimators produced by \texttt{chippr} from the same set of 
photo-$z$ PDFs.  Then I can calculate some metrics of the accuracy and 
precision of the error distributions relative to the true values of the 
cosmological parameters that produced the mock data.]}
		\label{fig:scgsr}
	\end{center}
\end{figure}


\subsection{Results}
\label{sec:results}

\begin{figure}
	\begin{center}
		\caption{[moneyplot of error ellipses resulting from different 
$n(z)$ estimators]}
		\label{fig:money}
	\end{center}
\end{figure}

\section{Discussion}
\label{sec:discussion}

The results of Fig \ref{fig:money} have significant implications for the 
developing data analysis pipelines of next-generation telescope surveys.  
However, the method presented here has its own limitations, which are 
reiterated to discourage the community from applying this work inappropriately.

We intend to pursue a number of extensions of the work presented in this paper 
in future work.

\section{Conclusion}
\label{sec:conclusion}

We now summarize answers to the questions posed in the introduction:

\begin{itemize}
	\item Existing $n(z)$ estimation methods produce biased estimators that 
propagate to inaccuracies in characterizing the cosmological parameters.
	\item Photo-$z$ PDFs are probabilistic data products so must be handled 
in a mathematically consistent manner such as the probabilistic graphical model 
outlined in this paper.
	\item In addition to coming with its own error distribution, the $n(z)$ 
estimator produced by \chippr is quantifiably more accurate than established 
estimators.
	\item Propagation of the \chippr result leads to a quantifiable 
improvement in the constraints on cosmological parameters.
\end{itemize}

In conclusion, we discourage the community from continuing to use the stacked 
estimator and reductions of photo-$z$ PDFs to redshift point estimates in 
obtaining estimators of $n(z)$.  \chippr is freely available to the community 
for incorporation into evolving data analysis pipelines.  


\appendix{Catalog Production}
\label{app:catalog}

Sec. \ref{sec:mockdata} outlines one way to make mock data, but it's not 
necessarily a logically consistent way to do it, a problem identified in trying 
to use it for the test cases in Sec. \ref{sec:likelihoods}.  If the goal is to 
emulate the traditional $z_{spec}$ vs. $z_{phot}$ scatterplots, we need to 
start from probability distributions in that space.  This means there are two 
parameters $z_{spec}$ and $z_{phot}$ drawn from a two-dimensional distribution 
$p(z_{spec}, z_{phot} | \textul{\sigma}, \vec{\phi})$, where $\textul{\sigma}$ 
contains parameters concerning the relationship between $z_{spec}$ and 
$z_{phot}$, and $\vec{\phi}$ contains hyperparameters concerning $n(z)$.  We 
know from basic probability that
\begin{align}
\label{eq:trivial}
p(z_{spec}, z_{phot} | \textul{\sigma}, \vec{\phi}) =& p(z_{phot} | z_{spec}, 
\textul{\sigma}, \vec{\phi})\ p(z_{sped} | \textul{\sigma}, \vec{\phi}).
\end{align}
Because $z_{spec}$ depends only on $\vec{\phi}$ and $z_{phot}$ is independent 
of $\vec{\phi}$, this becomes
\begin{align}
\label{eq:simpler}
p(z_{spec}, z_{phot} | \textul{\sigma}, \vec{\phi}) =& p(z_{phot} | z_{spec}, 
\textul{\sigma})\ p(z_{spec} | \vec{\phi}).
\end{align}
We know how to make reasonable choices for both terms on the righthand side of 
Eq. \ref{eq:simpler}.

In the $z_{spec}$ vs. $z_{phot}$ scatterplot, $p(z_{phot} | z_{spec}, 
\textul{\sigma}, \vec{\phi})$ would be represented by vertical slices, but we 
are interested in the horizontal slices that represent $p(z_{spec} | z_{phot}, 
\textul{\sigma}, \vec{\phi})$.  Now we must distinguish the true 
hyperparameters $\vec{\phi}_{true}$ from the interim hyperparameters 
$\vec{\phi}_{int}$.  If the horizontal slices carry information about 
$\vec{\phi}$, how can we get a catalog of $p(z_{spec} | z_{phot}, 
\textul{\sigma}, \vec{\phi}_{int})$ that's consistent with $p(z_{spec} | 
z_{phot}, \textul{\sigma}, \vec{\phi}_{true})$?  Somehow that information needs 
to be propagated through in order to recover it.

We achieve this goal by dividing out the true $n(z)$ and multiplying by the 
interim $n(z)$ as in
\begin{align}
\label{eq:transformation}
p(z_{spec} | z_{phot}, \textul{\sigma}, \vec{\phi}_{int}) =& p(z_{spec} | 
z_{phot}, \textul{\sigma}, \vec{\phi}_{true})\ \frac{p(z_{spec} | 
\vec{\phi}_{int})}{p(z_{spec} | \vec{\phi}_{true})},
\end{align}
since $p(z_{spec} | z_{phot}, \textul{\sigma}, \vec{\phi})$ is separable and 
the term $p(z_{spec} | z_{phot}, \textul{\sigma})$ changes independently of the 
$p(z_{spec} | \vec{\phi})$ term.  The information about $\vec{\phi}_{true}$ is 
retained by the elements of the catalog: only $(z_{spec}, z_{phot})$ pairs 
drawn from the true distribution will be represented in the catalog!

It appears as if this is the only fully self-consistent way to ensure that a 
draw from the posterior distribution $p(z_{spec} | z_{phot}, \textul{\sigma}, 
\vec{\phi})$ does in fact follow that distribution.  In the cases with fiducial 
posteriors $p(z_{spec} | z_{phot}, \sigma, \vec{\phi})$, the method of Sec. 
\ref{sec:mockdata} was equivalent to this, but Sec. \ref{sec:likelihoods} 
requires a more sophisticated approach.

\begin{acknowledgements}
AIM thanks <NAME> for assistance with the \texttt{CosmoLike} code, 
<NAME> for insightful input on statistics, <NAME> for 
advice on debugging, and <NAME> for helpful comments provided in the 
preparation of this paper.
\end{acknowledgements}


\end{document}
{'fraction_non_alphanumeric': 0.11290764968869141, 'fraction_numerical': 0.0024685262898049864, 'mean_word_length': 4.652197504069452, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7625810', 'n_tokens_mistral': 22971, 'n_tokens_neox': 21701, 'n_words': 8935}