\section{Introduction}

Let $(x_1, x_2, \ldots)$ be a positive sequence with finite sum $t = \sum_{i=1}^\infty x_i$. Its size-biased permutation (s.b.p) is the same sequence presented in a random order $(x_{\sigma_1}, x_{\sigma_2}, \ldots)$, where $\mathbb{P}(\sigma_1 = i) = \displaystyle \frac{x_i}{t}$, and for $k$ distinct indices $i_1, \ldots, i_k$, 
\begin{equation}\mathbb{P}(\sigma_k = i_k| \sigma_1 = i_1, \ldots, \sigma_{k-1} = i_{k-1}) = \frac{x_{i_k}}{t - (x_{i_1} + \ldots + x_{i_{k-1}})}. \label{eqn:sbp} \end{equation}
An index $i$ with bigger `size' $x_i$ tends to appear earlier in the permutation, hence the name size-biased.
Size-biased permutation of a random sequence is defined by conditioning on the sequence values.

\subsection*{A brief history}
One of the earliest occurences of size-biased permutation is in social choice theory. Luce \cite{luce} studied distributions on the permutation on $n$ letters given by (\ref{eqn:sbp}) as a function of the $x_i$'s. These are the relative scores or desirabilities of the candidates, to be inferred through observing multiple rankings. This ranking model is now known as the Plackett-Luce model and it has wide applications \cite{plackett, plmeier, volkmer}.

Around the same time, biologists in population genetics were interested in inferring the distribution of alleles in a population through sampling. In these applications, $x_i$ is the abundance and $x_i/t$ is the relative abundance of the $i$-th species \cite{en78}. Size-biased permutation models the outcome of successive sampling, where one samples without replacement from the population and records the abundance of newly discovered species in the order that they appear. To account for the occurrence of new types of alleles through mutation and migration, they considered random abundance sequences and did not assume an upper limit to the number of possible types. Species sampling from random infinite sequence is sometimes known as \emph{size-biased random permutation}, a term coined by Patil and Taillie \cite{pt77}. The earliest work along this vein is perhaps that of McCloskey \cite{mc65}, who obtained results on the size-biased permutation of ranked jumps in a certain Poisson point process (p.p.p). The distribution of this ranked sequence is now known as the Poisson-Dirichlet $PD(0, \theta)$, and the distribution of its size-biased permutation is the $GEM(0, \theta)$; see Section \ref{subsec:parallel} for their definitions. This work was later generalized by Perman, Pitman and Yor \cite{pitmanyor}, who studied size-biased permutation of ranked jumps of a subordinator; see Section \ref{sec:sbp}.

The finite combinatorial version of size-biased sampling from $PD(0, \theta)$ is the Ewens sampling formula \cite{ew72}. Kingman \cite{ki78a} wrote: `One of the most striking results of recent theoretical research in population genetics is the sampling formula enunciated by Ewens and shown by (others) to hold for a number of different population models'. In studying this formula, Kingman initiated the theory of partition structure \cite{ki78a, ki78b}. Kingman showed that the Ewens sampling formula defines a particular partition structure by deletion of type; see \cite{ghp} for recent developments. Subsequent authors have studied partition structures and their representations in terms of exchangeable random partitions, random discrete distributions, random trees and associated random processes of fragmentation and coalescence, Bayesian statistics and machine learning. See \cite{csp} and references therein.     
\subsection*{Organization}
This paper focuses on \emph{finite i.i.d size-biased permutation}, that is, the size-biased permutation of $n$ independent and identically distributed (i.i.d) random variables from some distribution $F$ on $(0, \infty)$. Our setting is a finite dimensional analogue of the size-biased permutation of ranked jumps of a subordinator studied in Perman-Pitman-Yor (PPY) \cite{pitmanyor}, as well as a special form of \emph{induced order statistics} \cite{david73, bhat74}; see Section \ref{sec:history} for a brief historical account of this field. We utilize this connection in our paper to arrive at new results. In Section \ref{sec:sbp} we study joint and marginal distribution of finite i.i.d size-biased permutation through a Markov chain. We draw connections between our settings and that of PPY \cite{pitmanyor}, and prove a converse of the stick-breaking representation of Patil and Taillie when $F$ is gamma (Proposition \ref{prop:pt}). Section \ref{sec:ios} we show that finite i.i.d size-biased permutation is a form of induced ordered statistics, and use this fact to derive previous distributional results. Comparisons to corresponding statements in Section \ref{sec:sbp} lead to a new beta-gamma identity when F is gamma (Corollary \ref{cor:newgb2}). As the sequence length tends to infinity, we derive asymptotics of the last $u$ fraction of finite i.i.d size-biased permutation in Section \ref{sec:asymp}, and that of the first few terms in Section \ref{sec:ppp}.

\subsection*{Notations} We shall write $gamma(a,\lambda)$ for a Gamma distribution whose density at $x$ is $\lambda^a x^{a-1}e^{-\lambda x}/\Gamma(a)$ for $x > 0$, and $beta(a,b)$ for the Beta distribution whose density at $x$ is $\displaystyle \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}x^{a-1}(1-x)^{b-1}$ for $x \in (0,1)$. For an ordered sequence $Y_n(k), k = 1, \ldots, n$, let $Y^{rev}_n(k) = Y_n(n-k+1)$ be the same sequence presented in reverse. For order statistics, we write $(Y^{\uparrow}(1), Y^{\uparrow}(2), \ldots)$ for the increasing sequence, and $(Y^{\downarrow}(1), Y^{\downarrow}(2), \ldots)$ for the decreasing sequence.  Throughout this paper we use $(X_1, \ldots, X_n)$ for the underlying i.i.d sequence, and $(X_n[1], \ldots, X_n[n])$ for its size-biased permutation.

\section{Connections to the work of Perman-Pitman-Yor\cite{pitmanyor}.}\label{sec:sbp}
PPY considered the size-biased permutation of the sequence of ranked jumps of a \emph{subordinator} $\tilde{T}$, that is, a non-decreasing process with no drift component, right continuous paths, stationary independent increments, $\tilde{T}(0) = 0$, and
$$ \tilde{T}(s) = \sum_i X_i \mathbf{1}_{\{\sigma_i \leq s\}} $$ 
for $ 0 \leq s \leq 1$, where $\{(\sigma_i, X_i), i \geq 1\}$ is the countable random set of points of a Poisson point process (p.p.p) on $(0, \infty)^2$ with intensity measure $ds \Lambda (dx)$. The $X_i$'s are the jumps of the subordinator $\tilde{T}$. 
They assumed $\Lambda(dx) = \rho(x)\, dx$ for some density $\rho$, $\Lambda(0, \infty) = \infty, \Lambda(1, \infty) < \infty$, $\int_0^1 x \Lambda(dx) < \infty$. So $T > 0$ is an infinitely divisible random variable with L\'{e}vy measure $\Lambda$ and no drift component in its L\'{e}vy-Khitchine representation (see, for example, \cite[\S 15]{kallenberg}). 

Proposition \ref{prop:converge} states that the PPY setup is the limit in distribution of finite i.i.d s.b.p. This is an application of the principle that convergence of order statistics is equivalent to convergence of size-biased permutation; see \cite{gnedin} for a precise statement.
\begin{prop}[Gnedin \cite{gnedin}]\label{prop:converge}
Let $X = (X^{\downarrow}(1), X^{\downarrow}(2), \ldots)$ be the infinite sequence of ranked jumps of a subordinator, whose sum $T = \sum_{i\geq 1}X^\downarrow(i)$ is a positive, infinitely divisible, a.s. finite random variable whose L\'evy-Khitchine representation has no drift component. Let $(X^n, n \geq 1)$ be an i.i.d positive triangular array, that is, $X^n = (X^n_1, \ldots, X^n_n)$ where $X^n_i$ are i.i.d, a.s. positive, and $T_n = \sum_{i=1}^nX^n_i \stackrel{d}{\rightarrow} T$ as $n \to \infty$. Then as $n \to \infty$,
$$ \mbox{s.b.p of } X^n \,\, \stackrel{f.d.d}{\longrightarrow} \mbox{ s.b.p of } X $$
meaning convergence in distributoin of finite dimensional distributions.
\end{prop}
\begin{proof} Recall that the sequence of decreasing order statistics $(X^{n, \downarrow}(1),\ldots, X^{n, \downarrow}(n))$ converges in distribution to $X$ \cite[\S 15]{kallenberg}. Since $T_n, T > 0$ a.s. and $T_n \stackrel{d}{\rightarrow} T$, the normalized sequence $(X^{n, \downarrow}(1)/T, \ldots, X^{n, \downarrow}(n)/T)$ converges in distribution to $X/T$. The result follows from
\cite[Theorem 3]{gnedin}. 
\end{proof}

One can obtain another finite version of the PPY setup by setting $\Lambda(0, \infty) <~ \infty$, but this can be reduced to finite i.i.d s.b.p by conditioning. Specifically, $\tilde{T}$ is now a compound Poisson process, where the subordinator waits for an exponential time with rate $\Lambda(0, \infty)$ before making a jump $X$, whose length is independent of the waiting time and distributed as $\mathbb{P}(X \leq t) = \Lambda(0, t]/\Lambda(0, \infty)$ \cite{bertoin99}. If $(X_1, X_2, \ldots)$ is the sequence of successive jumps of $(\tilde{T}_s, s \geq 0)$, then $(X_1, X_2, \ldots, X_N)$ is the sequence of successive jumps of $(\tilde{T}_s, 0 \leq s \leq 1)$, where $N$ is a Poisson random variable with mean $\Lambda(0, \infty)$, independent of the jump sequence $(X_1, X_2, \ldots)$. For $N > 0$, properties of the size-biased permutation of $(X_1, \ldots, X_N)$ can be deduced from those of a finite i.i.d size-biased permutation by conditioning on $N$. 

\subsection{Joint distribution, Markov property and stick-breaking}\label{subsec:joint.ppy}
Recall that $(X_1, \ldots, X_n)$ is an i.i.d sequence from distribution $F$, and $(X_n[1], \ldots, X_n[n])$ is its size-biased permutation. Let $T_{n-k} = X_n[k+1] + \ldots + X_n[n]$ denote the sum of the last $n-k$ terms in an i.i.d size-biased permutation of length $n$. We now derive the joint distribution of the first $k$ terms $X_n[1], \ldots, X_n[k]$ when $F$ has density $\nu_1$. The parallels in PPY are discussed in Section ~\ref{subsec:parallel}. 

\begin{prop}[Barouch-Kaufman \cite{bkauf}] For $1 \leq k \leq n$, let $\nu_k$ be the density of $S_k$, the sum of $k$ i.i.d random variables with distribution $F$. Then
\begin{align}
& \mathbb{P}(X_n[1] \in dx_1, \ldots, X_n[k] \in dx_k) \nonumber \\
=& \frac{n!}{(n-k)!}\left(\prod_{j=1}^kx_j\, \nu_1(x_j)\, dx_j\right)\, \int_0^\infty \nu_{n-k}(s) \prod_{j=1}^k(x_j + \ldots + x_k+s)^{-1}\, ds. \label{eqn:joint.xnk}
\end{align}
\end{prop}
\begin{proof} 
Let $\sigma$ denote the random permutation on $n$ letters defined by size-biased permutation as in (\ref{eqn:sbp}). Then there are $\frac{n!}{(n-k)!}$ distinct possible values for $(\sigma_1, \ldots, \sigma_k)$. By exchangability of the underlying i.i.d random variables $X_1, \ldots, X_n$, it is sufficient to consider $\sigma_1 = 1$, $\ldots$, $\sigma_k = ~k$. Note that
$$\mathbb{P}\left((X_1, \ldots, X_k) \in dx_1\, \ldots dx_k, \sum_{j=k+1}^nX_j \in ds\right) = \nu_{n-k}(s) \, ds \prod_{j=1}^k \nu_1(x_j) \, dx_j.$$
Thus, restricted to $\sigma_1 = 1$, $\ldots$, $\sigma_k = ~k$, the probability of observing $(X_n[1], \ldots, X_n[k]) \in dx_1\, \ldots dx_k$ and $T_{n-k} \in ds$ is precisely
$$\frac{x_1\, dx_1}{x_1 + \ldots + x_k + s}\, \frac{x_2\, dx_2}{x_2 + \ldots + x_k + s} \cdots \frac{x_k\, dx_k}{x_k + s}\, \nu_{n-k}(s)\, \left(\prod_{j=1}^k\nu_1(x_j)\, dx_j\right) ds. $$
By summing over $\frac{n!}{(n-k)!}$ possible values for $(\sigma_1, \ldots, \sigma_k)$, and integrating out the sum $T_{n-k}$, we arrive at (\ref{eqn:joint.xnk}).
\end{proof}
Note that $X_n[k] = T_{n-k+1} - T_{n-k}$ for $k = 1, \ldots, n-1$. Thus we can rewrite (\ref{eqn:joint.xnk}) in terms of the joint law of $(T_{n}, T_{n-1}, \ldots, T_{n-k})$:
\begin{equation}\label{eqn:joint.T}
\mathbb{P}(T_n \in dt_0, \ldots, T_{n-k} \in dt_k) = \frac{n!}{(n-k)!}\left(\prod_{i=0}^{k-1} \frac{t_i - t_{i+1}}{t_i}\nu_1(t_i - t_{i+1})\right)\, \nu_{n-k}(t_k)\, dt_0\ldots dt_k.
\end{equation}
Rearranging (\ref{eqn:joint.T}) yields the following result, which appeared as an exercise in \cite[\S 2.3]{cyor}.
\begin{cor}[Chaumont-Yor \cite{cyor}]
The sequence $(T_n, T_{n-1},\ldots, T_1)$ is an inhomogeneous Markov chain with transition probability
\begin{equation}\label{eqn:T.tp}
 \mathbb{P}(T_{n-k} \in ds| T_{n-k+1} = t) = (n-k+1)\, \frac{t-s}{t}\cdot\nu_1(t-s)\, \frac{\nu_{n-k}(s)}{\nu_{n-k+1}(t)}\, ds,
\end{equation}
for $k = 1, \ldots, n-1$. Together with $T_n \stackrel{d}{=} S_n$, equation (\ref{eqn:T.tp}) specifies the joint law in (\ref{eqn:joint.T}), and vice versa. 
\end{cor}

An equivalent way to state (\ref{eqn:T.tp}) is that for $k \geq 1$, conditioned on $T_{n-k+1} = t$, $X_n[k]$ is distributed as the first size-biased pick out of $n-k+1$ i.i.d random variables conditioned to have sum $S_{n-k+1} = t$. This provides a recursive way to generate a finite i.i.d s.b.p: first generate $T_n$ (which is distributed as $S_n$). Conditioned on the value of $T_n$, generate $T_{n-1}$ via (\ref{eqn:T.tp}), let $X_n[1]$ be the difference. Now conditioned on the value of $T_{n-1}$, generate $T_{n-2}$ via (\ref{eqn:T.tp}), let $X_n[2]$ be the difference, and so on. Let us explore this recursion from a different angle by considering the ratio $\displaystyle W_{n,k} := \frac{X_n[k]}{T_{n-k+1}}$ and its complement, $\overline{W}_{n,k} = 1 - W_{n,k} = \displaystyle \frac{T_{n-k}}{T_{n-k+1}}$. For $k \geq 2$, note that
\begin{equation}\label{eqn:stickbreak}
 \frac{X_n[k]}{T_n} = \frac{X_n[k]}{T_{n-k+1}}\, \frac{T_{n-k+1}}{T_{n-k+2}}\cdots \frac{T_{n-1}}{T_n} = W_{n,k}\prod_{i = 1}^{k-1}\overline{W}_{n,i}.
\end{equation}
The variables $\overline{W}_{n,i}$ can be interpreted as residual fractions in a \emph{stick-breaking} scheme: start with a stick of length 1. Choose a point on the stick according to distribution $W_{n,1}$, `break' the stick into two pieces, discard the piece of length $W_{n,1}$ and rescale the remaining half to have length 1. Repeating this procedure $k$ times, and (\ref{eqn:stickbreak}) is the fraction broken off at step $k$ relative to the original stick length. 

Together with $T_n \stackrel{d}{=} S_n$, one could use (\ref{eqn:stickbreak}) to compute the marginal distribution for $X_n[k]$ in terms of the ratios $\overline{W}_{n,i}$. In general the $W_{n,i}$ are not necessarily independent, and their joint distributions need to be worked out from (\ref{eqn:T.tp}). However, when $F$ has gamma distribution, $T_n, W_{n,1}, \ldots, W_{n,k}$ are independent, and (\ref{eqn:stickbreak}) leads to the following result of Patil and Taillie \cite{pt77}.

\begin{prop}[Patil-Taillie \cite{pt77}]\label{prop:pt} If $F$ has distribution $gamma(a, \lambda)$ for some $a, \lambda > 0$, then $T_n$ and the $W_{n,1}, \ldots, W_{n,n-1}$ in (\ref{eqn:stickbreak}) are mutually independent. In this case,
\begin{align*}
X_n[1] &= \gamma_0 \beta_1 \\
X_n[2] &= \gamma_0 \bar{\beta}_1 \beta_2 \\
& \ldots \\
X_n[n-1] &= \gamma_0 \bar{\beta}_1 \bar\beta_2\ldots \bar\beta_{n-2}\beta_{n-1} \\
X_n[n] &= \gamma_0 \bar{\beta}_1 \bar\beta_2\ldots \bar\beta_{n-1}
\end{align*}
where $\gamma_0$ has distribution $gamma(an,\lambda)$, $\beta_k$ has distribution $beta(a+1, (n-k)a)$, $\bar\beta_k = 1 - \beta_k$ for $1 \leq k \leq n-1$, and the random variables $\gamma_0, \beta_1, \ldots, \beta_n$ are independent. 
\end{prop}
\begin{proof} By assumption $S_k$ has distribution $gamma(ak, \lambda)$. One substitutes the density of $gamma(ak, \lambda)$ for $\nu_k$ in (\ref{eqn:T.tp}), and the result follows by direct computation.
\end{proof}

In the subordinator setting, McCloskey \cite{mc65} proved the analogue of Proposition \ref{prop:pt}, and Perman-Pitman-Yor \cite{pitmanyor} proved the converse, see Proposition \ref{prop:mccloskey} below. Using the same idea, we obtain the following converse to Proposition \ref{prop:pt}, which appears to be new.
\begin{cor}\label{cor:pt} If $T_n$ is independent of $X_n[1]/T_n$, then $F$ is $gamma(a, \lambda)$ for some $a, \lambda > 0$. 
\end{cor}
\begin{proof} Lukacs \cite{lukacs} proved that if $X, Y$ are non-degenerate, positive independent random variables, then $X+Y$ is independent of $\frac{X}{X+Y}$ if and only if both $X$ and $Y$ have gamma distributions with the same scale parameter. Note that
$$
\mathbb{P}(X_n[1]/T_n \in du, T_n \in dt) = n\, u\, \mathbb{P}\left(\frac{X_1}{X_1 + (X_2+\ldots+X_n)} \in du, T_n \in dt\right)
$$
and 
$$\mathbb{P}(X_n[1]/T_n \in du) = n\, u\, \mathbb{P}\left(\frac{X_1}{X_1 + (X_2+\ldots+X_n)} \in du\right).$$ 
 Since $X_n[1]/T_n$ and $T_n$ are independent, $\frac{X_1}{X_1 + (X_2+\ldots+X_n)}$ is independent of $T_n$. The conclusion follows by applying Lukacs' theorem to the pair $X_1$ and $(X_2 + \ldots + X_n)$. 
\end{proof}

\subsection{Parallels in PPY \cite{pitmanyor}} \label{subsec:parallel}
In the subordinator setting of PPY, let $\tilde{T}_k$ denote the remaining sum after removing the first $k$ terms of the size-biased permutation of the infinite sequence $(X^\downarrow(1), X^\downarrow(2), \ldots)$. The sequence $(\tilde{T}_0, \tilde{T}_1, \ldots)$ is a Markov chain with stationary transition probabilities \cite[equation 2.a]{pitmanyor}
$$
 \mathbb{P}(\tilde{T}_{1} \in dt_1| \tilde{T}_0 = t) = \frac{t-t_1}{t}\cdot\rho(t-t_1)\, \frac{\nu(t_1)}{\nu(t)}\, dt_1.
 $$

Conditionally given $\tilde{T}_0 = t_0, \tilde{T}_1 = t_1, \ldots, \tilde{T}_n = t_n$, the sequence of remaining terms in the s.b.p $(X[n+1], X[n+2], \ldots)$ is distributed as $(X^\downarrow(1), X^\downarrow(2), \ldots, )$ conditioned on $\sum_{i\geq 1}X^\downarrow(i) = t_n$, independent of the first $n$ size-biased picks \cite[Theorem 4.2]{pitmanyor}. The stick-breaking representation in (\ref{eqn:stickbreak}) now takes the form
\begin{equation}\label{eqn:stick2}
 \frac{X[k]}{\tilde{T}_0} = W_{k}\prod_{i=1}^{k-1}\overline{W}_{i},
\end{equation}
where $X[k]$ is the $k$-th size-biased pick, and $W_{i} = \frac{X[i]}{\tilde{T}_{i-1}}$, $\overline{W}_{i} = 1 - W_{i} = \frac{\tilde{T}_i}{\tilde{T}_{i-1}}$. Proposition ~\ref{prop:pt} and Corollary \ref{cor:pt} parallel the following result.
\begin{prop}[McCloskey \cite{mc65} and PPY \cite{pitmanyor}]\label{prop:mccloskey}
The random variables $\tilde{T}_0$ and $W_{1}, W_2, \ldots$ in (\ref{eqn:stick2}) are mutually independent if and only if $\tilde{T}_0$ has distribution $gamma(a, \lambda)$ for some $a, \lambda > 0$. In this case, the $W_{i}$ are i.i.d with distribution $beta(1, a)$ for $i = 1, 2, \ldots$.
\end{prop}

We take a small detour to explain some results related to Proposition \ref{prop:pt} and \ref{prop:mccloskey} 
on characterization of size-biased permutations. For a random discrete distribution prescribed by its probability mass function $(P_k) = (P_1, P_2, \ldots)$ with $\sum_iP_i = 1$, $P_i \geq 0$, let $(\tilde{P}_k)$ be the s.b.p. of $(P_k)$. One may ask when is a given distribution $(Q_k)$ the s.b.p of some distribution $(P_k)$. Clearly $(\widetilde{\tilde{P}}_k) \stackrel{d}{=} (\tilde{P}_k)$ for any $(P_k)$, thus this question is equivalent to characterizing random discrete distributions on $\mathbb{N}$ which are invariant under size-biased permutation (ISBP). 
Pitman \cite[Theorem 4]{pitman96} gave a complete answer in terms of symmetry of a certain function of the finite dimensional distributions. Furthermore, Pitman proved a complete characterization of ISBP when $P_k$ can be written as the right hand side of (\ref{eqn:stick2}) with $W_{1}, W_2, \ldots$ independent. In this situation, apart from some limiting cases, $(P_k)$ is ISBP if and only if $W_{i}$ is distributed as $beta(1-\alpha, \theta + i\alpha)$, $i = 1, 2, \ldots$ for certain pairs of real numbers $(\alpha, \theta)$. The two main cases are $(0 \leq \alpha < 1, \theta > -\alpha)$ and $(\alpha = -a < 0, \theta = na)$ for some $n = 1, 2, \ldots$. In both settings, $(P_k)$ is known as the $GEM(\alpha, \theta)$ distribution. The abbreviation GEM was introduced by Ewens, which stands for Griffiths-Engen-McCloskey. The McCloskey case of Proposition \ref{prop:mccloskey} is $GEM(0, \theta)$, and the Patil-Taillie case of Proposition \ref{prop:pt} is $GEM(-a, na)$. The sequence obtained by ranking a $GEM(\alpha, \theta)$ sequence is called a Poisson-Dirichlet distribution with parameters $(\alpha, \theta)$ \cite{pitmanyor}. The $GEM$ distribution has a generative description known as the Chinese restaurant process \cite[\S 3]{csp} and has applications in Bayesian statistics and machine learning, see \cite{csp}.

\section{Connections to induced order statistics}\label{sec:ios}
When $n$ i.i.d pairs $(X_i, Y_i)$ are ordered by their $Y$-values, the corresponding $X_i$'s are called the \emph{induced order statistics} of the vector $Y$, or its \emph{concomitants}. Gordon \cite{gordon83} first proved the following result for finite $n$ which shows that finite i.i.d s.b.p is a form of induced order statistics. Here we state the infinite sequence version, which is a special case of PPY\cite[Lemma 4.4]{pitmanyor}.

\begin{prop}[PPY \cite{pitmanyor}]\label{prop:coupling} Let $(x_1, x_2, \ldots)$ be a fixed positive sequence with finite sum $t = \sum_{i=1}^\infty x_i$. For $i = 1, 2,\ldots$, let $Y_i = \epsilon_i/x_i$ where $\epsilon_i$ are i.i.d standard exponentials. Let $(Y^{\uparrow}(1),Y^{\uparrow}(2), \ldots)$ be the increasing order statistics of the $Y_i$'s, and let $X^\ast(k)$ be the value of the $x_i$ such that $Y_i$ is $Y^{\uparrow}(k)$. Then $(X^\ast(k), k = 1, 2, \ldots)$ is a size-biased permutation of $(x_1, x_2, \ldots)$. In particular, the size-biased permutation of a positive i.i.d sequence $(X_1, \ldots, X_n)$ is distributed as the induced order statistics of the sequence $(Y_i = \epsilon_i/{X_i}, 1 \leq i \leq n)$ for an independent sequence of i.i.d standard exponentials $(\epsilon_1, \ldots, \epsilon_n)$.
\end{prop}
\begin{proof} Note that the $Y_i$'s are independent exponentials with rates $x_i$. Let $\sigma$ be the random permutation such $Y_{\sigma(i)} = Y^\uparrow(i)$. Note that $X^\ast(k) = x_{\sigma(k)}$. Then
$$\mathbb{P}(\sigma(1) = i) = \mathbb{P}(Y_i = \min\{Y_j, j = 1, 2, \ldots\}) = \frac{x_i}{t},$$
thus $X_n^\ast(1) \stackrel{d}{=} x[1]$. In general, for distinct indices $i_1, \ldots, i_k$, by the memoryless property of the exponential distribution, 
\begin{align*}
& \, \mathbb{P}(\sigma(k) = i_k  | \sigma(1) = i_1, \ldots, \sigma(k) = i_{k-1} ) \\
=& \, \mathbb{P}\left(Y_{i_k} = \min\{Y_{\sigma(j)}, j \geq k \} | \sigma(1) = i_1, \ldots, \sigma(k) = i_{k-1}\right) = \frac{x_{i_k}}{t - \sum_{j=1}^{k-1} x_{i_j}}.
\end{align*}
Induction on $k$ completes the proof.
\end{proof}
 
This proposition readily supplies simple proofs for joint, marginal and asymptotic distributions of i.i.d s.b.p, as explicitly computed in this section. For instance, the proof of the following nesting property of i.i.d s.b.p, which can be cumbersome, amounts to i.i.d thinning.
\begin{cor}Consider a finite i.i.d s.b.p $(X_n[1], \ldots, X_n[n])$ from a distribution $F$. For $1 \leq m \leq n$, select $m$ integers $a_1 < \ldots < a_m$ by uniform sampling from $\{1, \ldots, n\}$ without replacement. Then the subsequence $\{X_n[a_j], 1 \leq j \leq m\}$ is jointly distributed as a finite i.i.d s.b.p of length ~$m$ from $F$.
\end{cor}  
 
\subsection{Joint and marginal distribution revisited} \label{subsec:ios.joint}

Proposition \ref{prop:coupling} immediately yields the joint distribution of the $X_n[k]$'s, which is a different formula to (\ref{eqn:joint.xnk}).
\begin{prop}\label{lem:joint.x.u}
$(X_n[k], k = 1, \ldots, n)$ is distributed as the first coordinate of the sequence of pairs $((X^{\ast}_n(k), U_n^\downarrow(k)), k =1, \ldots, n)$, where $U_n^\downarrow(1) \geq \ldots \geq U^\downarrow_n(n)$ is a sequence of uniform order statistics, and conditional on $(U^\downarrow_n(k) = u_k, 1 \leq k \leq n)$, the $X^{\ast}_n(k)$ are independent with distribution $(G_{u_k}(\cdot), k = 1, \ldots, n)$, where
\begin{equation} \label{eqn:Gu}
 G_u(dx) = \frac{xe^{-\phi^{-1}(u)x}\, F(dx)}{-\phi'(\phi^{-1}(u))}.
\end{equation}
Here $\phi$ is the Laplace transform of $X$, that is, $\phi(y) =~ \int_0^\infty e^{-yx}\, F(dx)$, $\phi'$ its derivative and $\phi^{-1}$ its inverse function.
\end{prop}
\begin{proof} Let $(X_1, \ldots, X_n)$ be $n$ i.i.d draws from $F$, $Y_i = \epsilon_i/X_i$ for $n$ i.i.d standard exponentials $\epsilon_i$, independent of the $X_i$'s. Note that $\{(X_i, Y_i), 1 \leq i \leq n\}$ is an i.i.d sample from the joint distribution $ F(dx)[xe^{-yx}\, dy].$ Thus $Y_i$ has marginal density
\begin{equation}
P(Y_i \in dy) = -\phi'(y)\, dy, \hspace{1em} 0 < y < \infty, \label{1.b}
\end{equation}
and its distribution function is $F_Y = 1 - \phi$.  Given $\{Y_i =  y_i, 1 \leq i \leq n\}$, the $X_n^\ast(i)$ defined in Proposition \ref{prop:coupling} are independent with conditional distribution $\widetilde{G}(y_i, \cdot)$ where
\begin{equation}
\widetilde{G}(y, dx) = \frac{xe^{-yx}\, F(dx)}{-\phi'(y)}. \label{eqn:Gy}
 \end{equation}

Equation (\ref{eqn:Gu}) follows from writing the order statistics as the inverse transforms of ordered uniform variables
\begin{equation}\label{eqn:ynun}
(Y_n^\uparrow(1), \ldots, Y^\uparrow_n(n)) \stackrel{d}{=} (F_Y^{-1}(U^\downarrow_n(n)), \ldots, F_Y^{-1}(U^\downarrow_n(1))) \stackrel{d}{=} (\phi^{-1}(U^\downarrow_n(1)), \ldots, \phi^{-1}(U^\downarrow_n(n)))
\end{equation}
where $(U^\downarrow_n(k), k = 1, \ldots, n)$ is an independent decreasing sequence of uniform order statistics. Note that the minus sign in (\ref{1.b}) results in the reversal of the sequence $U_n$ in the second equality of (\ref{eqn:ynun}).
\end{proof}
\begin{cor} \label{cor:marginal}
For $1 \leq k \leq n$ and $0 < u < 1$, let 
\begin{equation}\label{2.a}
f_{n,k}(u) = \frac{\mathbb{P}(U_{n}(k) \in du)}{du} = n\binom{n-1}{k-1}u^{n-k}(1-u)^{k-1}
\end{equation}
be the density of the $k$-th largest of the $n$ uniform order statistics $(U_n(i), i = 1, \ldots, n)$. Then
\begin{equation}\label{2.c}
\frac{\mathbb{P}(X_n[k] \in dx)}{xF(dx)} = \int_0^\infty e^{-xy}f_{n,k}(\phi(y))\, dy.
\end{equation}
\end{cor}
In particular, for the first and last values,
\begin{align*}
\frac{\mathbb{P}(X_n[1] \in dx)}{xF(dx)} &= n\int_0^\infty e^{-xy}\phi(y)^{n-1}\, dy \\
\frac{\mathbb{P}(X_n[n] \in dx)}{xF(dx)} &= n\int_0^\infty e^{-xy}(1-\phi(y))^{n-1}\, dy.
\end{align*}

\begin{ex}\label{ex:gamma1}
Suppose $F$ is $gamma(a,1)$. Then $\phi(y) = (\frac{1}{1+y})^a$, and $\phi^{-1}(u) =~ u^{-1/a} - 1.$ 
Hence $G_u$ in (\ref{eqn:Gu}) is
$$G_u(dx) = \frac{x}{au^{(a+1)/a}}\, e^{-(u^{-1/a} - 1)x}\, F(dx) = \frac{x^a}{\Gamma(a+1)}u^{-(a+1)/a}e^{-xu^{-1/a}}.$$
That is, $G_u$ is $gamma(a+1, u^{-1/a})$. 
\end{ex}

\subsection{A new beta-gamma identity}\label{subsec:ios.gamma}
When $F$ is $gamma(a,\lambda)$, Lemma \ref{lem:joint.x.u} gives the following result, which is a remarkable complement to the Patil-Taillie representation in Proposition \ref{prop:pt}. 

\begin{prop}\label{prop:newgb}
Suppose $F$ is $gamma(a,\lambda)$. Then $G_u$ is $gamma(a+1, \lambda u^{-1/a})$, and
\begin{equation}\label{eqn:xrev.g}
 (X_n[k], k = 1, \ldots, n) \stackrel{d}{=} ([U^\downarrow_n(k)]^{1/a}\gamma_k, k = 1, \ldots, n)
 \end{equation}
where $\gamma_1, \ldots, \gamma_n$ are i.i.d $gamma(a+1, \lambda)$ random variables, independent of the sequence of decreasing uniform order statistics $(U^\downarrow_n(1), \ldots, U^\downarrow_n(n))$. Alternatively, jointly for $k =~ 1, \ldots, n$
\begin{align*}
X_n^{rev}[1] &= \gamma_1 \beta_{an, 1} \\
X_n^{rev}[2] &= \gamma_2 \beta_{an, 1}\beta_{an-a, 1} \\
& \ldots \\
X_n^{rev}[n-1] &= \gamma_{n-1} \beta_{an, 1}\beta_{an-a, 1}\cdots \beta_{2a, 1} \\
X_n^{rev}[n] &= \gamma_{n} \beta_{an, 1}\beta_{an-a, 1}\cdots \beta_{a, 1}.
\end{align*}
where the $\beta_{an - ia,1}$ for $i = 0, \ldots, n-1$ are distributed as $beta(an-ia, 1)$, and they are independent of each other and the $\gamma_k$'s.
\end{prop}
\begin{proof} The distribution $G_u$ is computed in the same way as in Example \ref{ex:gamma1} and (\ref{eqn:xrev.g}) follows readily from Proposition \ref{lem:joint.x.u}.
\end{proof}
A direct comparison of the two different representations in Proposition \ref{prop:pt} and \ref{prop:newgb} creates $n$ distributional identities. For example, the equality $X_n[1] = X^{rev}_n[n]$ shows that the following two means of creating a product of independent random variables produce the same result in law:
\begin{equation}\label{eqn:beta.ini}
 \beta_{a+1, (n-1)a} \, \gamma_{an, \lambda} \stackrel{d}{=} \beta_{an, 1} \, \gamma_{a+1, \lambda}
 \end{equation}
 where $\gamma_{r,\lambda}$ and $\beta_{a,b}$ denote random variables with distributions $gamma(r, \lambda)$ and $beta(a,b)$, respectively. In fact, by comparing the consecutive ratios $\frac{X_n[i+1]}{X_n[i]}$ for $i = 1, \ldots, n-1$, the other $n-1$ identities all reduce to the following single equation after setting $ia - a = b$.
\begin{cor}\label{cor:newgb2} For $a, b, \lambda > 0$,
\begin{equation}\label{eqn:beta.ratio}
\frac{1 - \beta_{a+1,a+b}}{\beta_{a+1,a+b}}\beta_{a+1, b} \stackrel{d}{=} \frac{\gamma_{a+1,\lambda}}{\gamma'_{a+1,\lambda}}\beta_{a+b,1},
 \end{equation}
where the random variables $\beta_{a+1,a+b}, \beta_{a+1,b}, \beta_{a+b,1}, \gamma_{a+1,\lambda}, \gamma'_{a+1,\lambda}$ are mutually independent, beta and gamma distributed with parameters indicated by subscripts.
\end{cor}
The validity of this identity can be checked by comparing moments. Conversely, (\ref{eqn:beta.ini}) and (\ref{eqn:beta.ratio}) allow one to go between Proposition \ref{prop:newgb} and the Proposition \ref{prop:pt}.


\section{Asymptotics of the last $u$ fraction of the size-biased permutation}\label{sec:asymp}
In this section we derive Glivenko-Cantelli and Donsker-type theorems for the distribution of the last $u$ fraction of terms in a finite i.i.d s.b.p. These are special cases of more general statements which hold for arbitrary induced order statistics in $d$ dimensions (see Section \ref{sec:history}).
Features pertaining to i.i.d s.b.p are presented in the following lemma, which has an interesting successive sampling interpretation; see Section \ref{subsec:heuristic}.
\begin{lem}\label{lem:ode} Suppose $F$ has support on $[0, \infty)$ and finite mean. For $u \in (0,1)$, define 
\begin{equation}\label{eqn:Fu}
F_u(dx) = \frac{e^{-x\phi^{-1}(u)}}{u}\, F(dx)
\end{equation}
and extend the definition to $\{0, 1\}$ by continuity, where $\phi$ is the Laplace transform of $F$ as in Proposition \ref{lem:joint.x.u}. Then $F_u$ is a probability distribution on $[0, \infty)$ for all $u \in [0,1]$, and $G_u$ in (\ref{eqn:Gu}) satisfies
\begin{equation}\label{eqn:Gu.is.sb}
G_u(dx) = xF_u(dx)/ \mu_u
\end{equation}
where $\mu_u = \int x F_u(dx) = \frac{-\phi'(\phi^{-1}(u))}{u}$. Furthermore,
\begin{equation}\label{eqn:integral} \int_0^u G_s(dx)\, ds = F_u(dx) \end{equation}
for all $s \in [0,1]$. In other words, the density
$$f(u,x) = F_u(dx) / F(dx) = u^{-1}e^{-x\phi^{-1}(u)}$$
of $F_u$ with respect to $F$ solves the differential equation 
\begin{equation}\label{eqn:evo}
 \frac{d}{du}[uf(u,x)] = \frac{-xf(u,x)}{\mu_u} \end{equation}
with boundary condition $f(1,x) \equiv 1$. 
\end{lem}
\begin{proof} By direct computation. 
\end{proof}

We now state a Glivenko-Cantelli-type theorem which applies to size-biased permutation of finite deterministic sequences. Versions of this result are known in the literature \cite{bnw, holst, gordon83, rosen}, see discussions in Section \ref{sec:history}. We offer an alternative proof using induced order statistics.

\begin{thm}\label{thm:gc} Let $((x_i^n, 1 \leq i \leq n), 1 \leq n)$ be a deterministic triangular array of positive numbers  with corresponding c.d.f sequence $(E_n, 1 \leq n)$. Suppose
\begin{equation}
\sup_x|E_n(x) - F(x)| \to 0 \mbox{ as } n \to \infty \label{eqn:sup.condition}
\end{equation}
for some distribution $F$ on $(0, \infty)$. Let $u \in (0, 1]$. Let $E_{n,u}(\cdot)$ be the empirical distribution of the first $\lfloor nu \rfloor$ terms in a size-biased permutation of the sequence $(x_i^n, 1 \leq i \leq n)$. Then for each $\delta \in (0,1)$,
\begin{equation}\label{eqn:gc.main}
\sup_{u \in [\delta, 1]} \sup_\mathbf{I} | E_{n,u}(\mathbf{I}) - F_u(\mathbf{I})| \stackrel{a.s.}{\longrightarrow} 0 \mbox{ as } n \to \infty,
\end{equation}
where $\mathbf{I}$ ranges over all subintervals of $(0, \infty)$.
\end{thm}

\begin{proof} Define $Y^n_i = \epsilon_i/x^n_i$ for $i = 1, \ldots, n$ where $\epsilon_i$'s are i.i.d standard exponentials as in Proposition \ref{prop:coupling}. Let $H_n$ be the empirical distribution function (e.d.f) of the $Y^n_i$'s, 
$$H_n(y) := \frac{1}{n}\sum_{i=1}^n\mathbf{1}_{\{Y^n_i < y\}}.$$ 
Let $J_n$ denote the e.d.f of $(x^n_i, Y^n_i)$. By Proposition \ref{lem:joint.x.u},
\begin{equation}\label{eqn:eemp}
E_{n,u}(\mathbf{I}) = \frac{n}{\lfloor nu \rfloor}\, \frac{1}{n} \sum_{i=1}^n\mathbf{1}_{\{x^n_i \in \mathbf{I}\}}\mathbf{1}_{\{Y^n_i < H_n^{-1}(1-u)\}} = \frac{n}{\lfloor nu \rfloor}J_n(\mathbf{I} \times [0, H_n^{-1}(1-u)]).
\end{equation}
Fix $\delta \in (0, 1)$, and let $u \in [\delta, 1]$. Let $\phi$ be the Laplace transform of $F$ and $J$ the joint law of $(X, \epsilon/X)$, where $X$ is a random variable with distribution $F$, and $\epsilon$ is an independent standard exponential. Note that $\frac{1}{u}J(\mathbf{I} \times [0, \phi^{-1}(u)]) = F_u(\mathbf{I}).$
Thus
\begin{align}
 \mathbb{E}_{n, u}(\mathbf{I}) - F_u(\mathbf{I}) = & \left(\frac{n}{\lfloor nu \rfloor}J_n(\mathbf{I} \times [0, H_n^{-1}(1-u)]) - \frac{n}{\lfloor nu \rfloor}J_n(\mathbf{I} \times [0, \phi^{-1}(u)])\right) \nonumber \\
& + \left(\frac{n}{\lfloor nu \rfloor}J_n(\mathbf{I} \times [0, \phi^{-1}(u)]) - \frac{1}{u}J(\mathbf{I} \times [0, \phi^{-1}(u)])\right). \label{eqn:split}
\end{align}
Let us consider the second term. Note that
\begin{align*}
J_n(\mathbf{I} \times [0, \phi^{-1}(u)]) &= \int_0^\infty e^{-t\phi^{-1}(u)}\mathbf{1}_{\{t \in \mathbf{I}\}}\, E_n(dt).
\end{align*}
Since $E_n$ converges to $F$ uniformly and $e^{-t\phi^{-1}(u)}$ is bounded for all $t \in (0, \infty)$ and $u \in [\delta, 1]$, 
$$\sup_{u \in [\delta, 1]}\sup_\mathbf{I}\left|\frac{n}{\lfloor nu \rfloor}J_n(\mathbf{I} \times [0, \phi^{-1}(u)]) - \frac{1}{u}J(\mathbf{I} \times [0, \phi^{-1}(u)])\right| \stackrel{a.s.}{\longrightarrow} 0 \mbox{ as } n \to \infty.$$
Let us consider the first term.
Since $J_n$ is continuous in the second variable, it is sufficient to show that  
\begin{equation}\label{eqn:hncont}
 \sup_{u \in [\delta, 1]}|H_n^{-1}(1-u) - \phi^{-1}(u)| \stackrel{a.s.}{\longrightarrow} 0 \mbox{ as } n \to \infty.
 \end{equation}
 To achieve this, let $A_n$ denote the `average' measure 
$$A_n(y) := \frac{1}{n}\sum_{i=1}^n\mathbb{P}(Y^n_i < y) = 1 - \int_0^\infty e^{-xy}\, dE_n(x).$$ 
A theorem of Wellner \cite[Theorem 1]{wellner} states that if the sequence of measures $(A_n, n \geq 1)$ is tight, then the Prohorov distance between $H_n$ and $A_n$ converges a.s. to $0$ and $n \to \infty$. In this case, since $E_n$ converges to $F$ uniformly, $A_n$ converges uniformly to $1 - \phi$. Thus $H_n$ converges uniformly to $1-\phi$, and (\ref{eqn:hncont}) follows.
\end{proof}

In particular, when $E_n$ is the e.d.f of $n$ i.i.d picks from $F$, then (\ref{eqn:sup.condition}) is satisfied a.s. by the Glivenko-Cantelli theorem. Thus Theorem \ref{thm:gc} implies that $F_u$ is the limit of the empirical distribution function (e.d.f) of the last $u$ fraction in a finite i.i.d s.b.p.
\begin{cor}\label{cor:emp}For $u \in (0,1]$, let $F_{n,u}(\cdot)$ denote the empirical distribution of the last $\lfloor nu \rfloor$ values of an i.i.d s.b.p. with length $n$. For each $\delta \in (0,1)$, as $n \to \infty$,
\begin{equation}\label{emp.u}
\sup_{u \in [\delta, 1]} \sup_\mathbf{I} | F_{n,u}(\mathbf{I}) - F_u(\mathbf{I})| \stackrel{a.s.}{\longrightarrow} 0,
\end{equation}
where $\mathbf{I}$ ranges over all subintervals of $(0, \infty)$.
\end{cor}

Under the settings in Theorem \ref{thm:gc}, one can also derive a Donsker-type result via entropy methods described in \cite[\S 2]{vaart2}, complementing similar results in the literature \cite[\S 5]{bnw}. In particular, the first term of (\ref{eqn:split}) depends on the flunctuation of empirical uniform order statistics around its limit. Its scaling depends on the rate of convergence of $E_n$ to $F$. The second term of (\ref{eqn:split}) converges in distribution after rescaled by $\sqrt{n}$, by a version of Donsker's theorem, to a Brownian bridge. We omit the details. A similar decomposition applies in Theorem ~\ref{thm:fclt}.

\subsection{A heuristic interpretation} \label{subsec:heuristic}
Since $X_n^{rev}[\lfloor u \, n \rfloor]$ converges in distribution to $G_u$ for $u \in [0, 1]$, Corollary \ref{cor:emp} lends a sampling interpretation to Lemma \ref{lem:ode}. Equation (\ref{eqn:evo}) has the heuristic interpretation as characterizing the evolution of the mass at $x$ over time $u$ in a successive sampling scheme. To be specific, consider a successive sampling scheme on a large population of $N$ individuals, with species size distribution $H$. Scale time such that at time $u$, for $0 \leq u \leq 1$, there are $Nu$ individuals (from various species) remaining to be sampled. Let $H_u$ denote the distribution of species sizes at time $u$, and fix the bin $(x, x + dx)$ of width $dx$ on $(0, \infty)$. Then $NuH_u(dx)$ is the number of individuals whose species size lie in the range $(x, x + dx)$ at time $u$. Thus $\frac{d}{du}NuH_u(dx)$ is the rate of individuals to be sampled from this range of species size at time $u$. The probability of an individual whose species size is in $(x, x+dx)$ being sampled at time $u$ is $\frac{xH_u(dx)}{\int_0^\infty xH_u(dx)}$. As we scaled time such that $u \in [0,1]$, in time $du$ we sample $N \, du$ individuals. Thus
$$ \frac{d}{du}NuH_u(dx) = -N\, \frac{xH_u(dx)}{\int_0^\infty xH_u(dx) }.$$
Let $f(u,x) = H_u(dx) / H(dx)$, then as a function in $u$, the above equation reduces to (\ref{eqn:evo}). 


\subsection{Functional central limit theorem}\label{subsec:fclt}
We now state a functional central limit theorem for i.i.d s.b.p. This is a special case of the general result of Davydov and Egorov \cite{egorov89} for induced order statistics. 

\begin{thm}[Davydov and Egorov, \cite{egorov89}]\label{thm:fclt} Suppose the first two moments of $F$ are finite. For a distribution $H$, let $\mu(H)$, $\sigma(H)$ denote its mean and standard deviation. For $u \in (0, 1]$, define
\begin{align*}
\xi_n(u) &= \frac{1}{n}\sum_{j=1}^{\lfloor nu \rfloor} X^{rev}_n[j] \\
\eta_n(u) &= \frac{1}{n}\sum_{j=1}^n X_j \mathbf{1}(Y_j < \phi^{-1}(u)) \\
m(u) &= \int_0^u \mu(G_s) \, ds = \mu(F_u).
\end{align*}
Then $\xi_n \to m$, $\eta_n \to m$ uniformly on $[0,1]$, and
$$\sqrt{n}(\xi_n - m) \Rightarrow \alpha, \hspace{2em} \sqrt{n}\, (\eta_n - m) \Rightarrow \eta $$
in the Skorokhod topology, where 
$$\alpha(u) = \int_0^u \sigma(G_s) dW(s), \hspace{1em} \eta(u) = \alpha(u) + \int_0^u m(s) \, dV(s).$$ 
$W$ is a standard Brownian motion, and $V$ is a Brownian bridge, independent of $W$. 
\end{thm}
The difference in $\xi_n$ and $\eta_n$ is the fluctuation of empirical uniform order statistics around its limit. Recall the analogue discussed at the end of Theorem \ref{thm:gc}. The proof of Theorem \ref{thm:fclt} can be found in \cite{egorov89}, together with similar statements on functional law of iterated logarithm for the processes $\eta_n, \xi_n$. 
\subsection{Historical notes on induced order statistics and successive sampling} \label{sec:history}
Induced order statistics were first introduced by David \cite{david73} and independently by Bhattacharya \cite{bhat74}. Typical applications stem from modeling an indirect ranking procedure, where subjects are ranked based on their $Y$-attributes although the real interest lies in ranking their $X$-attributes, which are difficult to obtain at the moment where the ranking is required.\footnote{One often uses $X$ for the variable to be ordered, and $Y$ for the induced variable, with the idea that $Y$ is to be predicted. Here we use $X$ for the induced order statistics since $X_n[k]$ has been used for the size-biased permutation. The role of $X$ and $Y$ in our case is interchangeable, as evident when one writes $X_iY_i = \epsilon_i$.} For example in cattle selection, $Y$ may represent the genetic makeup, for which the cattle are selected for breeding, and $X$ represents the milk yields of their female offspring. Thus a portion of this literature focuses on comparing distribution of induced order statistics to that of usual order statistics \cite{dayang77, yang77, nagada94, hall94}. 
The most general statement on asymptotic distributions is obtained by Davydov and Egorov \cite{egorov89}, who proved functional central limit theorem and functional law of the iterated logarithm for the process $S_{n,u}$ under tight assumptions. Their theorem translates directly into Theorem \ref{thm:fclt} for finite i.i.d s.b.p as discussed in Section \ref{subsec:fclt}. 
Various versions of results in Section \ref{sec:asymp}, including Theorem \ref{thm:fclt} are also known in the successive sampling community \cite{holst, rosen, gordon83, bnw, sen}. For example, Bickel, Nair and Wang \cite{bnw} proved Theorem \ref{thm:gc} with convergence in probability when $E_n$ and $F$ have the same discrete support on finitely many values.

\section{Poisson coupling of size-biased permutation and order statistics}\label{sec:ppp}

Comparisons between the distribution of induced order statistics and order statistics of the same sequence have been studied in the literature \cite{dayang77, yang77, nagada94, hall94}. However, finite i.i.d s.b.p has the special feature that there exists an explicit coupling between these two sequences as described in Proposition ~\ref{prop:coupling}. Using this fact, we now derive Theorem \ref{thm}, which gives a Poisson coupling between the \emph{last} $k$ size-biased terms $X_n^{rev}[1], \ldots, X_n^{rev}[k]$ and the $k$ smallest order statistics $X_n^{\uparrow}(1), \ldots X^{\uparrow}_n(k)$ as $n \to \infty$. \footnote{The superscript `rev' indicates that the order statistics in consideration are arranged in increasing order, as opposed to decreasing, which has been the convention of this paper.}
The existence of a Poisson coupling should not be surprising, since it is well-known that the increasing sequence of order statistics $(X_n^{\uparrow}(1), X_n^{\uparrow}(2), \ldots)$ converges to points in a Poisson point process (p.p.p) whose intensity measure depends on the behavior of $F$ near the infimum of its support, which is $0$ in our case. This standard result in order statistics and extreme value theory dates back to R\'{e}nyi \cite{renyi}, and can be found in \cite{evtHaan}. 

Our theorem is closely related to a result in PPY \cite[\S 4]{pitmanyor}. When $(X^\downarrow(1), X^\downarrow(2), \ldots)$ are ranked jumps of a subordinator, these authors noted that one can couple the size-biased permutation with the order statistics via the following p.p.p
\begin{equation}\label{eqn:ppois}
 N(\cdot) := \sum_{k\geq 1}\mathbf{1}[(X[k],Y^{\uparrow}(k)) \in \cdot] = \sum_{k\geq1}\mathbf{1}[(X^{\downarrow}(k),Y_k) \in \cdot]. \end{equation}
 $N(\cdot)$ has measure $ m(dxdy) = xe^{-xy}\, \Lambda(dx) dy. $ The first expression in (\ref{eqn:ppois}) defines a scatter of $(x,y)$ values in the plane listed in increasing $y$ values, and the second represents the same scatter listed in decreasing $x$ values. Since $\sum_{i\geq 1}X^\downarrow(i)$ is a.s. finite, the $x$-marginal of the points in  (\ref{eqn:ppois}) has the distribution of the size-biased permutation of the infinite sequence $(X^\downarrow(1), X^\downarrow(2), \ldots)$ since it prescribes the joint distribution of the first $k$ terms $X[1], \ldots, X[k]$ for any finite $k$. 
PPY use this p.p.p representation to generalize size-biased permutation to $h$-biased permutation, where the `size' of a point $x$ is replaced by an arbitrary strictly positive function $h(x)$; see \cite[\S 4]{pitmanyor}.

More generally, one obtains a random permutation of $\mathbb{N}$ from ranking points $(x_i, y_i)$ in a Poisson scatter $N(\cdot)$ on $(0, \infty)^2$ according to either their $x$ or $y$ values. For $i = 1, 2, \ldots$, let $(x^\ast_i)$ and $(y^\ast_i)$ denote the induced order statistics of the sequence $(x_i)$ and $(y_i)$ obtained by ranking points by their $y$ and $x$ values in increasing order, respectively. For $j, k = 1, 2, \ldots$, define sequences of integers $(K_j), (J_k)$ such that $x^{\uparrow}(J_k) = x^{\ast}_k$, $y^{\uparrow}(K_j) = y^{\ast}_j$; see Figure \ref{fig:scatter}.

\begin{figure}[hbtp]
\vskip-1.5em
\includegraphics[width = 0.7\textwidth]{./pic1.pdf}
\vskip-4.5em
\caption{A point scatter on the plane. Here $J_1 = 5, J_2 = 3, J_3 = 1, J_4 = 4, J_5 = 2$, and $K_1 = 3, K_2 = 5, K_3 = 2, K_4 = 4, K_5 = 1$. The permutations $J$ and $K$ are inverses. Conditioned on $x^\uparrow(2) = a$ and $y_2^\ast = b$, the number of points lying in the shaded region determines $J_2 - 1$.}
\label{fig:scatter}
\end{figure}
Suppose $N(\cdot)$ has intensity measure $m$ such that for all $x, y \in (0, \infty)$, $m((0, x) \times (0, \infty)) < \infty, m((0, \infty) \times (0, y)) < \infty$. For $j \geq 1$, conditioned on $x(j) = x, y^{\ast}_j = y$, 
\begin{equation}\label{eqn:kj}
K_j - 1 \stackrel{d}{=} Poisson\left(m((x, \infty) \times (0, y))\right) + Binomial\left(j-1, \frac{m((0, x) \times (0, y))}{m((0, x) \times (0, \infty))}\right),
\end{equation}
where the two random variables involved are independent. Similarly, for $k \geq 1$, conditioned on $x^\ast_k = x$, $y(k ) = y$, 
\begin{equation}\label{eqn:jk}
J_k - 1 \stackrel{d}{=} Poisson\left(m((0,x) \times (y,\infty))\right) + Binomial\left(k-1, \frac{m((0, x) \times (0, y))}{m((0, \infty) \times (y, \infty))}\right),
\end{equation}
where the two random variables involved are independent. When $m$ is a product measure, as is the case of i.i.d s.b.p, it is possible to compute the marginal distribution of $K_j$ and $J_k$ explicitly for given $j, k \geq 1$. We demonstrate such computations in Proposition ~\ref{prop:compute.general}. 

Before stating the theorem we need some technical results. The distribution of the last few size-biased picks depends on the behavior of $F$ near $0$, the infimum of its support. We shall consider the case where $F$ has `power law' near $0$, like that of a Gamma distribution. 
\begin{lem} \label{lem:approx.G.by.Gamma}
Suppose $F$ is supported on $(0, \infty)$ with Laplace transform $\phi$. Let $u = \phi(y)$, $X_u$ a random variable distributed as $G_u(dx)$ defined in (\ref{eqn:Gu}). For $\lambda, a > 0$,
\begin{equation}
F(x) \sim \frac{\lambda^a x^a}{\Gamma(a + 1)} \mbox{ as } x \to 0, \label{eqn:assump.F}
\end{equation}
if and only if,
\begin{align}
\phi(y) & \sim \lambda^a/y^{a} \hspace{1em} \mbox{ as } y \to \infty \label{eqn:phi.asymp}.
\end{align}
Furthermore, (\ref{eqn:assump.F}) implies
\begin{align}
u^{-1/a}X_u & \stackrel{d}{\longrightarrow}  gamma(a+1,\lambda) \hspace{1em} \mbox{ as } u \to 0.
\label{eqn:xu.gamma}
\end{align}

\end{lem}
\begin{proof} The equivalence of (\ref{eqn:assump.F}) and (\ref{eqn:phi.asymp}) follows from a version of Karamata Tauberian Theorem \cite[\S 1.7]{regvar}. Assume (\ref{eqn:assump.F}) and (\ref{eqn:phi.asymp}). We shall prove (\ref{eqn:xu.gamma}) by looking at the Laplace transform of the non-size-biased version $X_u'$, which has distribution $F_u$. For $\theta \geq 0$, 

\begin{equation}\label{eqn:laplace.xpu}
 E(\exp(-\theta X_u')) = \int_0^\infty u^{-1}\exp(-yx - \theta x)\, F(dx) = u^{-1}\phi(y + \theta) = \frac{\phi(y+\theta)}{\phi(y)}.
 \end{equation}
Now as $y \to \infty$ and $u = \phi(y) \to 0$, for each fixed $\eta > 0$, (\ref{eqn:phi.asymp}) implies
$$
E(\exp(-\eta u^{-1/a}X_u')) = \frac{\phi(y + \eta\phi(y)^{-1/a})}{\phi(y)} 
\sim \frac{\lambda^a(y+\eta \lambda^{-1} y)^{-a} }{\lambda^ay^{-a}} = \left(\frac{\lambda}{\lambda+\eta}\right)^a.
$$
That is to say
\begin{equation}\label{eqn:xuprime.gamma}
u^{-1/a} X_u' \stackrel{d}{\longrightarrow} gamma(a, \lambda).
\end{equation}
Since $\phi$ is differentiable, (\ref{eqn:laplace.xpu}) implies $E(X'_u) = \phi'(y) / \phi(y)$. Now $\phi$ has an increasing derivative $\phi'$, thus (\ref{eqn:phi.asymp}) implies $\phi'(y) \sim a\lambda^a/y^{a+1}$ as $y \to \infty.$ Therefore,
$$u^{-1/a}E(X_u') = \frac{\phi'(y)}{\phi(y)^{1 + 1/a}} \to \frac{a}{\lambda},$$ 
which is the mean of a $gamma(a,\lambda)$ random variable. Thus the random variables $u^{-1/a}X'_u$ are uniformly integrable, so for any bounded continuous function $h$, we can compute
$$E(h(u^{-1/a}X_u)) = \frac{E[(u^{-1/a}X_u')h(u^{-1/a}X_u')]}{u^{-1/a}E(X_u')} \to \frac{E[\gamma_{a,\lambda}h(\gamma_{a,\lambda})]}{E(\gamma_{a,\lambda})} = E(h(\gamma_{a+1,\lambda})) $$
where $\gamma_{b,\lambda}$ is a $gamma(b,\lambda)$ random variable. This proves (\ref{eqn:xu.gamma}).
\end{proof}

Now suppose $F$ satisfies (\ref{eqn:assump.F}) for some $\lambda, a > 0$. By standard results in order statistics \cite[Theorem 2.1.1]{evtHaan}, as $n \to \infty$, 
\begin{equation}\label{eqn:fidi.order}
\left(n^{1/a} X^{rev}_n(k), 1 \leq k \leq n\right) \stackrel{f.d.d.}{\longrightarrow} \left(\xi^{rev}(k): 1 \leq k < \infty\right),
\end{equation}
where f.d.d. stands for finite dimensional distribution. Here
\begin{equation}\label{eqn:fidi.sbs}
\xi^{rev}(k) = S_{(k)}^{1/a}\Gamma(a+1)^{1/a}/\lambda
\end{equation}
for $S_{(k)} = \epsilon_1 + \ldots + \epsilon_k$, where $\epsilon_i$ are i.i.d standard exponentials. 
We now present the analogue of (\ref{eqn:fidi.sbs}) for the last few size-biased picks $X_n^{rev}[1], \ldots, X_n^{rev}[k]$ and the promised Poisson coupling.
\begin{thm} \label{thm}
Suppose that (\ref{eqn:assump.F}) holds for some $\lambda, a > 0$. \\
\textbf{a)} As $n \to \infty$, 
\begin{equation}\label{15.a}
\left(n^{1/a} X_n^{rev}[k], 1 \leq k \leq n\right) \stackrel{f.d.d}{\longrightarrow} \left(\xi^{rev}[k]: 1 \leq k < \infty\right),
\end{equation}
with 
\begin{equation} \label{15.b}
\xi^{rev}[k] = T_{(k)}^{1/a}\gamma_k / \lambda,
\end{equation}
 where $T_{(k)} = \epsilon'_1 + \ldots + \epsilon'_k$ for i.i.d standard exponential r.v. $\epsilon'_i$, $\gamma_k, k = 1, \ldots, n$ are i.i.d $gamma(a+1,1)$, and the $(\epsilon_k')$ and $(\gamma_k)$ are independent. 
\vskip12pt \noindent
\textbf{b)} Let $\xi^{rev}(1) < \xi^{rev}(2) < \ldots $ be the increasing order statistics of the $\xi^{rev}[k]$'s. Then the $\xi^{rev}(k)$'s are jointly distributed as (\ref{eqn:fidi.sbs}), for a sequence $(\epsilon_k)$ which is implicitly defined in terms of the $(\epsilon'_k)$ and $(\gamma_k)$ above, and the f.d.d. convergence in (\ref{eqn:fidi.order}) and (\ref{15.a}) hold jointly.
\vskip12pt \noindent
\textbf{c)} For each $n$, let $J_n = (J_{nk}, 1 \leq k \leq n)$ be the permutation of $\{1, \ldots, n\}$ defined by
$$ X^{rev}_n[k] = X_n(J_{nk}). $$
As $n \to \infty$,
\begin{equation}\label{15.g}
\left(J_{nk}, 1 \leq k \leq n\right) \stackrel{f.d.d.}{\longrightarrow} \left(J_k: 1 \leq k < \infty\right),
\end{equation}
where $J_k$ is the random permutation of $\{1, 2, \ldots, \}$, defined by
\begin{equation}\label{15.h}
 \xi^{rev}[k] = \xi^{rev}(J_k)
 \end{equation}
 for $k = 1, 2, \ldots,$ and the f.d.d. convergence in (\ref{eqn:fidi.order}), (\ref{15.a}), (\ref{15.g}) all hold jointly.
\end{thm}
\begin{proof} By Lemma \ref{lem:approx.G.by.Gamma}, it is sufficient to prove the theorem for the case $F$ is $gamma(a, \lambda)$. Part a follows immediately from Proposition \ref{prop:newgb} and law of large numbers, in the same way that part b follows from the proof of (\ref{eqn:fidi.sbs}). For the last part, note that in principle everything has been presented as a function of the variables $T_{(k)}$ and $\gamma_k$. For $a, \lambda > 0$, define
$$\Psi_{a,\lambda}(s) = s^{1/a}\Gamma(a+1)^{1/a}/\lambda.$$ 
Observe that $\Psi_{a,\lambda}$ has inverse function $\Psi_{a,\lambda}^{-1}(x) = \lambda^a x^a / \Gamma(a+1)$. Rewrite (\ref{eqn:fidi.sbs}) as
\begin{equation}\label{eqn:xik}
\xi^{rev}(k) = \Psi_{a,\lambda}(S_{(k)}),
\end{equation}
then (\ref{15.b}) implies
\begin{equation}\label{16.a}
 S_k = \Psi_{a,\lambda}^{-1}(\xi^{rev}[k]) = \lambda^a \xi^{rev}[k]^a/\Gamma(a+1) = T_{(k)}\gamma_k^a / \Gamma(a+1).
 \end{equation}
Comparing (\ref{eqn:xik}) and (\ref{16.a}) gives a pairing between $S_k$ and $S_{(k)}$ via $\xi(k)$ and $\xi^{rev}[k]$. Hence we obtain another definition of $J_k$ equivalent to (\ref{15.h}):
$$ S_k = S_{(J_k)}. $$
Let $T_j$ be the $T$ value corresponding to the order statistic $S_{(j)}$ of the sequence $S_k$. That is, 
$$T_j = T_{(K_j)}$$
where $(K_j)$ is a random permutation of the positive integers. By (\ref{16.a}), $(J_k)$ is the inverse of $(K_j)$. Together with (\ref{eqn:fidi.order}) and (\ref{15.a}), this implies (\ref{15.g}), completing the proof of part c.
\end{proof}
 Since the $T_{(k)}$'s are increasing ordered points of a p.p.p on $(0, \infty)$ with rate 1 and independent of the $\gamma_k$'s, the point process on the positive quadrant $(0, \infty)^2$
\begin{equation}\label{16.b}
 N(\cdot) = \sum_k\mathbf{1}[(S_k, T_{(k)}) \in \cdot] = \sum_j\mathbf{1}[(S_{(j)}, T_j) \in \cdot]
 \end{equation}
is a p.p.p for the measure $\mu$ with density
\begin{equation}\label{eqn:measure.mu}
 \frac{\mu(ds\, dt)}{ds\, dt} = \frac{1}{a} \Gamma(a+1)^{1/a}(s/t)^{1/a}\exp\{-(\Gamma(a+1)s/t)^{1/a}\} 
 \end{equation}
The random permutation $J_k$ and its inverse $K_j$ are obtained from the p.p.p in (\ref{16.b}) in the $(s,t)$ plane by ranking the points in two different orders according to $s$-value or $t$-value. Since the projection of a Poisson process is Poisson, the $s$ and $t$-marginal of $\mu$ is just Lebesgue measure. Thus, the ranked $s$-values form a p.p.p with rate 1 which determine the ranked $\xi(j)$ by the \emph{deterministic} increasing transformation 
$$ \xi^{rev}(k) = \Psi_{a,\lambda}(S_{(k)}).$$ 
Furthermore,
$$\xi^{rev}[k] = \xi^{rev}(J_k) = \Psi_{a,\lambda}(S_k)$$ 
comes from the same set of $s$-values listed in order of increasing $t$-values. 

Marginal distributions of the random permutation $(J_k)$ and its inverse $(K_j)$ are given in (\ref{eqn:kj}) and (\ref{eqn:jk}). Note that for $k = 1,2, \ldots$, $S_k = T_{(k)}\gamma_k^a/\Gamma(a+1)$ for i.i.d $\gamma_k$ distributed as $gamma(a+1,1)$, independent of the sequence $(T_{(k)})$, and $T_k = \Gamma(a+1)S_{(k)}\tilde{\epsilon}_k^{-a}$ for i.i.d standard exponentials $\epsilon_k$, independent of the sequence $(S_{(k)})$ but not of the $\gamma_k$'s. Thus by conditioning on either $S_{(k)}$ or $T_{(k)}$, one can evaluate (\ref{eqn:kj}) and (\ref{eqn:jk}) explicitly. In particular, by a change of variable $r = \Gamma(a+1)^{1/a}(s/t)^{1/a}$, one can write $\mu$ in product form. This leads to the following.
\begin{prop}\label{prop:compute.general}
For $j \geq 1$, conditioned on $S_{(j)} = s, T_j = \Gamma(a+1)\, s \, r^{-a}$ for some $r > 0$, $K_j-1$ is distributed as
\begin{equation}\label{eqn:marginal.kj}
Poisson(m(s,r)) + Binomial(j-1, p(s,r))
\end{equation}
with 
\begin{equation}\label{eqn:ppp.density.kj}
m(s,r) = asr^{-a}\, \int_r^\infty x^{a-1}e^{-x}\, dx,
\end{equation}
and
\begin{equation}\label{eqn:psr}
p(s,r) = as^{2/a-2}r^{a-2}\int_0^r x^{a-1}e^{-x}\, dx,
\end{equation}
where the Poisson and Binomial random variables are independent.
Similarly, for $k \geq 1$, conditioned on $T_{(k)} = t, S_k = t \, r^a / \Gamma(a+1)$ for some $r > 0$, $J_k - 1$ is distributed as
\begin{equation}\label{eqn:marginal.jk}
 Poisson(m'(t, r)) + Binomial(k-1, p'(t,r))
 \end{equation}
with
\begin{equation}\label{eqn:ppp.density.jk} 
m'(t,r) = t \, \left(\frac{r^a + a \, \int_r^\infty x^{a-1}e^{-x}\, dx}{\Gamma(a+1)} - 1\right),
\end{equation}
and 
\begin{equation}\label{eqn:ptr}
p'(t,r) = \Gamma(a+1)^{1-2/a}\, at^{2/a-2}r^{a-1/a}\int_0^r x^{a-1}e^{-x}\, dx,
\end{equation}
where the Poisson and Binomial random variables are independent.
\end{prop}
\begin{prop}[Marginal distributions of $K_1$ and $J_1$]\label{prop:compute} 
Suppose that (\ref{eqn:assump.F}) holds for some $\lambda > 0$ and $a = 1$. Then the distribution of $K_1$, the $k$ such that $\xi(1) = \xi^{rev}[k]$, is a mixture of geometric distributions, and so is that for $J_1$, the $j$ such that $\xi^{rev}[1] = \xi(j)$. In particular,
\begin{equation}\label{18.k1}
P(K_1 = k) = \int_0^\infty p_rq_r^{k-1}e^{-r}\, dr
\end{equation}
where $p_r = r/(r + e^{-r})$, $q_r = 1 - p_r$, and
\begin{equation}\label{18.j1}
P(J_1 = j) = \int_0^\infty \tilde{p}_r\tilde{q}_r^{j-1}re^{-r}\, dr
\end{equation}
where $\tilde{p}_r = 1/(r+e^{-r})$, $\tilde{q}_r = 1 - \tilde{p}_r$.
\end{prop}
\begin{proof} When $a = 1$, $\int_r^\infty t^{a-1}e^{-t}\, dt = e^{-r}$. Substitute to (\ref{eqn:ppp.density.kj}) and (\ref{eqn:ppp.density.jk}) give 
$$m(s,r) = sr^{-1}e^{-r}, \hspace{1em} m'(t,r) = t(r-1+e^{-r}).$$ 
By a change of variable, (\ref{eqn:measure.mu}) becomes 
$$\frac{\mu(ds\, dr)}{ds\, dr} = se^{-r}, \hspace{1em} \frac{\mu(dt\, dr)}{dt\, dr} = tre^{-r}.$$ 
Thus, \emph{conditioned on $s$ and $r$}, $K_1 - 1$ is distributed as the number of points in a p.p.p  with rate $r^{-1}e^{-r}$ before the first point in a p.p.p with rate 1. This is the geometric distributions on $(0, 1, \ldots)$ with parameter $p_r = 1/(1 + r^{-1}e^{-r})$. Since the marginal density of $r$ is $e^{-r}$, integrating out $r$ gives (\ref{18.k1}). The computation for the distribution of $J_1$ is similar.
\end{proof}
One can check that (\ref{18.k1}) and (\ref{18.j1}) sum to 1. We conclude with a `fun' computation. Suppose that (\ref{eqn:assump.F}) holds for some $\lambda > 0$ and $a = 1$. That is, $F$ behaves like an exponential c.d.f near $0$. By Proposition \ref{prop:compute}, $E(J_1) = 9/4$ and $E(K_1) = \infty$. That is, the last size-biased pick is expected to be almost the second smallest order statistic, while the smallest order statistic is expected to be picked infinitely earlier on in a successive sampling scheme(!). The probability that the last species to be picked in a successive sampling scheme is also the one of smallest species size is
\begin{align*}
\lim_{n \to \infty} P(X_n^{rev}[1] = X_n(1)) &= P(\xi^{rev}[1] = \xi(1)) = P(J_1 = 1) = P(K_1 = 1) \\
&= \int_0^\infty \frac{re^{-r}}{r + e^{-r}}\, dr = 1 - \int_0^1 \frac{u}{u - \log u}\, du \approx 0.555229
\end{align*}



\bibliographystyle{plain}

