text:
\section{Introduction}
\label{sec:intro}

Surveys often ask respondents to report income in brackets or \textit{bins}, such as \$0-10,000, \$10,000-20,000,\ldots,\$200,000+. Even in surveys where respondents report exact incomes, incomes may be binned before publication, either to protect privacy or to summarize the income distribution compactly with the number of incomes in each bin. Table \ref{tab:nantucket} gives a binned summary of household incomes in Nantucket, the richest county in the US.

\begin{table}[tbh]
\centering
\begin{threeparttable}
\caption{Household incomes in Nantucket}
\label{tab:nantucket}
\begin{tabular}{cccc}
\hline
{\ul Min} & {\ul Max} & {\ul Households} & {\ul Cumulative distribution} \\
\$0         & \$10,000  & 165              & 5\%                 \\
\$10,000  & \$15,000  & 109              & 8\%                 \\
\$15,000  & \$20,000  & 67               & 9\%                 \\
\$20,000  & \$25,000  & 147              & 13\%                \\
\$25,000  & \$30,000  & 114              & 17\%                \\
\$30,000  & \$35,000  & 91               & 19\%                \\
\$35,000  & \$40,000  & 148              & 23\%                \\
\$40,000  & \$45,000  & 44               & 24\%                \\
\$45,000  & \$50,000  & 121              & 28\%                \\
\$50,000  & \$60,000  & 159              & 32\%                \\
\$60,000  & \$75,000  & 358              & 42\%                \\
\$75,000  & \$100,000 & 625              & 59\%                \\
\$100,000 & \$125,000 & 338              & 69\%                \\
\$125,000 & \$150,000 & 416              & 80\%                \\
\$150,000 & \$200,000 & 200              & 86\%                \\
\$200,000 &           & 521              & 100\%    \\ \hline
\end{tabular}
 \begin{tablenotes}
    \small
    \item
      \textit{Note}. Each bin's population is estimated from a 1-in-8 sample of households who took the American Community Survey in 2006-10. Incomes are in 2010 dollars.
    \end{tablenotes}

\end{threeparttable}
\end{table}

Binning presents challenges to investigators who want to estimate simple summary statistics such as the mean, median, or standard deviation, or inequality statistics such as the Gini coefficient, the Theil index, the coefficient of variation, or the mean log deviation. Researchers have implemented several methods for calculating estimates from binned incomes.

The simplest and most popular approach is to assign each case to the midpoint of its bin---using a robust pseudo-midpoint for the top bin, whose upper bound is typically undefined  (e.g., Table \ref{tab:nantucket}). The weakness of the midpoint approach is that it treats income as a discrete variable, but the method also has several strengths. The midpoint method is easy to implement and runs quickly. Midpoint estimates are also ``bin-consistent'' \citep{vonHippel2015} in the sense that midpoint estimates will get arbitrarily close to their estimands if the bins are sufficiently numerous and narrow.

Another approach is to fit the bin counts to a continuous parametric distribution. Popular distributions include 2-, 3-, and 4-parameter distributions from the generalized beta family, which includes the Pareto, lognormal, Weibull, Dagum, and other distributions \citep{McDonald2008}. One implementation fits up to 10 distributions and selects the one that fits best. An alternative is to use the AIC or BIC to calculated a weighted average of income statistics across several candidate distributions \citep{vonHippel2015}.

A strength of the parametric approach is that it treats income as continuous. A weakness is that even the best-fitting parametric distribution may not fit the bin counts particularly well. If the fit is poor, the parametric approach is not bin-consistent; that is, even with an infinite number of bins, each infinitesimally narrow, a parametric distribution may produce poor estimates if it is not a good fit to the underlying distribution of income.

A practical weakness of the parametric approach is that it is typically implemented using iterative methods which can be slow. The speed of a parametric fit may be acceptable if you fit a single distribution to a single binned dataset, but runtimes of hours are possible if you fit several distributions to thousands of binned datasets---such as every county or school district in the US. Other computational issues include nonconvergence and undefined estimates. These issues are rare but inevitable when you run thousands of binned datasets \citep{vonHippel2015}.

Neither approach---midpoint or parametric---is uniformly more accurate. With many bins, midpoint estimates are better because of bin-consistency, while with fewer than 8 bins, parametric estimates can be better because of their smoothness \citep{vonHippel2015}. Empirically, the parametric and midpoint approaches produce similarly accurate estimates from typical US income data with 15 to 25 bins. Both methods typically estimate the Gini within a few percentage points of its true value. This is accurate enough for many purposes, but can lead to errors when estimating small differences or changes such as the 5\% increase in the Gini of US family income that occurred between 1970 and 1980 \citep{vonHippel2015}.

A potential improvement is to fit binned incomes to a flexible nonparametric continuous density. Like parametric densities, a nonparametric density treats income as continuous. Like the midpoint method, a nonparametric density can be bin-consistent and fit the bin counts as closely as we like.

Unfortunately, past nonparametric approaches have been disappointing. A nonparametric approach using kernel density estimation had substantial bias under some circumstances \citep{Minoiu2012}. A nonparametric approach using a spline to model the log of the density \citep{Kooperberg1992} had even greater bias \citep{vonHippel2012}, though it was not clear whether the bias came from the method or its software implementation.

In this paper, we implement and test a nonparametric continuous method that outperforms its predecessors in both speed and accuracy. The method, which we call \textit{CDF interpolation}, simply connects points on the empirical cumulative distribution function (CDF). The method can connect the points using line segments or cubic splines. When cubic splines are used, the method is similar to ``histospline'' or ``histopolation'' methods which fit a spline to a histogram \citep{Wahba1976,Morandi1989}. But histosplines are limited to histograms which have bins of equal width \citep{Wang2015}. CDF interpolation is more general approach which can handle income data where the bins have unequal width and the top bins has no upper bound (e.g., Table \ref{tab:nantucket}).

We have implemented CDF interpolation in our R package \textit{binsmooth} \citep{Hunter2016}, which is available for download from the Comprehensive R Archive Network (CRAN).

Our results will show that statistics estimated with CDF interpolation are slightly more accurate than estimates obtained using midpoints or parametric distributions. In addition, CDF interpolation is much faster than parametric estimation, though not as fast as the midpoint method.

We also show that the differences between methods are dwarfed by the improvement we get if we constrain a method to match the grand mean of income, which the US Census Bureau often reports alongside the bin counts. If we constrain either the interpolated CDF or the midpoint method to match a known mean, we get dramatically better estimates of the Gini. Our \textit{binsmooth} package can constrain an interpolated CDF to match a known mean, and our new version (2.0) of the \textit{rpme} command for Stata can constrain the midpoint method to match a known mean as well.

In the rest of this paper, we define the midpoint, parametric, and interpolated CDF methods more precisely, then compare the accuracy of estimates in binned data summarizing household incomes within US counties. We also show how much estimates improve if we have the mean as well as the bin counts.

\section{Methods}
\label{sec:meth}

\subsection{Binned data}
A binned data set, such as Table \ref{tab:nantucket}, consists of counts $n_1, n_2, \ldots, n_B$ specifying the number of cases in each of $B$ bins. The total number of cases is $T = n_1 + n_2 + \cdots + n_B $. Each bin $b$ is defined by an interval $[l_b,u_b), b=1, \ldots, B$, where $l_b$ and $u_b$ are the lower and upper bound of income for that bin. The bottom bin often starts at zero ($l_1=0$), and the top bin may have no upper bound ($u_B=\infty$).

\subsection{A midpoint method}

The oldest and simplest way to analyze binned data is the midpoint method, which within each bin $b$ assigns incomes to the bin midpoint $m_b=(l_b+u_b)/2$ (\citealt{}{Heitjan1989} for a review). Then statistics such as the Gini can be calculated by applying sample formulas to the midpoints $m_b$ weighted by the counts $n_b$.

When the top bin has no upper bound, we must define a pseudo-midpoint for it. The traditional choice is ${\mu}_B = l_B \alpha/(\alpha-1)$, which would define the arithmetic mean of the top bin if top-bin incomes followed a Pareto distribution with shape parameter $\alpha > 1$ \citep{Henson1967}. The problem with this choice is that the arithmetic mean of a Pareto distribution is undefined if $\alpha \leq 1$ and grows arbitrarily large as $\alpha$ approaches 1 from above.

A more robust choice is the harmonic mean $h_B = l_B (1+1/\alpha)$, which is defined for all $\alpha>0$ \citep{vonHippel2015}. We use the harmonic mean in this article. We estimate $\alpha$ by fitting a Pareto distribution to the top two bins and calculating the maximum likelihood estimate from the following formula \citep{Quandt1966}:
\begin{equation}
    \hat\alpha = \frac{ln((n_{B-1}+n_B)/n_B)}
                        {ln(l_B/l_{B-1})}
    \label{eqn:alpha}
\end{equation}

If the survey provides the grand mean $\mu$, then we need not assume that top-bin incomes follow a Pareto distribution. Instead, we can calculate
\begin{equation}
    \hat\mu_B = \frac{1}{n} (T \mu - \sum_{b=1}^{B-1} n_B m_B )
    \label{eqn:top_bin_mean}
\end{equation}

\noindent which would be the mean of the top bin if the means of the lower bins were the midpoints $m_b$. Then $\hat\mu_B$ can serve as a pseudo-midpoint for the top bin.

When estimated in this way, it occasionally happens (e.g., in 4 percent of US counties) that the top bin's mean $\hat\mu_B$ is slightly less than its lower bound $l_B$. This is infelicitous, but $\hat\mu_B$ can still be used, and the resulting Gini estimates are not necessarily bad.\footnote{It might be a little better to set $\hat\mu_B=l_B$ and move the other midpoints slightly the the left, but this would affect only 4 percent of counties, and those only slightly since the affected counties typically have few cases in the top bin.}

The midpoint methods described in this section are implemented by the \textit{rpme} command for Stata \citep{Duan2016}, where  \textit{rpme} stands for ``robust Pareto midpoint estimator'' \citep{vonHippel2016}. Except for mean-matching, the approach is also implemented in the \textit{binequality} package for R \citep{Scarpino2017}.

\subsection{Fitting parametric distributions}

The weakness of the midpoint method is that it treats income as discrete. An alternative is to model income as a continuous variable $X$ that fits some parametric CDF $F(X|\theta)$. Here $\theta$ is a vector of parameters, which can be estimated by iteratively maximizing the log likelihood:
\begin{equation}
    \ell(\theta|X) = ln \prod_{b=1}^{B} (P(l_b<X<u_b))^{n_b}
        \label{eqn:log_likelihood}
\end{equation}

\noindent where $P(l_b<X<u_b)=F(u_b)-F(l_b)$ is the probability, according to the fitted distribution, that an income is in the bin $[l_b,u_b)$.\footnote{This formula includes the top bin $B$, where $u_B=\infty$ and $F(u_B)=1$. Some older articles (e.g., \citep{McDonald1979,McDonald1995,Bandourian2002}) add the following constant to the log likelihood: $ln(T!)-\sum_{b=1}^{B}{ln(n_b!)}$. This is not wrong, but it is unnecessary since adding a constant does not change the parameter values at which the log likelihood is maximized \citep{edwards1972}.}

While any parametric distribution can be considered, in practice it is hard to fit a distribution unless the number of parameters is small compared to the number of well-populated bins. Most investigators favor 2-, 3-, and 4-parameter distributions from the generalized beta family \citep{McDonald2008}, which includes the following 10 distributions: the log normal, the log logistic, the Pareto (type 2), the gamma and generalized gamma, the beta 2 and generalized beta (type 2), the Dagum, the Singh-Maddala, and the Weibull.

\textit{A priori} it is hard to know which distribution, if any, will fit well. The fit of a distribution can be tested by the following goodness-of-fit likelihood ratio statistic \citep{vonHippel2015}:

\begin{equation}
G^2=-2 (\hat{\ell} - \sum_{b=1}^{B}{n_b ln(n_b/T)}
\end{equation}

\noindent where $\hat{\ell}$ is the maximized log likelihood. Under the null hypothesis that the fitted distribution is the true distribution of income, $G^2$ would follow a chi-square distribution with $B^*-k$ degrees of freedom, where $k$ is the number of parameters, $B^*=min(B_{>0},B-1)$, and $B_{>0}$ is the number of bins with nonzero counts. We reject the fit of a distribution if $G^2$ has $p<.05$ in the null distribution.

In empirical data it is common to reject every distribution in the generalized beta family \citep{Bandourian2002,vonHippel2015}. In addition, some distributions may fail to converge, or may converge on parameter values that imply that the mean or variance of the distribution is undefined \citep{vonHippel2015}.

A solution is to fit all 10 distributions, screen out any with undefined moments and, among the distributions remaining, select the one that fits best according to the Akaike or Bayes information criterion (AIC or BIC):
\begin{align*}
AIC &= 2k - 2 \hat{\ell} \\
BIC &= ln(T) k  - 2 \hat{\ell}
\end{align*}

Or, instead of selecting a single best-fitting distribution, one can average estimates across several candidate distributions weighted proportionately to a function of the AIC or BIC (specifically $exp(-AIC/2)$ or $exp(-BIC/2)$) an approach known as \textit{model averaging}  \citep{Burnham2004}. In general, model averaging yields better estimates than model selection, but when modeling binned incomes the advantage of model averaging is negligible \citep{vonHippel2015}.

Although it sounds broad-minded to fit 10 different distributions, there is limited diversity in the generalized beta family. All the distributions in the generalized beta family are unimodal and skewed to the right. Some distributions are quite similar (e.g., Dagum and generalized beta), and others rarely fit well (e.g., log normal, log logistic, Pareto). So in practice the range of viable and contrasting distributions in the generalized beta family is small; you can fit just 3 well-chosen distributions (e.g., Dagum, gamma, and generalized gamma) and get estimates almost as good as those obtained from fitting all 10 \citep{vonHippel2015}.

We use the fitted distributios to estimate income statistics such as the mean, variance, Gini, or Theil. Sometimes the income statistic is a simple function of the distributional parameters $\theta$, but other times the function is unknown, or hard to calculate. As a general solution, it is easier to calculate income statistics by applying numeric integration to appropriate functions of the fitted distribution \citep{McDonald2008}.

When the grand mean is available, the distributional parameters could in theory be constrained to match the grand mean as well as approximate the bin counts. This would be difficult, though, since for distributions in the generalized beta family the mean is a complicated nonlinear function of the parameters. We have not attempted to constrain our parametric distributions to match a known mean.

The parametric approaches described in this section are implemented in the \textit{binequality} package for R \citep{Scarpino2017} and the \textit{mgbe} command for Stata \citep{Duan2016}. Here \textit{mgbe} stands for ``multi-model generalized beta estimator'' \citep{vonHippel2015}.

\subsection{Interpolated CDFs}

Since parametric distributions may fit poorly, an alternative is to define a flexible nonparametric density that fits the bin counts exactly. Our nonparametric approach is called CDF interpolation.

To understand CDF interpolation, consider that binned data define $B$ discrete points on the empirical cumulative distribution function (CDF). The empirical CDF for Nantucket is given in the last column of Table \ref{tab:nantucket}, which shows that 5\% of Nantucket households make less than \$10,000, 8\% make less than \$15,000, and so on. Formally, at an income of 0 the empirical CDF is $\hat{F}(0)=0$, and at each bin's upper bound $u_b$ the empirical CDF is the observed fraction of incomes that are less than $u_b$---i.e., $\hat{F}(u_b) = (n_1 + n_2 + \cdots + n_b)/T$.\footnote{If the binned incomes come from a sample (as they do in Table \ref{tab:nantucket}), then the estimate $\hat{F}(u_b)$ may differ from the true CDF ${F}(u_b)$ because of sampling error.}

Now to estimate a continuous CDF $\hat{F}(x), x>0$, we just connect the dots. That is, we define a continuous nondecreasing function that interpolates between the $B$ discrete points of the empirical CDF.

Then the estimated probability density function (PDF) is just the derivative of the interpolated CDF $\hat{F}(x)$. Note that the estimated PDF ``preserves areas'' \citep{Morandi1989}---i.e., preserves bin counts. That is, within in each bin ($\hat{P}(l_b<x<u_b)=\hat{F}(u_b)-\hat{F}(l_b)$), the total of the estimated density is equal to the observed fraction of incomes that are in that bin ($n_b/T$).

The shape of the PDF depends on the function that interpolates the CDF.

\begin{itemize}
 \item If the CDF is interpolated by line segments, then the CDF is polygonal, and the PDF is a step function that is discontinuous at the bin boundaries.
 \item If the CDF is interpolated more smoothly, by a continuously differentiable monotone cubic spline, then the PDF is piecewise quadratic --- i.e., continuous at the bin boundaries and quadratic between them.
 \end{itemize}

There remains a question of how to shape the CDF in the top bin, which typically has no upper bound. In our implementation, the CDF of the top bin can be rectangular, exponential, or Pareto.\footnote{In our implementation, the exponential and Pareto tails are approximated by a sequence of rectangles of decreasing heights.} Each of these distributions has one parameter, which is estimated as follows.

\begin{itemize}
\item If the grand mean of income is known, then the parameter shaping the top bin is constrained so that the mean of the fitted distribution matches the grand mean. It occasionally happens (e.g., in 4 percent of US counties) that we cannot reproduce the known mean this way, because the known mean is already less than the mean of the lower $B-1$ bins without the tail. In that case, we make an \textit{ad hoc} adjustment by shrinking the bin boundaries toward the origin -- that is, by replacing $(l_b,u_b)$ with $(s l_b,s u_b)$, where the shrinkage factor $s<1$ is chosen so that a small tail can be added to reproduce the grand mean. The shrinkage factor is rarely less than .995.

\item If the mean income is not known, we substitute an \textit{ad hoc} estimate. We obtain that estimate by temporarily setting the upper bound of the top bin to $u_B=2l_B$ and calculating the mean of a step PDF fit to all $B$ bins. Then we unbound the top bin and proceed as though the mean were known.
\end{itemize}

Income statistics are estimated by applying numerical integration to functions of the fitted PDF or CDF.

The methods in this section are implemented by our \textit{binsmooth} package for R \citep{Hunter2016}. Within the \textit{binsmooth} package, the \textit{stepbins} function implements a step-function PDF (and polygonal CDF), while the \textit{splinebins} function implements a cubic spline CDF (and piecewise quadratic PDF).

\subsection{Recursive subdivision}

Another way to obtain a smooth PDF estimate that preserves bin areas is to subdivide the bins into smaller bins, and then adjust the heights of the subdivided bins to shorten the jumps at the bin boundaries. This method, \textit{recursive subdivision}, is implemented by the \textit{recbin} command in the \textit{binsmooth} package for R \citep{Hunter2016}. Recursive subdivision is slower and more computationally intensive than CDF interpolation and the resulting estimates are practically identical. We present the details of recursive subdivision in the Appendix.

\section{Data and Results}

Between 2006 and 2010, the American Community Survey (ACS) took a 1-in-8 sample of US households \citepalias{uscacs}. Household incomes were inflated to 2010 dollars and summarized in binned income tables for  each of the 3,221 US counties. The published bin counts are estimates of the population counts. We can approximate the sample counts by dividing the population counts by 8. Dividing counts by a constant makes no difference to any of our statistics, except for the BIC and $G^2$ statistics that are used when fitting parametric distributions.

The Census also published means and Gini coefficients for each county \citep{Bee2012}. These statistics were estimated from exact incomes before binning, and so are more accurate than any estimate that could be calculated from the binned data. They are sample estimates which may differ from population values, but they remain a useful standard of comparison for our binned-data estimates.

\subsection{Results for Nantucket}

Table \ref{tab:nantucket} summarized the binned incomes for Nantucket County. Figure~\ref{fig:comparepdfs} fits several distributions to the Nantucket data.

    \begin{figure}[hbt]
        \centering
         \caption{Different PDFs for the Nantucket data.}
        \includegraphics[width=6.5in]{nantucket.eps}
        \label{fig:comparepdfs}
        \floatfoot{The Dagum distribution is drawn in purple, the quadratic spline in blue, and the step function in black. Gray spikes illustrate the midpoint method.}
    \end{figure}

The midpoint method is illustrated by gray spikes at the bin midpoints; the spike heights are proportional to the bin counts. The black step function is the PDF implied by a linear interpolation of the CDF, and the blue curve is the piecewise quadratic PDF implied by a cubic spline of the CDF. Both CDF interpolations fit the bin counts perfectly; in fact, their jagged appearance suggests they may \textit{overfit} the data---a concern that we will revisit in the Conclusion. The step PDF looks slightly less volatile than the piecewise quadratic PDF, suggesting that the step PDF may be less overfit.

The purple curve is the Dagum distribution. The Dagum fits Nantucket better than other distributions from the generalized beta family, but it does not fit well. It fails the $G^2$ goodness-of-fit test, and visually it fits the bin counts poorly. For example, between \$70,000 and \$150,000 the Dagum curve suggests there should be substantially fewer households than there are, and above \$150,000 it suggests that there should be more.

Except for the Dagum distribution, all the methods in Figure~\ref{fig:comparepdfs} are calibrated to reproduce the grand mean.

Table~\ref{tab:nantucket_estimates} summarizes the Nantucket estimates. The true mean is \$137,000 and the true Gini is .547. All the methods underestimate these quantities. When fit without knowledge of the true mean, every method underestimates the mean by 12 to 20 percent and the Gini by 15 to 21 percent, with the simple midpoint method coming closer than its more sophisticated competitors.

\begin{table}[hbt]
\centering
\caption{Estimates for Nantucket}
\label{tab:nantucket_estimates}
\begin{tabular}{llllll}
\hline
                   & {\ul Estimate}      & {\ul Mean}      & {\ul Gini} \\
                   & True                     & \$137,811 & .547 &  &  \\
Without the true mean  & Midpoint (RPME)          & \$121,506 & .464 &  &  \\
                   & Parametric (MGBE)        & \$112,960 & .453 &  &  \\
                   & CDF interpolation (linear) & \$110,419 & .438 &  &  \\
                   & CDF interpolation (cubic spline)  & \$110,419 & .433 &  &  \\
Matching the true mean & Midpoint                 &           & .510 &  &  \\
                   & CDF interpolation (linear) &           & .537 &  &  \\
                   & CDF interpolation (cubic spline)  &           & .525 &  & \\
                   \hline
\end{tabular}
\end{table}

When given the true mean, the midpoint and CDF interpolation methods do much better. They still underestimate the Gini, but only by 2 to 7 percent. The closest estimate is obtained by linear interpolation of the CDF. A smoother cubic spline interpolation does a little worse, but still better than the midpoint method.

Although the estimates for Nantucket are less accurate overall than the estimates for most other counties, the relative performance of different methods in Nantucket is similar to what we will see elsewhere.

\subsection{Results for all US counties}
\label{sec:bakeoff}

Figure~\ref{fig:scatterplots} evaluates all county Gini estimates graphically by plotting the estimated Gini $\hat{\theta}_j$ of each county against the published Gini ${\theta}_j$. In the bottom row, where the methods are constrained to match the published mean, the Gini estimates are close to a diagonal reference line ($\hat{\theta}_j=\theta$) indicating nearly perfect estimation. In the top row, where the methods do not match the mean, the estimates are more scattered, indicating lower accuracy.

    \begin{figure}[hbt]
        \centering
        \includegraphics[width=6in]{Compare_accuracy.eps}
        \caption{Accuracy of Gini coefficients estimated by different methods, with and without mean-matching.}
        \label{fig:scatterplots}
    \end{figure}

We can summarize the accuracy of Gini estimates in several ways. For a single county $j$, the percent estimation error is $e_j = 100\times(\hat{\theta}_j - \theta_j)/\theta_j$. Then across all counties, the \textit{percent relative bias} is the mean of $e_j$, the \textit{percent root mean squared error} (RMSE) is the square root of the mean of $e_j^2$, and the \textit{reliability} is the squared correlation between $\theta_j$ and $\hat{\theta}_j$.

Table~\ref{tab:errorres} summarizes our findings. If the estimators ignore the published county means, then estimated Ginis have biases between 0\% and -3\%, RMSEs between 3\% and 4\%, and reliabilities between 82\% and 88\%. The interpolated CDF estimates have the best bias, the best RMSE, and the second best reliability, and they are just as good with linear interpolation as with cubic spline interpolation. The parametric estimates have the best reliability, but the worst bias, the worst RMSE, and by far the worst runtime, at 4.5 hours.

\begin{table}[tbh]
\begin{threeparttable}[tbh]
    \caption{Speed and accuracy of Gini estimates for 3,221 US counties.
    \label{tab:errorres}}
\centering
\begin{tabular}{llcccl}
\hline
                  &     {\ul Estimator}              & {\ul Bias} & {\ul RMSE} & {\ul Reliability} & {\ul Runtime} \\
Without true mean & Midpoint                        & -2\%       & 4\%            & 82\%              & 4 sec         \\
                  & Parametric                      & -3\%       & 4\%            & 89\%              & 4.5 hr        \\
                  & CDF interpolation (linear)     & 0\%        & 3\%            & 88\%              & 40 sec        \\
                  & CDF interpolation (cubic spline) & -1\%       & 3\%            & 88\%              & 2 min         \\
Matching true mean    & Midpoint                               & -1\%       & 2\%            & 98\%              & 7 sec         \\
                  & CDF interpolation (linear)     & 0\%        & 1\%            & 99\%              & 36 sec        \\
                  & CDF interpolation (cubic spline) & -1\%       & 1\%            & 99\%              & 2 min      \\
\hline

\end{tabular}
    \begin{tablenotes}
    \small
    \item
\textit{Note}. RMSE=root mean squared error.
Runtimes on an Intel i7 Core processor with a speed of 3.6 GHz.
    \end{tablenotes}
    \end{threeparttable}
\end{table}

When the methods are constrained to match the published county means, the estimates improve dramatically. The bias shrinks to 0-1\%, the RMSE shrinks to 1-2\%, and the reliability grows to 98-99\%. The midpoint estimates are excellent, and the interpolated CDF estimates are even better, and just as good with linear interpolation as with cubic spline interpolation.

The differences among the methods are much smaller than the improvement that comes from constraining any method to match the mean. Of course, this observation is only helpful when the mean is known.

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

CDF interpolation produces estimates that are at least a little better than midpoint or parametric estimates, whether the true mean is known or not. And CDF interpolation runs much faster than parametric estimation, thought not as fast as midpoint estimation.

We initially suspected that cubic spline interpolation would improve on simple linear interpolation, but empirically this turns out to be false. In estimating county Ginis, linear CDF interpolation was at least as accurate as cubic spline interpolation.

The accuracy of linear CDF interpolation is remarkable, since it implies a step function for the PDF. Step PDFs look clearly unrealistic, especially in the top and bottom bins where the step function is flat while the true distribution likely has an upward or downward slope \citep{cloutier1995}. Our step PDF permits a downward Pareto or exponential slope in the top bin, and our cubic spline CDF can fit an upward slope to the bottom bin. But neither of these refinements does much to improve the accuracy of Gini estimates.

The differences in accuracy among the methods are small, and they are dwarfed by the improvement in accuracy that comes from knowing the grand mean. By constraining binned-data methods to match a known mean, we can typically get county Gini estimates that are within 1-2\% of the estimates we would get if the data were not binned. Our \textit{binsmooth} package for R can constrain interpolated CDFs to match a known mean, and our \textit{rpme} command for Stata can constrain the top-bin midpoint to match a known mean as well. We have not constrained our parametric distributions to match a known mean; we believe it would be difficult to do so.

While the mean-constrained estimates are very accurate, there may be room for improvement when the mean is unknown. Perhaps the most promising idea is additional smoothing. As we noticed in Figure~{fig:comparepdfs}, interpolated CDFs can be a bit jagged and may ``overfit'' the sample in the sense that they find nooks and crannies that might not appear in another sample or in the population. Likewise interpolated CDFs may be overfit to a specific set of bin boundaries. If the fitted CDF were a little smoother and did not quite preserve the counts of the least populous bins, it might fit the population and other samples (perhaps with different bin boundaries) a little better.


text:
\section{Introduction}
X-ray emission from clusters is mainly due to bremsstrahlung emission
\citep{Felten66,Mitchell76} from the hot intracluster medium
(ICM). The X-ray flux is proportional to the square of density, with
some temperature dependence, and so is a sensitive tracer of
variations in the thermodynamic properties of the ICM.

Although relaxed clusters are largely hot atmospheres in hydrostatic
pressure equilibrium, density and temperature variations are important
probes of astrophysical processes within the cluster and of cluster-wide
perturbations such as mergers. AGN feedback in clusters
\citep{Fabian12} injects bubbles of radio plasma into the ICM,
displacing the X-ray emitting gas and creating cavities in X-ray
images \citep[e.g.][]{BohringerPer93,FabianPer00,McNamara00}. In
addition, AGN are observed to shock their surroundings
\citep[e.g.][]{Forman07,Randall15} and likely generate sound waves in
the ICM \citep{FabianPer06,SandersSound08,Blanton11}. These processes
are seen in simulations of AGN feedback
\citep[e.g.][]{Ruszkowski04,Sijacki06}.

Cold fronts --- discontinuities in temperature and density --- are
extremely common in clusters \citep{MarkevitchCFShock07}. These are
seen in both merging \citep[e.g.][]{Vikhlinin01} and relaxed
\citep[e.g.][]{PaternoMahler13} clusters. In the relaxed cases the fronts are
believed to be caused by gas sloshing in the potential well due to the
passage of a subcluster \citep{Ascasibar06}. Such sloshing may remain
for several Gyr and can give rise to several edges within a single
cluster. Merging subclusters additionally generate shocks, also seen
as surface brightness edges \citep[e.g.][]{Markevitch02,Russell10}.

Mergers and AGN feedback can inject turbulence in the ICM
\citep{NormanBryan99}. Turbulence within the ICM should be associated
with gas density and therefore surface brightness variations which can
be quantified
\citep{Churazov12,SandersAWM712,ZhuravlevaFluct14,Walker15}.

The radial variation in density in clusters, particularly in relaxed
objects, gives rise to steeply peaked surface brightness
profiles. Therefore it is often difficult to see variations in surface
brightness profile on these mountainous X-ray peaks.  Various
techniques have been developed to enhance the structure seen in images
of galaxy clusters, allowing variations to be observed. These include
dividing or subtracting a symmetric model, such as the $\beta$ model
\citep[e.g.][]{Arnaud01}, or the average at each radius
\citep[e.g.][]{Churazov99}. More complex models can also be used, such
as fits to surface brightness contours with ellipses
\citep{SandersAWM712}.

Another common technique is to use unsharp masking of cluster images
to highlight smaller-scale substructure
\citep[e.g.][]{FabianPer03}. The typical method is to smooth the image
first by a Gaussian with a large width and subtract (or divide) this
from the same image smoothed by a Gaussian with a smaller width. The
process results in an image which suppresses both large and small
scale structure. A further method to suppress
unwanted cluster signal
is a Fourier bandpass filter \citep{SandersSound08}. However, a risk
with Fourier techniques is that spurious circular ringing artifacts
from edges and point sources can be introduced if frequency cut-offs
in the applied filter are too sharp.

Here we apply a gradient measuring filter to X-ray images of galaxy
clusters. Finding edges and measuring gradients in surface brightness
are useful, because all the astrophysical processes we have previously
mentioned introduce density variations, and therefore surface
brightness gradients, into X-ray images. For example, shocks and cold
fronts produce edges (i.e. very steep gradients), while sound waves
should produce alternating flat and steep gradients. As we are not
interested in the total X-ray emission in a region, using the gradient
removes much of the X-ray peak. Gradient filtering has previously been
used in examining simulations of galaxy cluster to look for edges
associated with cold fronts and sloshing \citep{Roediger13}.

The use of the Gaussian Gradient Magnitude (GGM) filter was introduced
to X-ray analysis in our study of deep \emph{Chandra} observations of
the Centaurus cluster \citep{SandersCent16}.  The GGM filter
calculates the gradient of an image assuming Gaussian derivatives
(with a width $\sigma$). In comparison, the Sobel operator,
a type of
gradient filter, convolves two $3\times3$ matrices with the image and
can be used to compute the magnitude of the gradient in an image on a
pixel-by-pixel basis.  The advantage of the GGM filter over the Sobel
filter is that $\sigma$ can be adjusted to measure the gradient over
more or fewer pixels depending on the data quality and the scale of
the features of interest. In X-ray cluster images, a large $\sigma$
would be used for regions in the outskirts where there are few or no
X-ray counts per pixel and a small value in the centre where there are
many counts.

In this work we apply the GGM filter to other high quality data sets
from the \emph{Chandra} archive to demonstrate that the technique is a
powerful method for the identification of physical processes taking
place in the ICM. We examine two relaxed clusters with short central
cooling times and active AGN feedback and one disturbed system
undergoing a merger. In the Perseus cluster, A\,426, there are
multiple X-ray cavities \citep{Bohringer04,FabianPer00}, a weak shock
\citep{FabianPer03}, ripples \citep{FabianPer06} and uplifted high
metallicity material \citep{SandersNonTherm05}, which may be sound
waves from AGN feedback. M\,87 contains a bright jet, multiple
bubble-like cavities \citep{Young02} and a weak shock
\citep{Forman05,Forman07}. There are cool arms of metal rich material
being dragged out by the radio bubbles
\citep{Young02,Simionescu07,Million10}. A\,3667 is a system undergoing a
merger \citep{Knopp96} and hosts a sharp surface brightness
discontinuity, a cold front, indicating material is moving through the
ambient gas with a Mach number of $\sim 1$
\citep{Vikhlinin01}. Optically, the cluster has two distinct sets of
galaxies \citep{Sodre92} and a radio relic \citep{Rottgering97}.

\section{Data preparation}
\begin{table}
  \caption{\emph{Chandra} data sets examined for the Perseus cluster.
    For each observation we list the observation identifier, observation
    starting date, exposure (unfiltered and after filtering for flares in ks)
    and the ACIS CCDs from which the data were examined, with the ACIS mode
    (S or I). $^{*}$ marks the reference observations that others
    were reprojected to.}
  \centering
  \begin{tabular}{rrrrr}
    OBSID & Date & Exposure & Filtered & CCDs \\ \hline
      502 & 1999-09-20 &   5.1 &   2.4 &  I: 0,1,2,3,6,7 \\
      503 & 1999-11-28 &   9.0 &   8.8 &  S: 2,3,6,7,8 \\
     1513 & 2000-01-29 &  24.9 &  10.5 &  S: 2,3,6,7,8 \\
     3209 & 2002-08-08 &  95.8 &  94.0 &  S: 1,3,6,7 \\
     4289 & 2002-08-10 &  95.4 &  93.0 &  S: 1,3,6,7 \\
     6139 & 2004-10-04 &  56.4 &  53.2 &  S: 2,3,5,6,7 \\
     4946 & 2004-10-06 &  23.7 &  23.3 &  S: 2,3,5,6,7 \\
     4948 & 2004-10-09 & 118.6 & 111.3 &  S: 2,3,5,6,7 \\
     4947 & 2004-10-11 &  29.8 &  29.4 &  S: 2,3,5,6,7 \\
     4949 & 2004-10-12 &  29.4 &  29.2 &  S: 2,3,5,6,7 \\
     4950 & 2004-10-12 &  96.9 &  75.6 &  S: 2,3,5,6,7 \\
$^{*}$4952& 2004-10-14 & 164.2 & 147.3 &  S: 2,3,5,6,7 \\
     4951 & 2004-10-17 &  96.1 &  93.8 &  S: 2,3,5,6,7 \\
     4953 & 2004-10-18 &  30.1 &  29.7 &  S: 2,3,5,6,7 \\
     6145 & 2004-10-19 &  85.0 &  84.6 &  S: 2,3,5,6,7 \\
     6146 & 2004-10-20 &  47.1 &  42.4 &  S: 2,3,5,6,7 \\
    11716 & 2009-10-10 &  39.6 &  38.1 &  I: 0,1,2,3,6 \\
    12025 & 2009-11-25 &  17.9 &  11.7 &  I: 0,1,2,3,6,7 \\
    12033 & 2009-11-27 &  18.9 &  12.0 &  I: 0,1,2,3,6,7 \\
    11713 & 2009-11-29 & 112.2 &  75.9 &  I: 0,1,2,3,6,7 \\
    12036 & 2009-12-02 &  47.9 &  34.8 &  I: 0,1,2,3,6,7 \\
    11715 & 2009-12-02 &  73.4 &  68.6 &  I: 0,1,2,3,6 \\
    12037 & 2009-12-05 &  84.6 &  79.1 &  I: 0,1,2,3,6 \\
    11714 & 2009-12-07 &  92.0 &  80.1 &  I: 0,1,2,3,6 \\ \hline
    Total &            &       &1328.7 &              \\
  \end{tabular}
  \label{tab:per_obs}
\end{table}

\begin{table}
  \caption{\emph{Chandra} data sets examined for M\,87.
    The columns are the same as for Table \ref{tab:per_obs}.}
  \centering
  \begin{tabular}{rrrrr}
    OBSID & Date & Exposure & Filtered & CCDs \\ \hline
      352 & 2000-07-29 &  37.7 &  33.5 & S: 2,3,6,7 \\
     3717 & 2002-07-05 &  20.6 &  10.5 & S: 2,3,6,7 \\
$^{*}$2707& 2002-07-06 &  98.7 &  88.2 & S: 2,3,6,7 \\
     6186 & 2005-01-31 &  51.6 &  42.8 & I: 0,1,2,3 \\
     5826 & 2005-03-03 & 126.8 & 120.0 & I: 0,1,2,3 \\
     5827 & 2005-05-05 & 156.2 & 147.7 & I: 0,1,2,3 \\
     7212 & 2005-11-14 &  65.2 &  61.2 & I: 0,1,2,3 \\
     7210 & 2005-11-16 &  30.7 &  27.5 & I: 0,1,2,3 \\
     7211 & 2005-11-16 &  16.6 &  15.7 & I: 0,1,2,3 \\
     5828 & 2005-11-17 &  33.0 &  31.4 & I: 0,1,2,3 \\
    15180 & 2013-08-01 & 138.8 & 137.2 & I: 0,1,2,3 \\
    15178 & 2014-02-17 &  46.5 &  46.1 & I: 0,1,2,3 \\
    16585 & 2014-02-19 &  45.0 &  44.0 & I: 0,1,2,3 \\
    16586 & 2014-02-20 &  49.2 &  48.4 & I: 0,1,2,3 \\
    16587 & 2014-02-22 &  37.3 &  37.0 & I: 0,1,2,3 \\
    15179 & 2014-02-24 &  41.4 &  40.2 & I: 0,1,2,3 \\
    16590 & 2014-02-27 &  37.6 &  37.0 & I: 0,1,2,3 \\
    16591 & 2014-02-27 &  23.5 &  22.9 & I: 0,1,2,3 \\
    16592 & 2014-03-01 &  35.6 &  35.0 & I: 0,1,2,3 \\
    16593 & 2014-03-02 &  37.6 &  36.8 & I: 0,1,2,3 \\ \hline
    Total &            &       &1063.2 &            \\
  \end{tabular}
  \label{tab:m87_obs}
\end{table}

\begin{table}
  \caption{\emph{Chandra} data sets examined for A\,3667.
    The columns are the same as for Table \ref{tab:per_obs}.}
  \centering
  \begin{tabular}{rrrrr}
    OBSID & Date & Exposure & Filtered &  CCDs          \\ \hline
      513 & 1999-09-22 &  44.8 &  40.1 & I: 0,1,2,3,5,6 \\
      889 & 2000-09-09 &  50.3 &  49.7 & I: 0,1,2,3,6,7 \\
$^{*}$5751& 2005-06-07 & 128.9 & 125.0 & I: 0,1,2,3,6   \\
     6292 & 2005-06-10 &  46.7 &  45.7 & I: 0,1,2,3,6   \\
     5752 & 2005-06-12 &  60.4 &  59.4 & I: 0,1,2,3,6   \\
     6295 & 2005-06-15 &  49.5 &  48.9 & I: 0,1,2,3,6   \\
     5753 & 2005-06-17 & 103.6 &  74.0 & I: 0,1,2,3,6   \\
     6296 & 2005-06-19 &  49.4 &  48.6 & I: 0,1,2,3,6   \\
     7686 & 2007-06-23 &   5.0 &   5.0 & I: 0,1,2,3     \\ \hline
    Total &            &       & 496.5 &                \\
  \end{tabular}
  \label{tab:a3667_obs}
\end{table}

\begin{figure*}
  \includegraphics[width=\textwidth]{fig01_1_small.pdf}
  \includegraphics[width=\textwidth]{fig01_2_small.pdf}
  \includegraphics[width=\textwidth]{fig01_3_small.pdf}
  \caption{Exposure-corrected background-subtracted images of the
    Perseus cluster (top row), M\,87 (centre row) and A\,3667 (bottom
    row) in the 0.5 to 7 keV band. The left panels include the point
    sources, while they are cosmetically hidden in the right
    panels. The colour bar units are $\log_{10} \hbox{$\photon\cm^{-2}\s^{-1}\,$}
    \rm pixel^{-1}$. The Perseus and M\,87 images were smoothed by a
    Gaussian of 2 pixels (0.984 arcsec), while A\,3667 was smoothed by
    4 pixels. North is to the top and east is to the right in all the
    images in this paper. The white dashed boxes
      show the range in areas shown in the individual scale images.}
  \label{fig:expcorr}
\end{figure*}

We downloaded the data sets listed in Tables \ref{tab:per_obs},
\ref{tab:m87_obs} and \ref{tab:a3667_obs} from the \emph{Chandra}
archive. The data sets were reprocessed using
\textsc{ciao} \citep{Fruscione06} version 4.7
and \textsc{caldb} version 4.4.10. We filtered bad time periods using
X-ray lightcurves. For observations taken with ACIS-S (the
S subarray on the Advanced CCD Imaging Spectrometer), we extracted lightcurves from
CCD 5, if available, otherwise CCD 7. When analysing ACIS-I data, we
used CCDs 0, 1 and 2. The bad time periods were chosen using an
iterative $\sigma$ clipping algorithm, clipping 200s periods with rates
outside $2.5\sigma$ where $\sigma$ is the Poisson error on the mean
count rate.  The individual observations were reprojected to the
coordinate system of the marked reference observations,
indicated
with a $^{*}$ in Tables \ref{tab:per_obs}, \ref{tab:m87_obs} and
\ref{tab:a3667_obs}. Images were constructed for each
observation and CCD using single-detector pixel binning.

Spatial masks were applied to the data to improve the quality
of the output images. We excluded the outer edges of
the ACIS-I array and those of individual ACIS-S CCDs. It is unclear why
these steps should be necessary, but doing these removed structures
associated with the edges. For M\,87, we also applied narrow exclusion
regions along the CCD read-out direction to
remove the readout streak associated with out-of-time events. In the
Perseus cluster, we masked out some regions which were far off-axis in
some observations but covered by other observations with much better
PSFs. For example, we masked out the inner few arcmin in the offset
ACIS-I observations which were covered by on-axis ACIS-S observations.

We created exposure maps for each observation and CCD using
\textsc{mkexpmap}. As input spectra for the exposure map calculation,
we assumed for Perseus a $6$~keV plasma, with a metallicity of
$0.5\hbox{$\thinspace \mathrm{Z}_{\odot}$}$, absorption equivalent to a Hydrogen column of $10^{21}
\hbox{$\cm^{-2}\,$}$ and a redshift of $0.0183$ \citep{SandersPer07}. For M\,87 we
assumed a temperature of $2.2$~keV, a metallicity of $1.1\hbox{$\thinspace \mathrm{Z}_{\odot}$}$,
absorption of $1.93 \times 10^{20} \hbox{$\cm^{-2}\,$}$ \citep{Million10} and a
redshift of $0.004283$ \citep{Cappellari11}. For A\,3667 we took the
temperature to be $7$~keV and the ICM abundance to be $0.3\hbox{$\thinspace \mathrm{Z}_{\odot}$}$
\citep{Vikhlinin01}, we used redshift of $0.0556$
\citep{StrubleRood99} and the Galactic column of $4.44 \times 10^{20}
\hbox{$\cm^{-2}\,$}$ \citep{Kalberla05}.

Standard blank-sky background event files were used to remove the
instrumental and diffuse X-ray background. For each foreground
observation we identified the appropriate backgrounds for each of the
CCDs. We filtered out events which fell on bad pixels of the
respective foreground observations. The exposure time of each
background dataset was adjusted to match the foreground rate in the 9
to 12 keV band. If the exposure times for the different CCD
backgrounds were different, we reduced the exposure times to match the
shortest background, discarding random events to keep the rate the
same. The backgrounds for each CCD were then merged and reprojected to
match the aspect of the respective foreground observations. We then
adjusted these observation backgrounds to have the same ratio of
exposure time to total background as the respective foreground
observation had to the total foreground exposure, by reducing exposure
times and discarding events appropriately. The observation backgrounds
were reprojected to the reference foreground observation. Background
images for each observation and CCD were made using the same binning
as the foreground images.

Total exposure-corrected and background-subtracted images
(Fig.~\ref{fig:expcorr} left panels) were then calculated from the
foreground, background and exposure-map images. Firstly, the spatial
masks were applied to each image. The total background image was
then subtracted from the total foreground image after scaling by the
ratio of exposure times. The total background-subtracted image was
then divided by the total exposure map.  We identified point sources
in the image by eye (the usual \textsc{wavdetect} detection tool
became confused by the high surface brightness central structures in
Perseus and M\,87). To cosmetically remove these from our exposure
corrected images, we replaced pixels inside the point source regions
with random values selected from the immediately surrounding
pixels. These cosmetically-corrected images can be seen in
Fig.~\ref{fig:expcorr} (right panels).

\section{Gradient filtering}
\subsection{Single scale filtering}
The GGM filter calculates the gradient of the image assuming Gaussian
derivatives. The gradient of the image is computed along the two axes
by convolving the image by the gradient of a one-dimensional Gaussian
function. These two gradient images are combined to create a total
gradient image. We used the implementation from \textsc{scipy}
(\url{http://scipy.org/}). As inputs we used the unsmoothed
exposure-corrected background-subtracted images, with the point
sources cosmetically removed.

\subsubsection{The Perseus cluster}
\begin{figure*}
  \includegraphics[width=\textwidth]{fig02_small.pdf}
  \caption{GGM-filtered images of the Perseus cluster with $\sigma=1$,
    2, 4, 8, 16 and 32 pixels (1 pixel is 0.492 arcsec). The bar has a
    length of 4 arcmin (89 kpc). The smallest and
      largest regions shown are indicated in Fig.~\ref{fig:expcorr}.}
  \label{fig:per_scales}
\end{figure*}

Fig.~\ref{fig:per_scales} shows filtered images of the Perseus cluster
using six exponentially increasing size scales. For each value of
$\sigma$ we show different regions of the cluster as it becomes harder
to measure gradients using small $\sigma$ when the count rate becomes
lower (see Section \ref{sect:detgrads}). Using
$\sigma=1$ and $2$, the image is sensitive to the finest structures in
the centre, which include the edges of the inner cavities and the
shocked region surrounding them \citep{FabianPer06}. Furthermore, some
of the features associated with the absorbing high-velocity system and
cool X-ray emitting filaments \citep{FabianPerFilament03} are
visible. Increasing the scales to $2$ and $4$, the cool spiral looping
around the north of the cluster \citep{SandersPer04} can be seen. To
the north is the `fountain', a structure seen in soft X-rays and by
its H$\alpha$ emission \citep{FabianPer06}. In addition the ripples in
X-ray surface brightness can be seen to the eastern side of the
cluster. In the 4 and 8 maps, the outer north-west and south ghost
cavities become visible, which are pointed to by low frequency radio
spurs \citep{FabianPer02}. Further out, in the 8 and 16 maps, is the
western edge \citep{Churazov03,FabianPer11}, curving in the same
direction as the inner spiral and at a radius of
  $\sim 110$~kpc. The top of the curve is flattened and is associated
with a surface brightness depression, one of a series in that
direction \citep{SandersPer07,FabianPer11}. These features may be
ancient relics of AGN activity, indicated by the radio emission and
long H$\alpha$ filament pointing towards them. In the 8 and 16 maps to
the south is the inverted-edge known as the `bay', which may be
another accumulation of previous AGN activity. On the largest scales,
the core of the cluster appears to be contained within an 9 arcmin
radius (200 kpc) egg-shaped region.

\begin{figure*}
  \includegraphics[width=\textwidth]{fig03_small.pdf}
  \caption{GGM-filtered images of M\,87 with $\sigma=1$, 2, 4, 8, 16
    and 32 pixels (1 pixel is 0.492 arcsec). The bar has a length of 2
    arcmin. Assuming a distance of 16.1 Mpc \citep{Tonry01}, this
    corresponds to 9.4 kpc.  The smallest and
      largest regions probed are shown in Fig.~\ref{fig:expcorr}.}
  \label{fig:m87_scales}
\end{figure*}

\begin{figure*}
  \includegraphics[width=\textwidth]{fig04_1_small.pdf}
  \includegraphics[width=\textwidth]{fig04_2_small.pdf}
  \caption{(Top row) GGM-filtered images of A\,3667 with $\sigma=8$,
    16 and 32 pixels (1 pixel is 0.492 arcsec). The bar has a length
    of 4 arcmin (260 kpc).  The smallest and
      largest regions here are shown in Fig.~\ref{fig:expcorr}. The
      white box shows the region examined in the surface brightness
      and gradient profiles. (Bottom row) Filtered images of a
      Poisson realisation of an elliptical $\beta$ model fitted to the
      X-ray data, showing the large-scale gradient and noise.}
  \label{fig:a3667_scales}
\end{figure*}

\subsubsection{M\,87}
The X-ray-emitting jet in M\,87 \citep{Marshall02,WilsonYang02} is
clearly seen at the centre of the smallest scale map
(Fig.~\ref{fig:m87_scales}). Surrounding the jet lies a cocoon of
relativistic plasma which creates cavities in the X-ray emission in
the jet and counter-jet directions \citep{Young02,Forman07}. In the
$\sigma=4$ map are the two well-known arms of soft X-ray emission
\citep{Bohringer95} extending out from the cluster along the arms of
the radio source, likely metal rich material lifted by the radio
bubbles \citep{Young02,Simionescu07,Million10}.  The inner 10 kpc
contains a series of bubbly structures, particularly along the eastern
arm. The arm may be made up of a series of buoyant bubbles. In the
eastern arm the radio plasma and X-rays are cospatial, while along the
southwestern arm they are anticoincident \citep{Forman07} and appear
to spiral around each other \citep{Forman05}. Surrounding the nucleus
with a radius of $\sim 13 \rm\thinspace kpc$ is a circular structure
\citep{Young02,Forman05} which is a weak shock \citep{Forman07}. The
$\sigma=8$ pixel map shows that there are a series of edges which have
a similar curvature to the shock, particularly in the region towards
the southwest. There appears a second edge beyond the shock towards
the northeast. We also see linear structures where the arms cross the
shock. On the largest scales are outer cavities in the X-ray emission,
surrounded by an other ring \citep{Forman07}.

\begin{figure}
  \includegraphics[width=\columnwidth]{fig05.pdf}
  \caption{Surface brightness and gradient profiles in a strip across
    A\,3667 from the south-west to the north-east
    (Fig.~\ref{fig:a3667_scales}). The X-ray surface brightness in
    $\hbox{$\photon\cm^{-2}\s^{-1}\,$}$ is shown in the top panel. The other panels show
    the gradient magnitude for the scales shown in $\hbox{$\photon\cm^{-2}\s^{-1}\,$} \:
    \mathrm{pixel}^{-1}$, for 0.492 arcsec pixels.
    For comparison are plotted the surface brightness and
    filtered profiles of a Poisson
    realisation of a model (Fig.~\ref{fig:a3667_scales} bottom row).}
  \label{fig:a3667_prof}
\end{figure}

\subsubsection{A\,3667}
There is less dynamic range in the image of A\,3667 and so we only
display the filtered maps for three spatial scales
(Fig.~\ref{fig:a3667_scales} top row). Despite this, the filtering
enhances a number of structures, many of which are not obvious in the
original image. To help assess the significance of
features we show simulated images of the cluster in
Fig.~\ref{fig:a3667_scales} (bottom row). These are based on an
elliptical $\beta$ model fit to the X-ray data, which has no
structure except for a central core. A Poisson realisation of the
fit was filtered in the same way as the original data, showing the
overall gradient and noise.
To examine the significance of the features we calculated
profiles along the main axis of the cluster in surface brightness
and GGM-filtered images on different scales
(Fig.~\ref{fig:a3667_prof}), including the $\sigma=4$ scale. For
comparison we plot the profiles from the simulated model, which
reproduces the overall gradient profiles and level of noise. There
are edges seen in the data on all filtering scales which are not
seen in the filtered model.  There are significant sharp structures
seen in the $\sigma=4$ and 8 profiles which become washed out or
mixed with other structures in the larger scale maps. On larger
scales, the gradient filter reveals longer scale gradients which are
easily missed in the noise in the smaller scale maps.

The most prominent structure is the well-known cold front indicating
that the bright X-ray emitting region is moving through its
surroundings at approximately sonic speeds \citep{Vikhlinin01}. In
these deeper data we see that the edge is not perfectly smooth but
there are features along it, likely Kelvin-Helmholtz
instabilities. There is some faint X-ray emission associated with the
cD galaxy, but a bright point source lies 20 arcsec (22~kpc) to the
south-west. We were not able to fully remove this
point source given its brightness, size and complex surrounding
structure. To the south-east of the galaxy are bright well-defined
blobs, clearly seen in the original X-ray image,
one of which is a triangle-shaped region approximately 60~kpc in
size. These regions are lower in temperature in
their surroundings and have sharp edges. They appear to be in rough
pressure equilibrium with their surroundings. They could be material
stripped during the merger, but it is unclear where the material was
stripped from.

Running between the central galaxy and the triangle is a relatively
featureless region, labelled the inner plateau. Surrounding this
plateau are many linear structures, in particular behind the cold
front. These could be projected Kelvin-Helmholtz
  instabilities (see fig.~8 of \citealt{Roediger13}). We see the
edges of a large scale excess and depression, claimed to be
300-kpc-long Kelvin-Helmholtz instabilities \citep{Mazzotta02}. At the
north-west of our field, where the data quality is poorer, are long
linear edges (labelled Outer Edges). These structures can also be seen
in the deep \emph{XMM} data of \cite{Finoguenov10}.
They are not coincident with the edge of the radio
relic to the north of the cluster, but lie a few hundred kpc inside
them. They could mark the edge of a stripped tail of material, more
easily seen in the \emph{XMM} data.  To the south of the cold front
is another edge, likely the edge detected by \cite{Vikhlinin01}, which
was identified by them as a possible bow shock.

\subsection{Detecting gradients in Poisson noise images}
\label{sect:detgrads}

The ability of the GGM filter to detect a gradient in surface
brightness depends on the magnitude of the jump, the length it occurs over,
the surface brightness and
the value of $\sigma$. As the GGM filter computes a gradient
magnitude, the output must always be zero or positive. Poisson
fluctuations in the input image will therefore produce a noise signal
in the output which is not removed by averaging over area.

\begin{figure}
  \includegraphics[width=\columnwidth]{fig06.pdf}
  \caption{Fractional gradient value detected for a
    1D sharp jump in surface brightness.  The lines show the mean
    values and standard deviations as a function of surface
    brightness. (Top panel) Values using $\sigma=1$, 2, 4, 8, 16 and
    32 pixels for a 16 per cent jump. (Bottom panel) Values for jumps
    of 4, 8, 16 and 32 per cent using $\sigma=4$ pixels.}
  \label{fig:jump_grads}
\end{figure}

To assess this quantitatively, we repeatedly filtered Poisson realisations of a
simple model images with a one-dimensional jump in surface brightness
with zero width around a certain mean surface
brightness. Fig.~\ref{fig:jump_grads} (top panel) shows how the
resulting gradient value at the jump pixel
varies as a function of surface brightness
and $\sigma$ for a fixed fractional surface brightness jump.
At low count rates the gradient signal becomes
increasingly dominated by noise. At higher count rates the mean gradient
tends towards a constant value and the standard deviation decreases.
The noise signal at low count rates has the same Poisson error origin
as the noise in the real signal at high count rates, scaling as the
surface brightness to the power $-1/2$.  This noise component also
scales as $\sigma^{-1/2}$, as expected if sensitive to
the number of counts within the filter.
As the filter has a finite width and the jump
is narrow, the determined gradient varies as $1/\sigma$.
If a model with continuous gradient over a few $\sigma$ is examined,
its value is recovered for all values of $\sigma$.

The ability to detect a jump depends on the magnitude of the jump, shown
using $\sigma=4$ in Fig.~\ref{fig:jump_grads} (lower panel). Jumps become
more visible with increasing jump size and count rate.
In this one-dimensional case the count rate at which the gradient diverges
from the noise profile decreases
with the inverse of the fractional magnitude of the jump.

\subsection{Assessing the significance of structures}

\begin{figure}
  \centering
  \includegraphics[width=0.9\columnwidth]{fig07_1_small.pdf}
  \includegraphics[width=0.9\columnwidth]{fig07_2_small.pdf}
  \caption{Comparison between different methods applied to the same
    region of Perseus. Top panel: X-ray image, smoothed by a
    $\sigma=2$ Gaussian.  2nd panel: unsharp-masked image, dividing
    Gaussian-smoothed images with $\sigma=2$ and $16$ pixels. 3rd
    panel: ratio between data and average at each radius, smoothed
    using $\sigma=2$. 4th panel: GGM-filtered image ($\sigma=4$,
    hiding point sources), with scale in $\log_{10} \hbox{$\photon\cm^{-2}\s^{-1}\,$}
    \rm pixel^{-1}$, 5th panel: GGM-filtered simulated data with the
    average radial surface brightness as the real data and same colour
    scale as 4th panel.}
  \label{fig:comparereg}
\end{figure}

Although GGM filtering is a simple image analysis technique, it is
important to assess the significance of features in resulting
maps. Our ability to detect gradients depends on the
  count rate, size of gradient (in both magnitude and length) and
  $\sigma$ (Section \ref{sect:detgrads}). The simplest method to
assess significance is to compare the filtered image with raw
data. Often structures can be directly observed in the raw data by
using a colour scale well matched to the region in question, or by
blinking between the filtered and original images.
The significance of features can also be assessed by comparing them
with the strength of noise at the same radius, where the count rate
is likely to be similar.  Filtered images can also be compared to
unsharp-masked images or images with the radial average removed
to ensure that features are robust. A further
technique is to make a simulated cluster image, based on a smooth
surface brightness profile or other model, and to filter this in the
same way as the data. The real and simulated filtered images can then
be compared to assess the significance of structures
(see e.g. Figures \ref{fig:a3667_scales} and
\ref{fig:a3667_prof}).  Comparison of images with those at other
wavelengths can also confirm the existence of structures in filtered
images (see Sections \ref{sect:scales} and \ref{sect:shocks}).

We examine a small region to the north-east of the Perseus cluster
core (Fig.~\ref{fig:comparereg}), showing a smoothed X-ray image, a
GGM-filtered image, a filtered image of a simulation of the cluster
with the same radial profile, an unsharp masked image and an image
showing the fractional residuals to the average at
each radius. The magnitude of the structures in the real data is
much larger than those in the simulated filtered map, although we see
that the real data has noise in it at a similar level to the filtered
map. The noise patterns tend to be small linear structures which lie
perpendicular to the surface brightness gradient. The main gradient
component that the filter is measuring is radial and so fluctuations
due to noise are most strongly seen in the radial direction after
filtering, leading to these characteristic noise patterns. A key
method for assessing the level of noise is to look at the amount of
noise in other directions where the gradient and data quality are
similar.

The A\,3667 data (Fig.~\ref{fig:a3667_scales}) highlight that the
appearance of noise after the filtering process depends on the data
quality. In the north west parts of the 8 and 16 scale images there
are regions filled with fluctuating dark and bright structures. This
is a region in the cluster which is only covered using relatively
short exposures. Therefore one must be careful when assessing features
in filtered data where the data quality varies strongly over the
image.

\subsection{Combining images of different scales}
\label{sect:scales}

\begin{figure}
  \includegraphics[width=\columnwidth]{fig08.pdf}
  \caption{Radial profiles of the radial weights used when making the
    combined filtered images of the Perseus cluster, M\,87 and
    A\,3667.}
  \label{fig:scalefactors}
\end{figure}

Filtering on a certain length scale is applicable to a particular
region in the cluster, where the count rate is large enough to allow
the gradient to be measured. We have therefore implemented a scheme
where we add the images with different scales together, weighting each
of the images using a radial weighting scheme. We adjust a number of
control points (radius and weighting factor) and use linear
interpolation to calculate the values for intermediate radii. The
motivation for the radial weighting is that the number of counts, the
main quantity which determines our ability to measure the gradient,
mainly varies radially in clusters. The weighting procedure also
allows the magnitude of the features in the centre to be suppressed
relative to the outskirts, in order to plot them on the same image. We
constructed a graphical user interface in order to adjust the
different radial scaling factors. Fig.~\ref{fig:scalefactors} shows
the relative radial weights of the different filtered images for
the three clusters.

\begin{figure*}
  \includegraphics[width=\textwidth]{fig09_small.pdf}
  \caption{Combined GGM-filtered of the Perseus cluster, adding maps
    with $\sigma=1$ to 32 pixels with radial weighting.}
  \label{fig:per_combine}
\end{figure*}

\begin{figure*}
  \includegraphics[width=\textwidth]{fig10_small.pdf}
  \caption{Combined GGM-filtered of M\,87, adding maps with $\sigma=1$
    to 32 pixels with radial weighting.}
  \label{fig:m87_combine}
\end{figure*}

Fig.~\ref{fig:per_combine} shows the combined image for the Perseus
cluster. In the inner part of the cluster, most of the signal comes
from the $\sigma=1$ and $2$ maps, while in the outskirts the
$\sigma=8$, $16$ and $32$ maps are
combined. Fig.~\ref{fig:m87_combine} shows the combined results for
M\,87.

\begin{figure}
  \includegraphics[width=\columnwidth]{fig11_small.pdf}
  \caption{GGM-filtered X-ray image of M87 from
    Fig.\ref{fig:m87_combine} (red), overlayed with 90~cm radio
    emission from \protect\cite{Owen00} (blue/green).}
  \label{fig:m87_radio}
\end{figure}

Fig.~\ref{fig:m87_radio} compares the filtered X-rays with
the 90~cm radio emission \citep{Owen00}, highlighting the close
connection between the X-ray and radio-emitting plasmas.
The radio source has
two arms, one to the east with a `mushroom' appearance and the other
to the south west with a filamentary structure. The southwestern arm
has a twisted appearance, where the X-rays and radio appear to be
anticoincident, likely dominated by magnetic structures \citep{Forman07}.
The eastern arm, in contrast, has coincident X-ray and radio structures.
This arm may be made up of a series of small radio bubbles plus the
large radio torus which makes the cap of the mushroom \citep{Forman07}.
Much of the edge of large scale radio structure is coincident with edges
in the X-ray surface brightness.

\begin{figure}
  \includegraphics[width=\columnwidth]{fig12_small.pdf}
  \caption{Combined GGM-filtered image of A\,3667, adding maps with
    $\sigma=8$, 16 and 32 with radial weighting.}
  \label{fig:a3667_combine}
\end{figure}

Fig.~\ref{fig:a3667_combine} shows a map combining filtered images
with scales of $\sigma=8$, 16 and 32 pixels for A\,3667. It highlights
the sharp cold front edge and surface brightness plateau between the
triangle structure and the central galaxy. The combined multi-scale
images demonstrate that it is possible to create useful qualitative
maps which show the majority of the surface brightness edges in a
galaxy cluster image.

\subsection{Comparison with other methods}
The GGM-filtered image of the example region
(Fig.~\ref{fig:comparereg}) clearly reveals the surface brightness
edges in the data, particularly those associated with the ripple-like
structures.

These ripples are not so easily seen in the image
showing the fractional deviations from the radial average, likely
due to them being mostly azimuthal in morphology. When subtracting
models as here, the choice of a model can greatly influence the
produced residual image \citep{SandersAWM712}, so care must be taken
to not introduce or remove important structures.  It is difficult to
construct models which replicate the complex structure seen, for
example, as the spiral in the core of the Perseus cluster.

Unsharp masking does show the ripples here, but they have a relatively
unclear and patchy appearance. This is similar to the results of an
analysis of images of the Centaurus cluster, where the linear
structures present are revealed using GGM filtering, but not using
unsharp masking \citep{SandersCent16}. A disadvantage of gradient
filtering over unsharp masking, however, is that intrinsically narrow
filamentary structures are broadened, such as the filament to the
north (on the western side of the image). A filament is converted to a
double structure due to the sharp gradient on either side and the flat
gradient along its ridge.

\begin{figure}
  \includegraphics[width=\columnwidth]{fig13_small.pdf}
  \caption{Comparison of unsharp masking and gradient filtering in the
    central region of M\,87. Top left panel: X-ray image with point
    sources removed and smoothed by a Gaussian of $\sigma=1$
    pixel. Top right panel: unsharp masking, showing fractional
    difference between images smoothed by 1 and 4 pixels. Centre left
    panel: unsharp masking using 2 and 8 pixels. Other panels:
    gradient filtered image with $\sigma$ as value given in pixels.}
  \label{fig:m87_compare}
\end{figure}

To compare unsharp masking and gradient filtering in more detail, we
show in Fig.~\ref{fig:m87_compare} the highly-structured central
region of M\,87. Shown are the smoothed data, two unsharp-masked
images (using different large and small smoothing scales) and three
GGM-filtered images (using scales of 2, 4 and 8 pixels). Unsharp
masking commonly produces negative artifacts surrounding bright
sources (seen near the jet here) which are difficult to distinguish
from, or mask, AGN-generated cavities in the ICM
(labelled negative artifacts).
Immediately surrounding the jet are filamentary structures
(labelled inner filaments) visible in the GGM-filtered images which
are lost in the unsharp-masked data, due to the large negative
residual artifacts there.

Unsharp masking also produces similar negative artifacts at sharp
surface brightness edges such as shocks or cold fronts, seen here to
the north-west (labelled CF edge).
GGM-filtering does not produce similar signals which
could otherwise be misinterpreted as cavities. The negative residuals
around cold fronts in unsharp-masked images are an issue when trying
to detect cavities in more distant objects, where there can be large
changes in surface brightness. Unsharp masking relies on a smoothed
image being a good approximation for the underlying cluster emission,
which is likely not to be the case in the peaked central regions.
Some features are however clearer in the unsharp-masked images,
particularly those which are smaller than the smoothing scale on a
relatively flat background, for example the small depressions to the
north of the image (labelled outer cavities).
Nevertheless, the GGM-filtered $\sigma=2$ and $4$
images reveal a great deal of structure which is not obviously present
in the unsharp-masked data. It is clear that continuous structures are
better connected on this complex underlying cluster emission in the
GGM-filtered images than when using unsharp masking.

Despite the power of the gradient filtering method, there are relative
advantages and disadvantages of GGM filtering, unsharp masking and
model subtraction. Therefore, a combination of the various techniques is
likely to help reliably identify structures, particular at higher
redshifts.

\subsection{Logarithmic gradients}
The average surface brightness profile for a galaxy cluster can be
approximated by a powerlaw in radius, at least over large radial
regions. Ideally to identify where there are deviations from a smooth
profile, it would be better to examine the gradient of the logarithm
of the surface brightness in logarithmic radius, rather than in linear
coordinates. We can partially achieve this by computing the logarithm
of the surface brightness and applying a gradient filter to the
resulting image. As X-ray images are noisy and there are often pixels
with zero values, a better way of doing this is to first smooth the
X-ray image with the Gaussian of the scale required, take the
logarithm of the smoothed image, and then compute the pixel-by-pixel
gradient magnitude.

\begin{figure}
  \includegraphics[width=\columnwidth]{fig14_small.pdf}
  \caption{Gradient image of Perseus ($\sigma=8$ pixels) (top panel)
    compared to a gradient of the log value of an image smoothed using
    a Gaussian of $\sigma=8$ pixels (bottom panel).}
  \label{fig:perlog}
\end{figure}

This method can produce good results, as shown by the filtered image
of Perseus in Fig.~\ref{fig:perlog}. The component of the gradient
which comes from the cluster profile is significantly reduced. The
image shows the structure from the inner shock out to the large scale
spiral using a single filtering scale.

A disadvantage of the technique is that the value becomes noisy where
the count rate is low in the outskirts. When the gradient of a
non-logarithmic X-ray image is computed, the pixels with low number of
counts typically are in regions with low absolute gradients and so the
noise on the gradient is low compared to the gradient value in the
centre where the counts rate are high. However, using a logarithmic
image the gradient in the outskirts is similar to the value in the
centre. As the count rate in the outskirts is lower, the scatter in
the value is higher. This noise can be seen in the north-east and
south-west parts of Fig.~\ref{fig:perlog}, where the observation is
shallower and the cluster fainter.

Logarithmic gradient images are therefore likely preferred if there
is a sufficiently high count rate across the region of interest
so that the gradient can be measured to a high fractional accuracy.
However, this criterion is unlikely to be met using typical
photon-starved observations with \emph{Chandra} except in the core
region or using large spatial scales.

\subsection{Finding shocks}
\label{sect:shocks}
As pointed out by \cite{Forman07}, between temperatures of around 1
and 3 keV the \emph{Chandra} 3.5 to 7.5 keV band count rate is
approximately proportional to the pressure-squared integrated along
the line of sight. Therefore, by gradient filtering such images, we
are able to detect pressure discontinuities and shocks in clusters.

\begin{figure*}
  \includegraphics[width=\textwidth]{fig15_small.pdf}
  \caption{Pressure-jump-sensitive GGM-filtered hard-band 3.5 to 7.5
    keV images of M\,87.  The images use $\sigma=4$ to 32 pixels,
    showing a set of pressure discontinuities marked by arrows.  The
    bar has a length of 1.5 arcmin (7.1 kpc).}
  \label{fig:m87_hard}
\end{figure*}

Fig.~\ref{fig:m87_hard} shows filtered images from scales of 4 to 32
pixels of M\,87 in this hard X-ray band. In the centre is an
egg-shaped region previously identified by \cite{Young02}, marked by A
in the 4 and 8 scale maps. This is likely a high pressure region
created by the current AGN outburst. Surrounding this feature is a
second edge in pressure (labelled B), seen clearly
in the 8 map and at lower significance in the 4 map. At a radius of
13~kpc is the clearest shock (labelled C and D), believed to be driven
by an earlier AGN episode approximately 14~Myr ago \citep{Forman07}.
Spectral fitting shows it to have a Mach number of
1.25 \citep{Million10}.

The most interesting aspect of this image, in agreement with the full
band image and filtered image, is that the 13~kpc shock is not a
complete circle, but breaks up into multiple edges (C, D and E), with
a further edge at lower surface brightness levels
(F). Edges F and G are stronger than the noise
at the same radius and can be seen in the unfiltered image. F lies
at the edge of the south-west radio plume. The splitting up of the C,
D, and E structure may be due to varying temperature structure along
the line of sight, affecting the sound speed. Alternatively there
could have been multiple outbursts.

\begin{figure}
  \includegraphics[width=\columnwidth]{fig16.jpg}
  \caption{Filtered 3.5 to 7.5 keV X-ray image of Perseus (red) with
    330~MHz radio emission (blue; \citealt{FabianPer02}). The X-ray
    image is a linear combination of the $\sigma=2$, 4 and 8 pixels
    maps. The image measures 4.3 arcmin (96 kpc) across.}
  \label{fig:per_hard}
\end{figure}

In the Perseus cluster the intracluster medium is too hot for these
hard-band images to be solely sensitive to pressure
variations. However, in the central region around the inner cavities
the 2.7 to 4 keV temperatures are close to the preferred
range. Fig.~\ref{fig:per_hard} shows a filtered hard X-ray image of
Perseus (combining three different scales), overlaying the radio
emission. The image highlights the jumps in pressure at the edge of
the shocks surrounding the inner cavities. The image can be
interpreted as two shocks, one surrounding each bubble (as seen in a
spectral fitting pressure map; \citealt{FabianPer06}). The southern
rim of the northern shock appears to pass through the southern cavity
and the rim of the northern rim of the southern shock appears to pass
through the northern cavity. The straight feature across the southern
cavity appears to be the edge of the shock and is not related directly
to the radio source. To the north-west the outer edge of the shock is
much less clear, where the radio plasma extends from the northern
inner cavity to the outer north-western ghost cavity.

\section{Conclusions}
We examine X-ray images of the Perseus cluster, M\,87 and A\,3667 with
the Gaussian gradient magnitude filter to detect edges. We show that
the filter is able to detect a host of structures within these
clusters. The method is often more sensitive to features than
unsharp-masking or subtracting radial cluster models. It also does not
introduce negative residual artifacts commonly seen in unsharp-masked
images.  By the use of a radial weighting scheme we can produce a
multi-scale image which demonstrates that a wealth of physical
processes are occurring in these clusters. Using pressure-sensitive
hard-energy-band images it is possible to use the method to detect
shocks in clusters.

\section*{Software repository}
The code described and used in this paper can be found at
\url{https://github.com/jeremysanders/ggm}.

\section*{Acknowledgements}
ACF, HRR and SAW acknowledge support from the ERC Advanced Grant
FEEDBACK. The scientific results reported in this article are based on
data obtained from the Chandra Data Archive.

\bibliographystyle{mnras}

redpajama_set_name: RedPajamaArXiv