text:
\section{Introduction}
Variability is one of the most fundamental features in data analysis. Statistics provides us with many tools suitable for measuring variability, the most popular of which is the variance. However, it can be argued that in some cases, for example when dealing with a distribution that significantly deviates from the Gaussian curve, variance  misses to convey some important information. Therefore, different metric have been proposed, which also provide quantitative information if the variance is infinite.

A particular metric known as the Gini mean difference (abbreviated to GMD), first proposed in 1912 by the Italian statistician, Corrado Gini. Since then GMD and the parameters that can be derived from it (such as the Gini coefficient) have received quite some interest. However, they are mainly applied to quantify and compare the income disparities among countries while it is being suggested that this application does not reflect the entire potential of the GMD.

Out of the Gini family, the most well-known index is the Gini coefficient, again mostly used in the area of income inequalities. It was developed independently of the GMD and defined with respect to the Lorenz curve (explained in more detail further), which definition until today remains the most well-known formulation of the Gini coefficient. In 1914, Corrado Gini discovered and documented the relationship between the Gini coefficient and the Gini mean difference which turns out to be very simple: the Gini coefficient equals the GMD divided by twice the mean, assuming that the latter is positive. As a result, the Gini coefficient ranges from 0 to 1, where 0 corresponds to perfect equality and 1 to the situation where only one observation is positive.

So far at least 14 distinct alternative representations have been described, each of them leading to a different interpretation. The readers interested in investigating these various representations further should refer to \cite{Gini_methodology} or for a more concise guide to \cite{dozen_ways}.

Our motivation for considering the Gini coefficient arises from the field of epidemiology, where it has been realized that, apart from the mean reproduction number, variability can also play an important role. Viruses such as COVID-19 are considered to be highly overdispersed (see \cite{Taleb2} and \cite{evidence_fat_tailed}). So far, describing overdispersion was restricted to a particular parametric model (see \cite{Lloyd}). We want to make the case that the Gini coefficient provides a measure of overdispersion that is not limited to any particular probability distribution, and is consistent with the insights provided by existing parametric models.

As computing the Gini coefficient with its classic definition might by arduous for some probability distributions, we first add two new representations: a probabilistic and an analytical one. We then apply the newly invented methodology to compute the Gini coefficient for a few classic probability distributions. We link our findings with the potential application to measuring overdispersion by giving more attention to the negative binomial distribution, which is often used in epidemiology to model the number of secondary cases per infected. In particular, we examine its asymptotic behavior.

This paper is organised as follows: in Section \ref{section_new_rep} we derive two new representations of the Gini coefficient: a probabilistic and an analytical one. In Section \ref{gini_various_distr} we apply the newly invented methodology to compute the Gini coefficient for a few classic probability distributions: the Poisson distribution in Section \ref{Poi_distr}, the exponential distribution in Section \ref{exp_distr}, the geometric distribution in Section \ref{geom_distr}, the Pareto distribution in Section \ref{Pareto_distr}, the uniform distribution in Section \ref{unif_distr} and the negative binomial distribution in Section \ref{negbinom_distr}. When possible, we refer to other examples of computing the Gini coefficient for the mentioned distributions. In Section \ref{applications} we focus on applications of the Gini coefficient for the negative binomial distribution, which is very often used for epidemiological modeling. We explain why the Gini coefficient should be considered as a measure of overdispersion in infectious diseases. In Theorems \ref{lim_k_small} and \ref{lim_k_large}, we prove the asymptotic behavior of the Gini coefficient of the negative binomial distribution for small and large parameters, respectively. Lastly, in Section \ref{conclusion} we close with a conclusion.

\section{New representations of the Gini coefficient} \label{section_new_rep}

The Gini coefficient is originally defined based on the Lorenz curve, which plots the proportion of the total income of the population ($y$-axis) that is cumulatively earned by the bottom $x \%$ of the population. Thus, the line at $45$ degrees represents perfect equality of incomes. The Gini coefficient can then be thought of as the ratio of the area that lies between the line of equality and the Lorenz curve over the total area under the line of equality. Before its link with the Gini mean difference was discovered, the coefficient was actually called a `concentration ratio'. However, the graphical representation is not always that useful. In order to derive alternative expressions of the Gini index, we first recall its well-known \cite{Gini_methodology, dozen_ways} definition in terms of the mean difference.

\begin{definition}[Gini coefficient] \label{def_gini_wiki}
The stochastic representation of the Gini coefficient of a random variable $X$ with finite mean $\mathbf{E}[X]$ is
\begin{equation} \label{Gini_classic}
    G(X) = \frac{\mathbf{E}[{|X_1 - X_2|}]}{2\mathbf{E}[X]},
\end{equation}
where $X_1$ and $X_2$ are independent copies of $X$.
\end{definition}

\subsection{Representation in terms of a excess random variable.}

Formula (\ref{Gini_classic}) is the most popular representation of the Gini coefficient. However, it turns out that for many probability distributions it yields a rather unpleasant form. Hence, we propose a more instructive representation in terms of excess random variable:

\begin{definition}[excess random variable]
Let $X$ be a non-negative random variable with $0<\mathbf{E}[X]<\infty$. We define $X^*$, the excess version of $X$, as
\begin{align} \label{def_size_biased}
    \mathbf{P}(X^* \geq x) = \frac{1}{\mathbf{E}[X]} \int_x^{\infty} \mathbf{P}(X > s) \hspace{0.1cm} \text{d}s, \hspace{0.4cm} x \in \mathbb{R}_{\geq 0}
\end{align}
\end{definition}

\begin{theorem} [Gini coefficient in terms of excess variable] \label{thm_size_biased_representation}
For $X\geq0$ and $X^*$ independent,
\begin{align} \label{gini_formula}
    G(X) = \mathbf{P}(X^* \geq X),
\end{align}
where $X^*$ is a excess random variable defined in (\ref{def_size_biased}).
\end{theorem}

\begin{proof}
By the definition of the Gini coefficient stated in (\ref{Gini_classic}),
\begin{align} \label{gini_rearrange}
    G(X) & = \frac{1}{2 \mathbf{E}[X]} \int_{0}^{\infty} \int_{0}^{\infty} |y-x| \hspace{0.1cm} \text{d}F_X(y) \hspace{0.1cm} \text{d}F_X(x) \nonumber \\
    & = \frac{1}{\mathbf{E}[X]} \int_{0}^{\infty}  \int_{x}^{\infty} (y-x) \hspace{0.1cm} \text{d}F_X(y) \hspace{0.1cm} \text{d}F_X(x).
\end{align}
We write $\int_{x}^{\infty} (y-x) \text{d}F_X(y) = \int_{x}^{\infty} \int_x^y \text{d}z \hspace{0.1cm} \text{d}F_X(y)$ and since the integrands are non-negative we can apply Fubini's theorem to change the order of integration in the latter and rearrange (\ref{gini_rearrange}) as
\begin{align*}
    G(X) & = \frac{1}{\mathbf{E}[X]} \int_{0}^{\infty}  \int_x^{\infty} \int_z^{\infty} \text{d}F_X(y) \hspace{0.1cm} \text{d}z \hspace{0.1cm} \text{d}F_X(x) \\
    & = \int_{0}^{\infty}  \int_x^{\infty} \frac{\mathbf{P}(X > z)}{\mathbf{E}[X]} \hspace{0.1cm} \text{d}z \hspace{0.1cm} \text{d}F_X(x) = \mathbf{P}(X^* \geq X).
\end{align*}
\end{proof}

Denote $\mathbb{N}_0 = \mathbb{N} \cup \{0\}$. Theorem \ref{thm_size_biased_representation} is valid both for discrete and continuous random variables. However, for an $\mathbb{N}_0$-valued random variable, $X^*$ in (\ref{gini_formula}) can be replaced by a discrete version of excess random variable, $X^*_d$, with the probability distribution
\begin{align*}
    \mathbf{P}(X^*_d = x) = \frac{\mathbf{P}(X>x)}{\mathbf{E}[X]}, \hspace{0.5cm} x = 0, 1, 2, 3 \ldots
\end{align*}

\begin{proposition}
Take a random variable $X$ taking values in $\mathbb{N}_0$. Then, $X^* \stackrel{d}{=} X^*_d + Y$ with $Y \sim Unif[0,1]$, independent of $X_d^*$. Consequently, $G(X) = \mathbf{P}(X^*_d \geq X)$.
\end{proposition}

\begin{proof}
Observe first that $X^*_d + Y$ has a density function $f_{X^*_d+Y}$ such that
\begin{align*}
    f_{X^*_d+Y}(z) = \frac{\mathbf{P}(X>x)}{\mathbf{E}[X]} = p_x \hspace{0.5cm} \text{for} \hspace{0.2cm} z \in [x,x+1).
\end{align*}
We now differentiate (\ref{def_size_biased}) to obtain the density of $X^*$:
\begin{align*}
    f_{X^*}(z) & = \frac{\text{d}}{\text{d}z} \mathbb{P}(X^* \leq z) = \frac{\mathbf{P}(X>z)}{\mathbf{E}[X]} = p_x \hspace{0.3cm} \text{for} \hspace{0.1cm} z \in [x,x+1).
\end{align*}
Thus, $X^* \stackrel{d}{=} X^*_d + Y$. Next, observe that, for all $\mathbb{N}_0$-valued $X$,
\begin{align} \label{show_equality}
    \mathbf{P}(X^* \geq X) = \sum_{x=0}^{\infty} \mathbf{P}(X = x) \mathbf{P}(X^* \geq x) = \sum_{x=0}^{\infty} \mathbf{P}(X = x) \mathbf{P}(X^*_d + Y \geq x).
\end{align}
Since we only consider $x \in \mathbb{N}_0$, it is true that
\begin{align*}
    \mathbf{P}(X^*_d + Y \geq x) = \mathbf{P}(X^*_d \geq x) \cdot 1 + \mathbf{P}(X^*_d = x-1) \cdot \mathbb{P}(Y=1) = \mathbf{P}(X^*_d \geq x).
\end{align*}
Substituting this into (\ref{show_equality}) completes the proof of the second statement.
\end{proof}

Since using $X^*_d$ to compute the Gini index simplifies some calculations, throughout the paper we will use this definition when working with $\mathbb{N}_0$-valued random variables.

\subsection{Representation in terms of the Fourier Transform.}

By noticing that the Gini coefficient is a convolution at $0$ of the distribution of $X$ and the tail distribution of $-X^*$, one can arrive at a different representation using Fourier theory. Since Fourier theory is different in the discrete and continuous world, we make a distinction between the two cases.

\subsubsection{Discrete integer-valued case.}

The following theorem provides a representation of the Gini coefficient
of a non-negative integer-valued random variable.

\begin{theorem} [Fourier representation of Gini coefficient in discrete case] \label{thm_Fourier_discrete}
Let $X$ be a nonnegative integer-valued random variable with probability mass function $p_X(j) = \mathbf{P}(X=j)$ and Fourier transform $\hat{p}_X(\theta) = \sum_{j=0}^{\infty} \mathbf{P}(X=j) e^{i \theta j}$. Then
\begin{align} \label{formula_one}
    G(X) = \frac{1}{4\pi \mathbf{E}[X]} \int_{0}^{2\pi} \frac{[1 - \hat{p}_X(\theta) - \mathbf{E}[X](e^{-i\theta} - 1)] \cdot \hat{p}_X(-\theta) }{1 - \cos{\theta}} \hspace{0.1cm} \text{d}\theta.
\end{align}
\end{theorem}

\begin{remark}
Since $ \hat{p}_X(\theta) \hat{p}_X(-\theta) = ||\hat{p}_X(\theta)||^2 $, (\ref{formula_one}) can be alternatively written as
\begin{align} \label{formula_one_alternative}
    G(X) = \frac{1}{4\pi \mathbf{E}[X]} \int_{0}^{2\pi} \frac{[1 - \mathbf{E}[X](e^{-i\theta} - 1)] \cdot \hat{p}_X(-\theta) - ||\hat{p}_X(\theta)||^2 }{1 - \cos{\theta}} \hspace{0.1cm} \text{d}\theta.
\end{align}
\end{remark}

\begin{proof}
Write
\begin{align} \label{Gini_as_convolution}
    G(X) & = \sum_{j=0}^{\infty} \mathbf{P}(X=j) \mathbf{P}(X^*_d \geq j) = \sum_{j=0}^{\infty} \mathbf{P}(-X=-j) \mathbf{P}(X^*_d \geq j).
\end{align}
Now, if we define the functions $f_X$ and $g_X$ as
\begin{align*}
    f_X(j) = \mathbf{P}(-X = j) \hspace{0.3cm}  \text{and} \hspace{0.3cm} g_X(j) = \mathbf{P}(X^*_d \geq j),
\end{align*}
then we notice that (\ref{Gini_as_convolution}) is a convolution of these functions, say $(f_X \star g_X)(n)$, at $n=0$, so $G(X) = (f_X \star g_X)(0),$ which can be rearranged with the help of Fourier theory. Applying the inverse Fourier transform and the convolution theorem we obtain
\begin{align} \label{inverse_fourier}
    (f_X \star g_X)(n)
    & = \frac{1}{2\pi} \int_0^{2\pi} \hat{f_X}(\theta) \cdot \hat{g_X}(\theta) \cdot e^{-in \theta} \hspace{0.1cm} \text{d} \theta,
\end{align}
Substituting $n=0$ in (\ref{inverse_fourier}), we arrive at the new representation of $G(X)$ as
\begin{align} \label{gini_in_terms_of_Fourier}
    G(X) = \frac{1}{2\pi} \int_0^{2\pi} \hat{f}_X(\theta) \cdot \hat{g}_X(\theta) \hspace{0.1cm} \text{d}\theta.
\end{align}
It remains to compute the Fourier transforms $\hat{f}_X(\theta)$ and $\hat{g}_X(\theta)$. It holds that
\begin{align} \label{Fourier_f_theta_1}
    \hat{f}_X(\theta) = \sum_{j=-\infty}^{0} e^{i \theta j} \mathbf{P}(X=-j) = \hat{p}_X(-\theta).
\end{align}
Further,
\begin{align*}
    \hat{g}_X(\theta) = \sum_{j=0}^{\infty} e^{i \theta j} \sum_{m=j}^{\infty} \mathbf{P}(X^*_d = m) = \sum_{m=0}^{\infty} \mathbf{P}(X^*_d = m) \frac{e^{i(m+1)\theta}-1}{e^{i \theta}-1},
\end{align*}
where the last equality follows from changing the order of summation and substituting $ \frac{e^{i(m+1)\theta}-1}{e^{i \theta}-1}=\sum_{j=0}^m  e^{i \theta j}$. Now, after extracting $\frac{1}{e^{i \theta}- 1}$ and splitting the sum we arrive at
\begin{align} \label{Fourier_g_theta_helpline3}
    \hat{g}_X(\theta) & = \frac{e^{i \theta}}{e^{i \theta}- 1} \sum_{m=0}^{\infty} \mathbf{P}(X^*_d = m) e^{i \theta m} - \frac{1}{e^{i \theta}- 1} \sum_{m=0}^{\infty} \mathbf{P}(X^*_d = m) \nonumber \\
    & = \frac{e^{i \theta}\hat{p}_{X^*_d}(\theta) - \hat{p}_{X^*_d}(0)}{e^{i \theta}- 1}=\frac{e^{i \theta}\hat{p}_{X^*_d}(\theta) - 1}{e^{i \theta}- 1}.
\end{align}
We now compute $\hat{p}_{X^*_d}(\theta)$ as
\begin{align*}
    \hat{p}_{X^*_d}(\theta) & = \sum_{m=0}^{\infty} e^{im\theta} \frac{\mathbf{P}(X > m)}{\mathbf{E}[X]} = \frac{1}{\mathbf{E}[X]} \sum_{m=0}^{\infty} e^{im\theta} \sum_{l>m}^{\infty} \mathbf{P}(X = l) \\
    & = \frac{1}{\mathbf{E}[X]} \sum_{l=1}^{\infty} \mathbf{P}(X = l) \frac{e^{i \theta l}-1}{e^{i \theta}-1}=\frac{1}{\mathbf{E}[X]} \sum_{l=0}^{\infty} \mathbf{P}(X = l) \frac{e^{i \theta l}-1}{e^{i \theta}-1},
\end{align*}
where we again changed the order of summation and computed $\sum_{m=0}^{l-1} e^{im\theta}$ to obtain the last equality. After extracting $\frac{1}{e^{i \theta}-1}$ and splitting the sum we have
\begin{align} \label{Fourier*helpline2}
    \hat{p}_{X^*_d}(\theta) & = \frac{1}{\mathbf{E}[X] (e^{i \theta}-1)} \bigg( \sum_{l=0}^{\infty} \mathbf{P}(X = l)e^{i \theta l} - \sum_{l=0}^{\infty} \mathbf{P}(X = l) \bigg) \nonumber \\
    & = \frac{\hat{p}_X(\theta) - 1}{\mathbf{E}[X] (e^{i \theta}-1)}.
\end{align}
Substituting (\ref{Fourier*helpline2}) into (\ref{Fourier_g_theta_helpline3}) yields
\begin{align*}
    \hat{g}_X(\theta) & =  \frac{e^{i \theta}\hat{p}_{X^*_d}(\theta) - 1}{e^{i \theta}- 1} = \frac{e^{i \theta}\frac{\hat{p}_X(\theta) - 1}{\mathbf{E}[X] (e^{i \theta}-1)} - 1}{e^{i \theta}- 1} \\
    & = \frac{e^{i\theta}}{(e^{i\theta}-1)^2} \frac{\hat{p}_X(\theta) - 1}{\mathbf{E}[X]} - \frac{1}{e^{i\theta}-1}.
\end{align*}
Finally, we  rewrite
    \[
    \frac{e^{i\theta}}{(e^{i\theta}-1)^2}
    =\frac{1}{(e^{i\theta}-1)(1-e^{-i\theta})}
    =-\frac{1}{||e^{i\theta}-1||^2},
    \]
so that
\begin{align} \label{Fourier_g_theta_final_1}
    \hat{g}_X(\theta) = \frac{1-\hat{p}_X(\theta)-\mathbf{E}[X](e^{-i\theta}-1)}{||e^{i\theta}-1||^2 \mathbf{E}[X]}.
\end{align}
We derive the final formula from Theorem \ref{thm_Fourier_discrete} by substituting (\ref{Fourier_f_theta_1}) and (\ref{Fourier_g_theta_final_1}) into (\ref{gini_in_terms_of_Fourier}) and noticing that $||e^{i\theta}-1||^2=\sin(\theta)^2+(\cos(\theta)-1)^2=2[1-\cos(\theta)].$
\end{proof}

\subsubsection{Continuous case.}
The following theorem provides an alternative representation of the Gini coefficient of a continuous non-negative random variable:

\begin{theorem} [Fourier representation of Gini coefficient in continuous case] \label{gini_continuous_fourier}
Let $X$ be a non-negative random variable with probability density function $p_X(x)$ and Fourier transform $\hat{p}_X(\theta) = \int_{0}^{\infty} p_X(x) e^{i\theta x} \hspace{0.1cm} \text{d}x$. Then
\begin{align} \label{gini_continuous_fourier_formula}
    G(X) = \frac{1}{2\pi \mathbf{E}[X]} \int_{\mathbf{R}} \frac{ \hat{p}_X(\theta)  (1 - \hat{p}_X(-\theta) - i \theta \mathbf{E}[X])}{\theta^2} \hspace{0.1cm} \mathrm{d}\theta.
\end{align}
\end{theorem}

\begin{proof}
The main idea in this proof is exactly the same as it was in the discrete case, though we need to define functions $f_X$ and $g_X$ in a slightly different way. Consider a convolution of $f_X(x), g_X(x)$ with
\begin{align*}
f_X(x) = p_X(x) \hspace{0.2cm} \text{and} \hspace{0.2cm} g_X(x) = \mathbf{P}(-X^*\leq x) = \mathbf{P}(X^* \geq -x),
\end{align*}
where $g_X$ is defined for negative values of $x$. Applying analogous reasoning as in the proof of Theorem \ref{thm_Fourier_discrete}, we get
\begin{align} \label{convolution_formula}
    G(X) = \frac{1}{2\pi} \int_{\mathbf{R}} \hat{f}_X(\theta) \hspace{0.1cm} \hat{g}_X(\theta) \hspace{0.1cm} \text{d}\theta.
\end{align}
We again compute the corresponding Fourier transforms. We have
\begin{align*}
    \hat{g}_X(\theta) & = \int_{-\infty}^{0} \mathbf{P}(X^* \geq -x) e^{i\theta x} \hspace{0.1cm} \text{d}x = \int_{0}^{\infty} \mathbf{P}(X^* \geq y) e^{-i\theta y} \hspace{0.1cm} \text{d}y \nonumber \\
    & = \int_{0}^{\infty} e^{-i\theta y} \int_{y}^{\infty} p_{X^*}(z) \hspace{0.1cm} \text{d}z  \hspace{0.1cm} \text{d}y = \frac{1}{\mathbf{E}[X]} \int_{0}^{\infty} e^{-i\theta y} \int_{y}^{\infty} \mathbf{P}(X > z) \hspace{0.1cm} \text{d}z  \hspace{0.1cm} \text{d}y \nonumber \\
    & = \frac{i}{\theta \mathbf{E}[X]} \int_{0}^{\infty} \mathbf{P}(X > z)(e^{-i \theta z} - 1 ) \hspace{0.1cm} \text{d}z,
\end{align*}
where the last equality follows after changing the order of integration and computing $\int_{0}^{z} e^{-i\theta y} \hspace{0.1cm} \text{d}y = \frac{i}{\theta}(e^{-i \theta z} - 1 )$. We substitute $\mathbf{P}(X > z) = \int_{z}^{\infty} f_X(y) \hspace{0.1cm} \text{d}y$ and change the order of integration again to obtain
\begin{align*}
    \hat{g}_X(\theta) = \frac{i}{\theta \mathbf{E}[X]} \int_{0}^{\infty} f_X(y) \hspace{0.1cm} \int_{0}^{y} (e^{-i \theta z} - 1 ) \hspace{0.1cm} \text{d}z \hspace{0.1cm} \text{d}y.
\end{align*}
Note that $\int_{0}^{y} (e^{-i \theta z} - 1 ) \hspace{0.1cm} \text{d}z = \frac{i}{\theta}(e^{-i \theta y} - 1) - y $. Hence,
\begin{align*}
    \hat{g}_X(\theta) & = \frac{-1}{\theta^2 \mathbf{E}[X]} \int_{0}^{\infty} f_X(y) e^{-i \theta y} \hspace{0.1cm} \text{d}y + \frac{1}{\theta^2 \mathbf{E}[X]} \int_{0}^{\infty} f_X(y) \hspace{0.1cm} \text{d}y \nonumber \\
    &\qquad - \frac{i}{\theta \mathbf{E}[X]} \int_0^{\infty} y f_X(y) \hspace{0.1cm} \text{d}y,
\end{align*}
which simplifies to
\begin{align} \label{g_final}
    \hat{g}_X(\theta) = \frac{1 - \hat{p}_X(-\theta) - i \theta \mathbf{E}[X]}{\theta^2 \mathbf{E}[X]}.
\end{align}
We obtain (\ref{gini_continuous_fourier_formula}) by substituting (\ref{g_final}) in (\ref{convolution_formula}) noting that $\hat{f}_X(\theta) = \hat{p}_X(\theta).$
\end{proof}

\section{Gini coefficient for various distributions} \label{gini_various_distr}

In this section, we apply our expressions to the computation of $G$ for various distributions. In some cases, the examples are new, and sometimes, our expressions provide easier derivations.

\subsection{Poisson distribution} \label{Poi_distr}

Take $X \sim \text{Poisson}(\lambda)$. The characteristic function equals $\hat{p}_X(\theta)=\exp(\lambda (e^{i\theta}-1)).$ Consequently, $\hat{p}_X(-\theta) = \exp(\lambda (e^{-i\theta}-1))$ and we substitute these results into (\ref{formula_one}) to obtain
\begin{align}
    G(X) = \frac{1}{4 \lambda\pi} \int_{0}^{2\pi} \frac{[1 - \exp(\lambda (e^{i\theta}-1)) - \lambda (e^{-i\theta}-1)] \cdot \exp(\lambda (e^{-i\theta}-1)) }{1-\cos(\theta)} \hspace{0.1cm} \text{d}\theta. \nonumber
\end{align}
Alternatively, we could use (\ref{formula_one_alternative}). We compute
\begin{align*}
    ||\hat{p}_X(\theta) ||^2 & = ||\exp(\lambda (e^{i\theta}-1))||^2 \nonumber \\
    & = e^{-2\lambda} ||e^{\lambda(\cos{\theta}+i \sin{\theta})}||^2 = e^{-2\lambda(1-\cos{\theta})}.
\end{align*}
Substituting yields
\begin{align}
    G(X) = \frac{1}{4 \lambda\pi} \int_{0}^{2\pi} \frac{(1-\lambda(e^{-i\theta}-1)) e^{\lambda (e^{-i\theta}-1)} - e^{-2\lambda(1-\cos{\theta})}}{1-\cos(\theta)} \hspace{0.1cm} \text{d}\theta. \nonumber
\end{align}
The excess formula $G(X) = \mathbb{P}(X^*_d \geq X)$ does not seem to be convenient in this example.
An alternative representation of $G(X)$ in terms of Bessel functions can be found in \cite{Ramasubban}.

\subsection{Exponential distribution} \label{exp_distr}

\begin{proposition}
If $X \sim \text{Exp}(\lambda)$ then $G(X) = \frac{1}{2}$.
\end{proposition}

\begin{proof}
Note that $X^* \stackrel{d}{=} X$ and apply Theorem \ref{thm_size_biased_representation}.
\end{proof}

\subsection{Geometric distribution} \label{geom_distr}
Take $X \sim$ Geom$(p)$ with $\mathbf{P}(X=j) = (1-p)^j p$. It can be shown that $\mathbf{P}(X^*_d \geq j) = (1-p)^j$. We plug it in $G(X)$ to arrive at
\begin{align} \label{gini_geom}
    G(X) = p \sum_{j=0}^{\infty} (1-p)^j (1-p)^j = \frac{p}{p(2-p)} = \frac{1}{2-p}.
\end{align}
In case of geometric distribution Fourier representation is also easy to compute and yields the same result as (\ref{gini_geom}).

Alternatively, one can consider a shifted version of the Geometric distribution following \cite{gini_formulas_bibi}:
take $Y \sim$ Geom'$(p)$ with $\mathbf{P}(Y=j) = p(1-p)^{j-1}$ supported on $\{1,2,3,...\}$. It follows that $\mathbf{P}(Y^*_d \geq j) = (1-p)^j$ and thus $G(Y) = \frac{p}{1-p} \sum_{j=1}^{\infty} (1-p)^{2j} = \frac{1-p}{2-p}.$

\subsection{Pareto distribution} \label{Pareto_distr}
Take $X \sim$ Pareto$(\alpha)$ with $f_X(x) = \frac{\alpha x_m^{\alpha}}{x^{\alpha+1}} = $, $\mathbf{E}X = \frac{\alpha x_m}{\alpha-1}$ and assume $\alpha > 1$. We first compute $\mathbf{P}(X^* \geq x)$ for $x \geq x_m$:
\begin{align*}
    \mathbf{P}(X^* \geq x) & = \frac{1}{\mathbf{E}[X]} \int_{y=x}^{+\infty} \mathbf{P}(X > y) \hspace{0.1cm}  \text{d}y = \frac{\alpha - 1}{\alpha x_m} \int_{y=x}^{+\infty} \bigg(\frac{x_m}{y} \bigg)^{\alpha} \hspace{0.1cm}  \text{d}y \nonumber \\
    & =  \frac{\alpha - 1}{\alpha} x_m^{\alpha-1} \bigg[ \frac{x^{1-\alpha}}{1-\alpha} \bigg]_{x}^{+\infty} = \frac{1}{\alpha} \bigg( \frac{x_m}{x} \bigg)^{\alpha-1}.
\end{align*}
Note that $\log\big(\frac{X}{x_m}\big) \sim Exp(\alpha)$, $\log\big(\frac{X^*}{x_m}\big)|X^*\geq x_m \sim Exp(\alpha-1)$. Since $\mathbf{P}(Y \geq Y') = \frac{\beta'}{\beta+\beta'}$, where $Y\sim Exp(\beta)$ and $Y'\sim Exp(\beta')$, we obtain
\begin{align*}
    \mathbf{P}(X^*\geq X) & = \mathbf{P}(X^*\geq X|X^*\geq x_m) \mathbf{P}(X^*\geq x_m) \\
    & = \mathbf{P}(\log\bigg(\frac{X^*}{x_m}\bigg)\geq \log\bigg(\frac{X}{x_m}\bigg)|X^*\geq x_m) \mathbf{P}(X^*\geq x_m) \\
    & = \frac{\alpha}{\alpha+\alpha-1} \cdot \frac{1}{\alpha} = \frac{1}{2\alpha-1}.
\end{align*}
This is not a new result, but its derivation is quite simple-- compare e.g. \cite[Example 1]{gini_formulas_bibi}.

\subsection{Continuous uniform distribution} \label{unif_distr}

We have not been able to find this explicit example in the literature.

\begin{proposition}
If $X \sim \mathcal{U}_{[a,b]}$, then $G(X) = \frac{b-a}{3(b+a)}$ = $\frac{2}{b-a} \cdot \frac{{\mathrm{Var}} X}{\mathbf{E}X}$.
\end{proposition}

\begin{proof}
\begin{align*}
    G(X) & = \int \mathbf{P}(X^* \geq x) f_X(x) \hspace{0.1cm} \text{d}x = \int_a^b \mathbf{P}(X^* \geq x) \frac{1}{b-a} \hspace{0.1cm} \text{d}x.
\end{align*}
Applying (\ref{def_size_biased}), substituting $\mathbf{E}X = \frac{a+b}{2}$ and changing the order of integration yields
\begin{align*}
    G(X) & = \frac{2}{(a+b)(b-a)} \int_a^{\infty} (y-a) \mathbf{P}(X > y) \hspace{0.1cm} \text{d}y \\
    & = \frac{2}{(a+b)(b-a)^2} \int_a^{b} (y-a) (b-y) \hspace{0.1cm} \text{d}y \\
    & = \frac{2}{(a+b)(b-a)^2}  \cdot \frac{1}{6} (b-a)^3 = \frac{b-a}{3(b+a)}.
\end{align*}
\end{proof}

\subsection{Negative binomial distribution} \label{negbinom_distr}

Let $X \sim NB(k,p)$ with probability mass function $\mathbf{P}(X = j) = \binom{k+j-1}{j}p^k (1-p)^j$ and $\mathbf{E}X = \frac{k(1-p)}{p}$. The classic definition of the Gini index (\ref{Gini_classic}) yields
\begin{align} \label{gini_neg_binom_ugly}
    G(X) = \frac{1}{p} \sum_{i=0}^{\infty} (-1)^i \binom{k+i}{i} \Bigg( \frac{1-p}{p^2} \Bigg)^i \frac{(2i)!}{i!(i+1)!}.
\end{align}
For a derivation see \cite{Ramasubban}, formula (2.11) for  the mean absolute difference of the negative binomial; then divide by twice the mean to obtain (\ref{gini_neg_binom_ugly}). The complexity of this expression was in fact our main motivation to develop new representations. Because of the alternating sign the formula is not very stable numerically and thus not handy to work with. The representation in terms of the excess random variable does not seem very helpful either. The computation leads to a rather non-intuitive formula containing three infinite sums. Instead, we apply Theorem \ref{thm_Fourier_discrete} to obtain the following representation:
\begin{theorem}
Take $X \sim NB(k,p)$ with probability mass function $\mathbf{P}(X = j) = \binom{k+j-1}{j}p^k (1-p)^j$ and $\mathbf{E}X = \frac{k(1-p)}{p}$, where $p$ is the probability of success. Then
\begin{align} \label{Gini_NB}
    & G(X) \nonumber = \\
    & \frac{p}{k(1-p)} \frac{1}{2\pi} \int_{0}^{2\pi} \frac{[1 - \frac{k(1-p)}{p} (e^{-i\theta}-1)] \cdot \bigg( \frac{pe^{i\theta}}{e^{i\theta}+p-1} \bigg)^k - \bigg( \frac{p^2}{p^2 + 2(p-1)(\cos{\theta}-1)} \bigg)^k }{ ||e^{i\theta}-1||^2} \hspace{0.1cm} \text{d} \theta.
\end{align}
\end{theorem}

\begin{remark}
Note that since $\binom{k+j-1}{j}$ can be interpreted as $\frac{\Gamma(k+j)}{\Gamma(j+1)\Gamma(k)}$, (\ref{Gini_NB}) also applies to non-integer $k$.
\end{remark}

\begin{proof}
We want to apply Theorem \ref{thm_Fourier_discrete} to the negative binomial variable $X$. Since
\begin{align} \label{term1}
    \hat{p}_{X}(\theta) = p^k (1 - (1-p)e^{i\theta})^{-k} = \frac{p^k}{(pe^{i\theta}-e^{i\theta}+1)^k},
\end{align}
then
\begin{align} \label{term2}
    \hat{p}_{X}(-\theta) = \frac{p^k}{(pe^{-i\theta}-e^{-i\theta}+1)^k} = \bigg( \frac{pe^{i\theta}}{e^{i\theta}+p-1} \bigg)^k.
\end{align}
Substituting (\ref{term1}) and (\ref{term2}) into (\ref{formula_one}) yields
\begin{align} \label{gini_neg_binom}
    & G(X)\nonumber \\
    & =\frac{p}{k(1-p)} \frac{1}{2\pi} \int_{0}^{2\pi} \frac{[1 - \frac{p^k}{(pe^{i\theta}-e^{i\theta}+1)^k} - \frac{k(1-p)}{p} (e^{-i\theta}-1)] \cdot \bigg( \frac{pe^{i\theta}}{e^{i\theta}+p-1} \bigg)^k }{||e^{i\theta}-1||^2} \hspace{0.1cm} \text{d}\theta.
\end{align}
We investigate the product of $\hat{p}_{X}(\theta)$ and $\hat{p}_{X}(-\theta)$:
\begin{align*}
     \bigg( \frac{p}{pe^{i\theta}-e^{i\theta}+1} \cdot \frac{pe^{i\theta}}{e^{i\theta}+p-1} \bigg)^k = \bigg( \frac{p^2e^{i\theta}}{(pe^{i\theta}-e^{i\theta}+1)(e^{i\theta}+p-1)} \bigg)^k.
\end{align*}
However, we have that
\begin{align*}
    (pe^{i\theta}- & e^{i\theta}+1)(e^{i\theta}+p-1) = pe^{2i\theta} + p^2e^{i\theta} - 2pe^{i\theta} - e^{2i\theta} + 2e^{i\theta} + p -1 \\
    & = e^{i\theta} \big(p^2-2p+2 + p(e^{i\theta} + e^{-i\theta}) - (e^{i\theta} + e^{-i\theta}) \big) \\
    & = e^{i\theta} \big( p^2 + 2(p-1)(\cos(\theta)-1) \big).
\end{align*}
Therefore
\begin{align*}
    \frac{p^2e^{i\theta}}{(pe^{i\theta}-e^{i\theta}+1)(e^{i\theta}+p-1)} = \frac{p^2}{p^2 + 2(p-1)(\cos{\theta}-1)}.
\end{align*}
After plugging in the above computation in (\ref{gini_neg_binom}) and rearranging we arrive at the expression (\ref{Gini_NB}).
\end{proof}

In the next section, we discuss the asymptotic behavior of the Gini coefficient of the negative binomial distribution as $k \to 0$ and $k \to \infty$.

\section{Applications to epidemics} \label{applications}
In epidemiology, the basic reproduction number $R_0$ denotes the expected number of new infections directly generated by one case. If $R_0>1$, the epidemic is growing, and the epidemic
is dying our when $R_0<1$.
It has been determined that transmission patterns of SARS-COV-2 and many other viruses are highly overdispersed (see \cite{Taleb2} or \cite{evidence_fat_tailed}).  Looking at $R_0$ alone might therefore be misleading. In case of COVID-19 in particular, it is estimated that more than 80$\%$ of new infections resulted from the top 10$\%$ of most infectious individuals (see \cite{Endo}, \cite{overdispersion_Israel} or \cite{overdispersion_USA}). This means that the vast majority of infected individuals will not pass the infection on to anyone. $R_0$ should therefore be complemented  with another statistic that sheds light on possible overdispersion in the number of secondary cases.

To this end, the authors in \cite{Lloyd} assume a parametric setting, in which  the number of secondary cases follows a negative binomial distribution with parameters $k$ and $R_0$, which results in its variance being equal to $R_0(1 + R_0/k)$. In such a setting, the smaller the value of $k$, the greater heterogeneity in the distribution of secondary infections. Despite coming from a very particular model, this interpretation became so popular that multiple papers adopted the symbol $k$ and referred to it as an `overdispersion parameter' (see for instance \cite{Endo}, \cite{Chen}, \cite{Eilersen} or \cite{Gini_Indonesia}). Even non-mathematical press followed this trend, calling $k$ the key to overcoming the pandemic of COVID-19 (see \cite{OverlookedVariable}).

In this section, we investigate the relation between $k$ and the Gini index. We show numerically that the Gini index, when specialized to the negative binomial distribution, depends monotonically on $k$, and complement this with rigorous asymptotic estimates for small and large $k$. This indicates that the Gini index provides insights consistent with the ones obtained so far, without the need to rely on a parametric model. Thus, the Gini index can serve as candidate to measure the variability in the infectiousness of a disease.

This section is organised as follows: in Section \ref{Gini_asypmtotic} we first enable a continuous interpretation of parameter $k$ of the negative binomial distribution by expressing is as a composition of Poisson and Gamma processes. Next we investigate monotonicity of the Gini coefficient in $k$. To that end, we illustrate its behavior in plots and support them with statements about asymptotics of the Gini coefficient for $k \to 0$ and $k \to \infty$. In Section \ref{gini_data_example} we show for the real-world data set that the Gini coefficient provides insights consistent with using a parametric model based on $k$.

\subsection{Asymptotic behavior of the Gini coefficient} \label{Gini_asypmtotic}

In Section \ref{negbinom_distr} we showed that the Gini index for the negative binomial distribution cannot be easily obtained from (\ref{Gini_classic}), neither by its excess representation $\mathbb{P}(X^*_d \geq X).$ In this section, we investigate the Gini index for small and large values of $k$.
For that purpose it is useful to first rewrite the negative binomial distribution in terms of the composition of Poisson and Gamma processes, which has the added benefit that the resulting interpretation is natural for non-integer values of $k$.

Let $X$ follow a $Poisson(\nu)$ distribution, where the rate $\nu$ is also a random variable and it follows a $Gamma(k,\lambda)$ distribution. To maintain a clear connection to epidemiological notation, we denote $\mathbf{E}[\nu] = R_0$.

\begin{lemma} \label{neg_binom_vs_poisson_gamma}
For $X$ as above, it holds that $X \sim NB(k,p)$ with parameter $k$, corresponding to the number of `successes', and parameter $p=  \frac{\lambda}{\lambda +1}$ denoting the probability of success.
\end{lemma}
The proof of Lemma \ref{neg_binom_vs_poisson_gamma} is straightforward and is hence omitted. Note that $p= (1 + R_0/k)^{-1}$. In the following, we will write $\Gamma_k$ referring to a Gamma process determined by $\Gamma(k,\lambda)$ and $N(\Gamma_k)$ to denote a subordinated Poisson process with rate following this Gamma process.

With the new, continuous interpretation in mind we will have a closer look at the behaviour of the Gini coefficient as $k$ increases. For creating some visual intuition, Figure \ref{fig:gini_asympt_plots}  plots the Gini coefficient as a function of $k$ for different values of $p=\frac{\lambda}{\lambda+1}$.
\begin{figure}[h]
    \centering
    \begin{subfigure}[b]{0.3\textwidth}
        \centering
        \includegraphics[width=\textwidth]{gg_Gini_k_large_p=0.9.Rplot}
            \caption[]%
            {{\small $p = 0.9$}}

    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{0.3\textwidth}
        \centering
        \includegraphics[width=\textwidth]{gg_Gini_k_large_p=0.5.Rplot}
            \caption[]%
            {{\small $p = 0.5$}}

    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{0.3\textwidth}
        \centering
        \includegraphics[width=\textwidth]{gg_Gini_k_large_p=0.1.Rplot}
            \caption[]%
            {{\small $p = 0.1$}}

    \end{subfigure}
    \caption{Plots of the Gini coefficient of negative binomial distribution} \label{fig:gini_asympt_plots}
\end{figure}
Figure \ref{fig:gini_asympt_plots} suggests that the Gini coefficient of the negative binomial is decreasing in $k$. Unfortunately, a rigorous proof of the statement remains an open problem. However, we succeeded in showing that the limiting behaviour of the Gini coefficient of the negative binomial distribution aligns with the desired statement, i.e., $G(NB(k,p))$ approaches $1$ as $k$ approaches $0$ and approaches $0$ as $k$ approaches infinity. We formulate these results in the following two theorems:

\begin{theorem} [Asymptotics of the Gini coefficient for $k$ small] \label{lim_k_small}
Let $X \sim NB(k,p)$. Then, for $k \to 0$,
\begin{align} \label{as_formula_k_small}
    G(NB(k,p)) = 1 + c \cdot k + o(k),
\end{align}
with $c=2\log(p)+\frac{2-p}{1-p}\log\big(2p(1-p)\big)=2\log \bigg(\frac{\lambda}{\lambda+1}\bigg) - (\lambda+2)\log\bigg(\frac{\lambda^2+2\lambda}{(\lambda+1)^2} \bigg)$.
\end{theorem}

\begin{theorem} [Asymptotic behavior of the Gini coefficient for $k$ large] \label{lim_k_large}
Let $X \sim NB(k,p)$. Then, for $k \to \infty$,
\begin{align} \label{as_formula_k_large}
    \sqrt{k} \hspace{0.1cm} G(NB(k,p)) \longrightarrow \frac{\sqrt{1 + \lambda}}{\sqrt{\pi}}.
\end{align}
\end{theorem}

\begin{figure}[h!]
    \centering
    \begin{subfigure}[b]{0.3\textwidth}
        \centering
        \includegraphics[width=\textwidth]{gg_zoom_p=0.9.Rplot.png}
            \caption[]%
            {{\small $p = 0.9$}}

    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{0.3\textwidth}
        \centering
        \includegraphics[width=\textwidth]{gg_zoom_p=0.5.Rplot.png}
            \caption[]%
            {{\small $p = 0.5$}}

    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{0.3\textwidth}
        \centering
        \includegraphics[width=\textwidth]{gg_zoom_p=0.1.Rplot.png}
            \caption[]%
            {{\small $p = 0.1$}}

    \end{subfigure}
    \caption{Behavior of the Gini coefficient for $k \to 0$ vs.\ $y=1+ck$} \label{fig:gini for k small zoom in close}
\end{figure}

\begin{figure}[h!]
    \begin{subfigure}[b]{0.3\textwidth}
        \centering
        \includegraphics[width=\textwidth]{gg_sqrt_k_times_Gini_vs_const_for_p=0.9.Rplot.png}
            \caption[]%
            {{\small $p = 0.9$}}

    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{0.325\textwidth}
        \centering
        \includegraphics[width=\textwidth]{gg_sqrt_k_times_Gini_vs_const_for_p=0.5.Rplot.png}
            \caption[]%
            {{\small $p = 0.5$}}

    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{0.325\textwidth}
        \centering
        \includegraphics[width=\textwidth]{gg_sqrt_k_times_Gini_vs_const_for_p=0.1.Rplot.png}
            \caption[]%
            {{\small $p = 0.1$}}

    \end{subfigure}
    \caption{Behavior of the Gini coefficient times $\sqrt{k}$ for $k \to \infty$ vs.\ $y=\frac{\sqrt{1 + \lambda}}{\sqrt{\pi}}$} \label{fig:gini_blowup}
\end{figure}
Proofs of these theorems are given in the Appendix. We validate the asymptotic estimates in the theorems by plotting the Gini coefficient and its limiting functions against each other. See Figure \ref{fig:gini for k small zoom in close} for $k \to 0$ and Figure \ref{fig:gini_blowup} for $k \to \infty$.

\subsection{The Gini coefficient for COVID-19' secondary infections} \label{gini_data_example}

In this section we investigate how our formulas perform in a real-world example. Heterogeneity in the number of secondary cases receives significant attention during the currently ongoing pandemic of COVID-19. \cite{Gini_Indonesia} aims to quantify overdispersion by estimating $k$. We use the data provided there to first compute the Gini coefficient explicitly (by (\ref{Gini_classic})) and then calculate it using formulas we derived in this paper. For the latter, we substituted the parameters $p$ and $k$ obtained by the authors. The table below presents all results.
\begin{table}[h!]
\begin{center}
\begin{tabular}{ |c|c|c| }
 \hline
 \textbf{Location} & \textbf{Jakarta - Depok} & \textbf{Batan} \\
 \hline
 $R_0$ & 6,79 & 2,47 \\
 \hline
 $k$ & 0,06 & 0,2 \\
 \hline
 $p$ & 0,008 & 0,06 \\
 \hline
 $G(X)$ computed explicitly from data & 0,9213411 & 0,83191721 \\
 \hline
 $G(X)$ computed with (\ref{formula_one})& 0,9269144 & 0,8151066 \\
 \hline
 $G(X)$ computed with (\ref{as_formula_k_small}) & 0.9193061 & 0.7623808\\
 \hline
\end{tabular}
\caption{\label{tab1}Secondary cases data summary for Jakarta-Depok and Batan.}
\end{center}
\end{table}
These results show that our formulas come very close to the actual value of the Gini coefficient. It is also apparent that they reflect the dynamic of epidemic - lower values of $k$ assigned to higher overdispersion correspond to higher value of the Gini coefficient.

\section{Conclusion} \label{conclusion}

The models based on the negative binomial distribution of secondary infections have been used to estimate overdispersion in the spread of infectious diseases for a long time. However, this approach necessitates parametric models suitable only for specific situations. We argue that it is preferable to use the Gini coefficient to measure overdispersion. It correctly captures the dynamic of epidemic and is arguably more informative than popular so far parameter $k$ as it is solely known as the measure of variability and it is not bound to any specific probability distribution.

Thanks to the new representations that we presented in this paper, we can express the Gini coefficient as a simple function of the parameters of the distribution for multiple probability distributions, as long as they have finite expected value, making the Gini index a
versatile statistic not depending on parametric assumptions. In particular, it can be useful in the case of heavy-tail distributions where the variance might be infinite. In \cite{Taleb} the authors present methods that allow estimating the Gini coefficient for heavy-tailed data with infinite variance.

The applications of the newly developed representations do not need to be limited to epidemiology. The Gini coefficient was advocated as an indicator in many other areas and we hope our results could help widen its usage. For reference see for example \cite{heart_rate}, where the Gini coefficient is suggested as a measure of heart rate variability to assess the level of mental stress, or \cite{White}, where the Gini coefficient is applied to explain inequalities in resource use around the globe.

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}
Ground-layer adaptive optics (GLAO) provides partial correction of atmospheric blurring over a significantly larger field of view (several to 10 arcminutes) and over a broader wavelength range (optical and infrared) than classical AO systems \cite{Rigaut:2002}.  Science applications enabled by GLAO are thus complementary to classical AO and include extra-galactic spectroscopic surveys over a broad range of redshifts, intergalactic and circumgalactic medium studies with integral field and slit spectrographs, and stellar population studies in the Milky Way and nearby galaxies where crowding limits current sensitivities. Maunakea is ideally suited for GLAO with its thin ground layer and weak free-atmosphere turbulence \cite{Chun:2009}. Currently, none of the 8-10 m class telescopes on Maunakea have a GLAO system; however, several observatories are investigating the feasibility and designing GLAO systems for the future. In this proceeding, we explore the feasibility of a GLAO system deployed at the W. M. Keck Observatory, which could have the benefit of feeding Keck's existing seeing-limited instruments including the multi-object spectrographs, LRIS and DEIMOS in the optical and MOSFIRE in the infrared.

Nearly all seeing-limited science cases would benefit from the expected 2$\times$ better spatial resolution in the near-infrared and red-optical delivered by GLAO \cite{Abdurrahman:2018,Chun:2018}.  The ultimate question that must be addressed is whether the science gains enabled by GLAO system at Keck are timely and significant enough to warrant the cost. The science gains delivered by a Keck GLAO system can be quantified with the following ingredients. First, we need to estimate the image quality delivered by a notional GLAO system. We utilize multiple atmospheric modeling codes and mean Maunakea atmospheric conditions to simulate point spread functions (PSFs) for Keck GLAO. We also use results from the on-sky 'imaka experiment to validate these simulations. Second, we need to determine how the existing instruments would further reduce the GLAO-delivered image quality due to internal aberrations, distortions, and flexure. Thirdly, we need to propagate these image quality estimates through to science metrics for a wide range of current and future science cases. Finally, science gains must always be weighed against cost and risk. In this proceeding, we present initial results from the first 2$-$3 steps of this process.

\section{GLAO Performance: On-sky Validation and Keck Simulations}
One of the primary motivations for GLAO at Keck has been results from on-sky experiments that indicated that the turbulence above Maunakea is typically dominated by a significant, but thin, ground layer \cite{Chun:2009}. In particular, the 'imaka GLAO experiment commissioned on-sky on the UH 2.2 m telescope on Maunakea \cite{Chun:2016} has shown that the FWHM can be improved by a factor of 1.5$-$1.7$\times$ and the Noise Equivalent Area is improved by $>$2$\times$ at R-band and I-band with a guide star constellation spread over 18' \cite{Abdurrahman:2018}. Preliminary results from 2018 indicate that improvements of $>$1.3$\times$ are achievable even at B-band and V-band over an 11' field of view \cite{Chun:2018}. The GLAO PSFs are also more stable by a factor of 7-10$\times$.  The results of 'imaka are promising for a Keck GLAO system.
Furthermore, GLAO performance has been validated at the VLT, showing 2$\times$ improved and very uniform FWHM over a 10' radius \cite{Madec:2018}.

Simulations of the expected performance of the a Keck GLAO system were made using the Multi-threaded Adaptive Optics Simulator (MAOS) \cite{Wang:2012}.  The initial system configuration is for a 20$\times$20 actuator system with four sodium laser guide stars placed on the periphery of a 10'$\times$10' field of view.  A single low-order natural-guide star wavefront sensor was placed at the center of the field to provide tip, tilt, and focus sensing.  Point spread functions were generated across the 10' field for wavelengths from 0.45 to 2.2 $\mu$m from a sequence of ten independent 5000-iteration simulations.  The input optical turbulence profile is derived from the Thirty Meter Telescope project's site testing at the TMT site on Maunakea \cite{Schock:2009}. The results shown in Fig.~\ref{fig:glao_sims} were run to provide an initial estimate of the gains across the science wavelengths of interest and represent averages across the field of view.  One-standard deviation in this variation are indicated on the plot.  The improvement in the FWHM is 2$\times$ at infrared wavelengths and $\sim$1.5$\times$ at optical wavelengths ($\sim$500 nm).  Analysis of other metrics, such as the encircled energy and the noise-equivalent area, is underway.  Future simulations will focus on optimizing the image quality for the specific science instrument fields of view and wavelengths and including the effect of the conjugation altitude the adaptive secondary.

\begin{figure}
\begin{center}
\includegraphics[height=8cm]{keck_glao_sim.png}
\end{center}
\caption{ \label{fig:glao_sims}
Simulated image quality, expressed in FWHM, for Keck in seeing-limited and GLAO modes at different wavelengths ({\em top}). The ratio of the GLAO FWHM to the seeing-limited FWHM is also shown ({\em bottom}).
}
\end{figure}

\section{Examples of GLAO Science Cases}

In this section, we present a brief summary of science cases that would most benefit from GLAO fed multi-object spectroscopy at Keck. More detailed science cases for GLAO at other telescopes have previously been published for Subaru \cite{SubaruGLAO}, Gemini \cite{GeminiGLAOstudy,GeminiGLAOsci}, CTIO \cite{Tokovinin:2015}, and VLT \cite{Hibon:2016,Siebenmorgen:2011,Bacon:2006}.

Measurements of galaxy properties test our physical understanding of galaxy formation and evolution.
The sensitivity of ground-based extra-galactic spectroscopy is limited by the sky brightness.  The improved spatial resolution delivered by GLAO is sufficient to
resolve galaxies that were previously unresolved with seeing-limited observations. As a result, spectroscopic slits can be narrower, thereby reducing the solid angle of blank sky that adds noise to the galaxy spectrum.
Cutting the spectroscopic aperture down to the size of resolved star-forming regions within galaxies
also increases spectral resolution, which improves sensitivity to spectral lines. It follows that GLAO
enables spectroscopic observations of fainter galaxies.

One of the great advantages of GLAO is how the improved spatial resolution allows us to resolve finer detail in distant galaxies (Fig.~\ref{fig:exgal_scale}). The power of spatially-resolved spectroscopy in large galaxy surveys has recently been demonstrated in the MaNGA survey using seeing-limited, integral-field, multi-object spectroscopy. A Keck GLAO system can push to higher redshift and lower-mass galaxies than current IFU surveys such as MaNGA.

\begin{figure}
\begin{center}
\includegraphics[height=8cm]{angular_scale.pdf}
\end{center}
\caption{ \label{fig:exgal_scale}
Interesting physical spatial scales in galaxies (e.g. disks, bars, bulges, clumps) as a function of redshift. Lines of constant angular resolution on the sky are over-plotted. GLAO at Keck could enable new spatially-resolved surveys of galaxies 3-5$\times$ higher redshift than seeing-limited surveys.
}
\end{figure}

There are many additional GLAO science cases including:
(1) detailed studies of galactic nuclei,
(2) resolving stellar populations and star formation histories in nearby galaxies,
(3) mapping and characterizing the interstellar medium, supernovae remnants, and outflows from stars in the Milky Way,
(4) detection of exoplanet from the astrometric wobble of their host star,
(5) complete imaging of Mars, Jupiter, and other planets simultaneously with their ring and moon systems.
Most of these science cases utilize the improved spatial resolution of GLAO to resolve further detail on an object or overcome stellar confusion rather than the improved sensitivity.

\section{Instrument Feasibility}

Many of Keck's existing multiplexed spectrographs and wide-field imagers could benefit from a GLAO feed.  However, the feasibility of using existing instruments depends on many factors. First and foremost is whether the intrinsic optical quality of the telescope and instrument are sufficiently good to take advantage of the higher spatial resolution that GLAO would deliver. Figure \ref{fig:instr_fov} shows the locations and field sizes for different instruments relative to the optical axis of the telescope. We also consider other factors that might limit the science gains from a GLAO feed such as atmospheric refraction, flexure, and slit construction. These turn out to be the limiting factors in most cases for existing instruments. In sections below, we consider the feasibility of the two multi-object optical spectrographs, LRIS and DEIMOS, as well as the infrared multi-object spectrograph, MOSFIRE.

\begin{figure}
\begin{center}
\includegraphics[height=8cm]{Instrument_FoV.png}
\end{center}
\caption{\label{fig:instr_fov}
The Keck instrument fields-of-view relative to the optical axis of the telescope. Both LRIS and DEIMOS fields extend out to 10' from off the optical axis and their image quality would likely be limited by off-axis telescope aberrations in the outer parts of the field. MOSFIRE is centered on the telescope optical axis.
}
\end{figure}

\subsection{LRIS}

The multi-object, optical slit spectrograph, LRIS, is located on the optical axis of the telescope (Figure \ref{fig:instr_fov}) and thus has access to the best image quality directly from the telescope. It's 6' $\times$ 6' FoV, 0.135"/pixel plate scale in the imaging direction, and designed image FWHM of $\mathrel{\hbox{\rlap{\hbox{\lower4pt\hbox{$\sim$}}}\hbox{$<$}}} 0.24$''  should make it an ideal GLAO instrument \cite{Oke:1995}. A GLAO feed for LRIS would allow the use of smaller slits for higher spectral resolution and lower background. Typical slits in seeing-limited mode are 0.8'' - 1.0'' and we anticipate a slit width reduction to 0.5''.

First, the as-delivered image quality on LRIS has been verified with on-sky images from both the LRIS Blue and Red camera. We collected a random sample of $\sim 400$ LRIS images from observations taken after 2011 Jan from the Keck Observatory Archive. Stars were identified using DAOPHOT with a number of criteria to ensure that resolved sources and artifacts were excluded. PSF fitting was performed on each star and the median stellar PSF FWHM for the data set was calculated (Figure \ref{fig:lris_fwhm}).  The median image FWHM for all cameras and filters is 0.836''.  The discrepancy between that value and the expected value of 0.67'' (atmospheric free seeing of 0.65'' and contribution from telescope astigmatism of 0.18'') suggests the presence of 0.49'' of unaccounted-for wavefront error during typical LRIS usage. The telescope contributes $\sim$0.18'' averaged over the field due to off-axis field aberrations. Thus, an unexplained term of 0.5'' remains and is likely due to slowly time-varying focus or astigmatism and should be correctable by the GLAO system.

\begin{figure}
\begin{center}
\includegraphics[scale=0.4]{lris_fwhm.png}
\end{center}
\caption{The distribution of measured FWHM from LRIS images of dense stellar fields. Measurements come from a random sample of 126 images obtained after 2011 Jan.
\label{fig:lris_fwhm}}
\end{figure}

Using realistic GLAO PSFs with the appropriate ASM conjugate altitude, we calculate exposure time improvement factors of 10\% to 140\% (i.e., t/1.1 $-$ t/2.4) for stars within a 4' $\times$ 4' field at 8000 $\AA$.  We assumed the GLAO system parameters include 3-guide stars with a constellation diameter of 4'.  These improvement factors depend on a number of parameters, including the line width, assumptions about the distribution of correctable and non-common path wavefront error, and other factors not explored here.  The most important factor is the contribution of non-common path wavefront error.  The median image quality obtained with LRIS suggests that this factor could be as high as 0.49''.  If it is, the improvement in required exposure time seen with GLAO would be small.  However, if the error is common-path and could be sensed and corrected by the GLAO system, the improvement in required exposure time with GLAO could be as high as 150\%.

In addition to the image quality, we consider other sources of aberration and distortion. LRIS has an atmospheric dispersion compensator, so differential atmospheric refraction should not contribute significantly for a GLAO PSF. On the other hand, LRIS is a Cassegrain instrument and is know to suffer from time-variable distortion and focus. Previous studies of the astrometric distortion \cite{Cohen:2007} indicate that in seeing limited mode, the flexure amounts to 0.1'' - 0.2'' over the field. However, the same document states that flexure of more than 10\% of a slit width is not tolerable. Thus if the LRIS slit sizes are reduced from 0.8'' to 0.5'' then the LRIS flexure would not keep the science object on its slit. These flexures need to be more extensively quantified in order to determine if some real time corrections could be applied when using GLAO.

\subsection{DEIMOS}

We investigated whether ground-layer adaptive optics (GLAO) can benefit spectroscopy with the DEIMOS spectrograph on the Keck telescope. DEIMOS is a wide field multi-object spectrograph designed for long exposure extra-Galactic work, particularly observation of distant, early, galaxies. The DEIMOS science field is 16.7' $\times$ 5' with a pixel scale of 0.1185'' pixel$^{-1}$ for both the spectrograph and imager. This would be sufficient to Nyquist-sample a PSF that has been sharpened by GLAO. The intrinsic image quality within DEIMOS is has been characterized using pinhole mask observations (Figure \ref{fig:deimos_iq}) and is 0.26'' on average from the instrument. However, the telescope has an increasing astigmatism that increases with distance from the telescope's optical axis and would yield and additional 0.3'' (added in quadrature) on average over the DEIMOS field. Finally, DEIMOS has not been equipped with an atmospheric dispersion compensator. During long exposures, this could cause objects to shift outside the slit due to atmospheric refraction.

\begin{figure}
\begin{center}
\includegraphics[height=3in]{deimos_iq.png}
\end{center}
\caption{\label{fig:deimos_iq}
DEIMOS image quality as measured by the FWHM of pinhole mask spots over the field of view. }
\end{figure}

\subsection{MOSFIRE}

MOSFIRE is a infrared imager and multi-object spectrograph located at the Cassegrain focus of the Keck I telescope. It has a 6.1'$\times$6.1' field of view and a plate scale of 0.18'', which is a bit coarse for sampling a 0.3''-0.4'' FWHM PSF from a GLAO feed. Ideally, MOSFIRE could be upgraded with a new detector (i.e. H4RG) with more and smaller pixels to preserve the field of view; but improve the sampling. The image quality inside MOSFIRE is excellent -- typically smaller than a single pixel \cite{McLean:2012}.  However, MOSFIRE lacks an atmospheric dispersion corrector, similar to DEIMOS, and would likely need one added to take full advantage of a GLAO feed.

\section{Conclusions}
We study the feasibility of deploying a GLAO system at the W.~M.~Keck Observatory. A GLAO system at Keck would likely deliver images with a 2$\times$ better FWHM, resulting in improved point source sensitivity, better image resolution, and possibly higher spectral resolution through the use of narrower slits on GLAO-fed spectrographs. While Keck's existing multi-object spectrographs (LRIS, DEIMOS, MOSFIRE) have sufficient instrumental image quality to take advantage of a GLAO feed, they each have a shortcoming that requires further investigation before we can determine the final delivered image quality. First, LRIS suffers from flexure as it is a Cassegrain mounted instrument. The degree of flexure may be large enough to prevent long exposures with a narrow (0.3'') slit. Second, DEIMOS and MOSFIRE both lack an atmospheric dispersion compensator. Again, narrowing the slit without an ADC could cause objects to shift out of the slit at certain wavelengths at low elevations due to atmospheric refraction. Installation of an ADC for DEIMOS and MOSFIRE will be explored. Investigations are also underway to determine the feasibility of implementing laser guide star wavefront sensors in front of each of the instruments, and of implementing a 1.4 m diameter adaptive secondary mirror at Keck.

\acknowledgments
We acknowledge funding from the W.~M.~Keck Observatory and the University of California Observatories to support the development of a GLAO concept for the Keck telescopes.