# Sequential Unequal Probability Sampling For Stream Population

Bardia Panahbehagh111Department of Mathematics, Kharazmi University, Tehran,
Iran (E-mail<EMAIL_ADDRESS>, Raphaël Jauslin222Institute of
statistics, University of Neuchâtel, Bellevaux 51, 2000 Neuchâtel, Switzerland
and Yves Tillé22footnotemark: 2

###### Abstract

A new unequal probability sampling method is proposed. This method is
sequential. The decision to select or not each unit is made based on the order
in which the units appear. A variant of this method allows to select a sample
from a stream. At each step, the decision to take the units successively
according to the order of appearance in the stream is made. This method
involves using a sliding window that is as small as possible. The method also
allows the sample to be spread and even the level of spreading to be adjusted.

Keywords algorithms, inclusion probability, flow sampling, window

##  1 Introduction

There are dozens of methods of unequal probability sampling. Most of them are
described in the following books Hanif and Brewer, (1980); Gabler, (1990);
Tillé, (2006). Nevertheless, very few methods allow you to decide to select
the units in the order in which they appear. Systematic sampling with unequal
probabilities (Madow,, 1949) allows sequential sampling from a stream.
However, once the first units are examined, this method quickly becomes
deterministic, making it predictable, which can be problematic. Indeed,
sampling cannot be predictable if it is used to perform a control (Busnel and
Tillé,, 2020).

The ordered pivotal method (Deville and Tillé,, 1998) is an alternative to
systematic sampling. It is a special case where pairs of units are processed
according to their order in the data list. Pivotal method and Deville’s
systematic sampling (Deville,, 1998; Tillé,, 2006) are two different
implementations of the same sampling design (Chauvet,, 2012). A disadvantage
of the pivotal method is that the decision to take or not take the first unit
may be made after examining a very large number of units in the sampling
frame.

There are also so-called reservoir methods that allow sampling from a stream.
Reservoir’s size could be the sample size (Chao,, 1982). It is updated each
time a unit enters in the stream (see also Cohen et al.,, 2009). Reservoir
methods were generalized to select balanced samples (Tillé,, 2019). However,
all these methods suffer from the same problem: the decision about a unit that
is at the beginning of the stream can be made very late.

The proposed method allows us to make decisions about the selection of units
according to their order of appearance in a list. The update of the
probabilities is realized only on a window containing a very small number of
units which follows the unit for which the decision is made. The method is
therefore particularly interesting for selecting units from a stream. The
great advantage of this new method is that it forces the decision on the units
according to their order of appearance in the stream.

Section 2 introduces the notations and the principal concept of sampling
theory. In Section 3, we present the proposed sampling method while Section 4
extends the method to stream sampling. Section 5 and Section 6 discusses
different limit cases of the method and how the problems are solved. Section 7
exemplifies the method on a real data set and shows that the method is comple-
tely comparable in terms of variance. The manuscript ends with a conclusion on
the proposed method.

##  2 Notation

Consider a finite population of size $N$ labeled by $U=\\{1,2,\dots,N\\}$, a
vector of values taken by a variable of interest
$\bm{y}=(y_{1},y_{2},\dots,y_{N})^{\top}$ and a vector of first-order
inclusion probabilities $\bm{\pi}=(\pi_{1},\pi_{2},\dots,\pi_{N})^{\top}$,
each associated with the respective labels in $U$. Consider also the total of
the variable of interest,

$Y=\sum_{k\in U}y_{k},$

as the main parameter. A sample $s$ is a subset of $U$ and a sampling design
$p(.)$ is a probability distribution on all the subsets of $U$. A random
sample $S$ is a vector that maps all subsets $s$ to an $N$ vector of 0 or 1
such that

$P(S=s)=p(s),~{}p(s)\geq 0,\mbox{ and }\sum_{s\subset U}p(s)=1,\mbox{ for all
}s\subset U.$

In order to estimate this parameter, two approaches can be used

  1. 1.

Consider a sampling design $p(.)$ and find a sampling method (or algorithm) to
implement this design. Based on this approach, according to $p(.)$, all the
first and second-order inclusion probabilities can be specified using

$\begin{array}[]{lll}\pi_{k}=\textrm{Pr}(k\in S)=\sum_{s\ni k}p(s)\mbox{ and
}\\\ \pi_{k\ell}=\textrm{Pr}(k,\ell\in S)=\sum_{s\ni\\{k,\ell\\}}p(s),\mbox{
for all }k,\ell\in U.\end{array}$

  2. 2.

Consider a vector of first-order inclusion probabilities $\bm{\pi}$ and find a
sampling method to select a random sample $S$ such that $\textrm{Pr}(k\in
S)=\pi_{k},$ for all $k\in U$.

In this article, the second approach is considered. If $\pi_{k}>0$, for all
$k\in U$, $Y$ can be estimated unbiasedly by means of the first-order
inclusion probabilities using the expansion estimator (Narain,, 1951; Horvitz
and Thompson,, 1952)

$\widehat{Y}=\sum_{k\in S}\frac{y_{k}}{\pi_{k}}.$

Its precision depends on the second-order inclusion probabilities

$\textrm{Var}(\widehat{Y})=\sum_{k\in U}\sum_{\ell\in
U}(\pi_{k\ell}-\pi_{k}\pi_{\ell})\frac{y_{k}}{\pi_{k}}\frac{y_{\ell}}{\pi_{\ell}}.$

Indeed, the values of $\bm{\pi}$ indicates the desire chance (probability) for
the units to be selected in the final sample. These probabilities can be
calculated based on some proper auxiliary variables. If $\pi_{k\ell}>0,$ for
all $k,\ell\in U,$ $\textrm{Var}(\widehat{Y})$ can be unbiasedly estimated by

$\widehat{\textrm{Var}}(\widehat{Y})=\sum_{k\in S}\sum_{\ell\in
S}\frac{\pi_{k\ell}-\pi_{k}\pi_{\ell}}{\pi_{k\ell}}\frac{y_{k}}{\pi_{k}}\frac{y_{\ell}}{\pi_{\ell}},$

however this does not guaranty to have a precise estimation. According to the
sampling method, the precision of the expansion estimator can be affected. As
the last point in this section, it is notable that precision is not always the
most important criteria to evaluate a strategy (set of a sampling method and
an estimator). Other important criteria can be applicability of the strategy
according to accessibility of the data over time that will be discussed more
in Section 4.

##  3 One-Step One-Decision Sampling Method

In this section, we present the procedure called One-Step One-Decision (OSOD)
sampling method. The general idea of the method is that at each step a
complete decision is made on the considered unit. This means that the unit is
selected in the sample or not. After that, the inclusion probabilities are
updated according to the decision made on this unit. This procedure is then
repeated on the following units and it takes at most $N$ steps to obtain the
final sample.

Consider a finite population $U$. Define the sum of the inclusion
probabilities

$n=\sum_{k\in U}\pi_{k},$

that can be non-integer. For simplicity, we first consider the method where
the decision is made for the first unit in the list. Suppose also that

$\sum_{k=1}^{N}\pi_{k}\geq 1\mbox{ and }\sum_{k=1}^{N}(1-\pi_{k})\geq 1.$

The first unit is selected with probability $\pi_{1}$. Let $\pi_{1}^{1}$ be
the updated inclusion probability of the unit 1. It can be 0 or 1 and can be
summarize as follow

$\displaystyle\pi^{1}_{1}=\left\\{\begin{array}[]{ll}1&\mbox{ with probability
}\pi_{1}\\\ 0&\mbox{ with probability }1-\pi_{1}.\end{array}\right.$ (3)

After deciding about the first unit, the remaining inclusion probabilities
need to be updated. The whole idea of the method is to increase (respectively
decrease) inclusion probabilities of the following units depending if
$\pi_{1}^{1}$ has been changed to 0 (respectively 1). The inclusion
probabilities of the remaining units are updated as follow

$\pi^{1}_{k}=\left\\{\begin{array}[]{ll}\pi^{1}_{k}(0)=\min(c_{1}\pi_{k},1)&\mbox{
if }\pi^{1}_{1}=0\\\\[5.69054pt]
\displaystyle\pi^{1}_{k}(1)=\frac{\pi_{k}-\pi^{1}_{k}(0)(1-\pi_{1})}{\pi_{1}}&\mbox{
if }\pi^{1}_{1}=1,\end{array}\right.\text{ for }k=2,\dots,N,$

where constant $c_{1}$ is defined by

$\sum_{k=2}^{N}\min(c_{1}\pi_{k},1)=n.$ (4)

Equation (4) shows that the inclusion probabilities of the remaining units are
increased proportionally such that the sample size is still respected. After
the first step, the decision is irrevocably made for unit 1. Simply, this
operation can be repeated for other steps $t=2,\dots,N.$ At step $t$,
decisions are made about the first $t-1$ units and $\bm{\pi}$ has been updated
$t-1$ times. The updated vector is denoted by

$\bm{\pi}^{t-1}=(\pi^{t-1}_{1},\pi^{t-1}_{2},\dots,\pi^{t-1}_{k},\dots,\pi^{t-1}_{N})^{\top},$

where its first $t-1$ units are in $\\{0,1\\}$. Again, it is supposed that

$\sum_{k=t}^{N}\pi^{t-1}_{k}\geq 1\mbox{ and
}\sum_{k=t}^{N}(1-\pi^{t-1}_{k})\geq 1,$ (5)

and a decision is taken for the $t$th unit with probability

$\displaystyle\pi^{t}_{t}=\left\\{\begin{array}[]{ll}1&\mbox{ with probability
}\pi^{t-1}_{t}\\\ 0&\mbox{ with probability
}1-\pi^{t-1}_{t}.\end{array}\right.$ (8)

Then, the inclusion probabilities are updated by

$\pi^{t}_{k}=\left\\{\begin{array}[]{ll}\pi^{t}_{k}(0)=\min(c_{t}\pi^{t-1}_{k},1)&\mbox{
if }\pi^{t}_{t}=0\\\\[5.69054pt]
\displaystyle\pi^{t}_{k}(1)=\frac{\pi^{t-1}_{k}-\pi^{t}_{k}(0)(1-\pi^{t-1}_{t})}{\pi^{t-1}_{t}}&\mbox{
if }\pi^{t}_{t}=1,\end{array}\right.\text{ for }k=t+1,t+2\dots,N,$

where $c_{t}$ is defined by

$\sum_{k=t+1}^{N}\min(c_{t}\pi^{t-1}_{k},1)=n-n_{t},$

and $n_{t}$ is the number of selected units up to step $t$.

Conditions (5) are necessary to leave the possibility of updating the other
inclusion probabilities. For example at the first step it is required to have
enough space in other inclusion probabilities to distribute $\pi_{1}$ as
$\pi^{1}_{k}(0)=\pi_{k}+\alpha_{k}\pi_{1}$ (if unit $1$ is not selected) or to
remove $1-\pi_{1}$ from them as $\pi^{1}_{k}(1)=\pi_{k}-\alpha_{k}(1-\pi_{1})$
(if unit $1$ is selected) for some $0<\alpha_{k}\leq 1$ on the following
inclusion probabilities. Then, hereafter, OSOD stands for the design with a
population that satisfies Conditions (5).

Result 1 shows that the sampling procedure respects the inclusion
probabilities and thus the fixed sample size. The procedure is presented more
formally in Algorithm 1. Currently, the sampling procedure takes into account
the whole population to modify the inclusion probabilities. In the next
section, a great improvement for sampling from a stream is proposed. Indeed,
only a small part of the following units can be considered to update the
inclusion probabilities. This completely modifies the application of the
method.

######  Result 1.

With OSOD, for all $t=1,2,\dots$ we have

$\textrm{E}(\pi^{t}_{k})=\pi_{k},\;\mbox{ for all }k\in
U,\;\;\;\;\;\;\mbox{and}\;\;\;\;\;\;\sum_{k=1}^{N}\pi^{t}_{k}=\sum_{k=1}^{N}\pi_{k}.$

The proof is given in Appendix.

Algorithm 1 One-Step One-Decision sampling algorithm with parameter
$\bm{\pi}$.

Initialize with $n_{0}=0$, $\bm{\pi}^{0}=\bm{\pi}$ and $n=\sum_{k\in
U}\pi_{k}$.

for $t=1,2,3,\dots,N$ do

Calculate $c_{t}$ such that
$\displaystyle\sum_{k=t+1}^{N}\min(c_{t}\pi^{t-1}_{k},1)=n-n_{t}$ where
$\displaystyle n_{t}=\sum_{k=1}^{t-1}\pi^{t-1}_{k}$,

Generate $u$, a realization of a uniform random variable in $[0,1]$,

if $u>\pi^{t}_{t}$ then

Set $\pi^{t}_{t}=0$ and $\pi^{t}_{k}=\min(c_{t}\pi^{t-1}_{k},1)$ for all
$k>t$.

else

Set $\pi^{t}_{t}=1$ and
$\pi^{t}_{k}=\frac{\pi^{t-1}_{k}-\min(c_{t}\pi^{t-1}_{k},1)(1-\pi^{t-1}_{t})}{\pi^{t-1}_{t}}$
for all $k>t$.

end if

if $\pi^{t}_{k}\in\\{0,1\\}$ for all $k=t+1,t+2,...,N$ then

Stop

end if

end for

##  4 Extending the Method to Stream Sampling

A data stream is a sequence of information arriving sequentially in time. When
sampling streams, usually a very large set of data coming in continuously.
Because of the limited space needed to store the data, it is desirable to
decide on the selection of units at the right time.

As pointed in Section 1, in the methods proposed so far, the decision about a
unit that is at the beginning of a stream can be made very late. Ideally, the
decision to select a unit or not from a stream must be taken as soon as
possible, which the OSOD method allows. The idea to extend the OSOD method to
stream sampling is to consider a window of units instead of the entire
population to update inclusion probabilities. Result 2 shows that the number
of units to consider depends on the inclusion probabilities on the following
units. Essentially, it is enough to wait until having a set of units such that
there exists a constant $c_{t}$ that satisfies $\pi_{t}\geq(1-1/c_{t})$.
Moreover, Result 3 shows also that a sufficient condition to update inclusion
probabilities is that the inclusion probabilities within the window sum up to
integer.

######  Result 2.

In OSOD, a necessary and sufficient condition for all $\pi^{t}_{k}(1)$ to be
non-negative is that

$\pi_{t}\geq 1-1/c_{t}.$

######  Result 3.

In OSOD, a sufficient condition for all $\pi^{t}_{k}(1)$ to be non-negative is
that $n=\sum_{k\in U}\pi_{k}$ is integer.

The proofs are given in Appendix.

Hence after making decision about the first unit, enough units with a $c_{1}$
that satisfies $\pi_{1}\geq(1-1/c_{1})$ is required. Also if there is an
integer window then it is possible to apply the method before. Still, it could
be possible that there is no window that sums up to integer and no constant
$c_{1}$ that satisfies the condition until the end of the stream. In this
case, it is possible to use Result 3 and the notion of phantom unit (Grafström
et al., 2012b, ) to artificially transform the sum of the inclusion
probabilities to an integer and finish the procedure. This step is well
explained in the Section 6. The method has many advantages for performing
sampling in streams, including

  * •

The decision about the units can be made based on their order,

  * •

To update the population inclusion probabilities after selecting a unit, the
whole population is not needed and only a small part of it (based on Result 2
or Result 3) is enough. In fact, in the case of streams, there is usually no
finite set called population and the data is generated online,

  * •

It is possible to make decision about one part of the population indepen-
dently of the other parts.

Algorithm 2 gives a detailed explanation of the proposed method. The next two
Sections discuss example and show intuitively why the proposed method spread
the units according to their order and the case where no constant $c_{t}$ can
be found for the entire population.

Algorithm 2 One-Step One-Decision Stream Sampling algorithm.

Initialize with $\pi_{1}^{0}=\pi_{1}$,

for $t=1,2,3,\dots$ do

Consider unit $t$,

if There is a window of units of size $m$ for which
$\pi^{t-1}_{t}\geq(1-1/c_{t})$ then

1) Set $n-n_{t}=\displaystyle\sum_{k=t}^{t+m-1}\pi^{t-1}_{k}$, $N=t+m-1$,

2) Implement Algorithm 1, just for one step at $t$.

else

1) Consider the largest possible window (all the available data), say of size
$M$,

2) Add a phantom unit $k=t+M$ with
$\pi^{t-1}_{t+M}=1-\left(\sum_{k=t}^{t+M-1}\pi^{t-1}_{k}-\left\lfloor\sum_{k=t}^{t+M-1}\pi^{t-1}_{k}\right\rfloor\right),$

3) Set $n_{t}=0$, $n=\displaystyle\sum_{k=t}^{t+M}\pi^{t-1}_{k}$ and $N=t+M$,

4) Implement Algorithm 1 from step $t$ for $M$ steps,

5) Delete the phantom unit with $\pi^{t+M-1}_{t+M}$,

6) Set $t=t+M$.

end if

end for

##  5 Spread the Units According to their Order

For some reasons, such as a high correlation between the main variable and the
indices, it is often desirable to select the sample spread on the order of the
indices. Let $m$ be the number of units in the window satisfying the condition
of Result 2. Except for very extreme inclusion probabilities, it appears easy
to find small $m$ in Algorithm 2.

To spread the sample, the parameter $m$ in Algorithm 2 is set as small as
possible. Select the smallest window in the sampling procedure is somehow
evocative of local pivotal method (Grafström et al., 2012a, ). Indeed with
such strategy, the respective unit fights with some closest units around and
has its inclusion probability modified to 1 or 0. This idea of creating
repulsion between nearest units is exactly the same if a small $m$ is chosen.
Indeed, for a particular unit, if the inclusion probability is changed to 1,
the procedure will decrease inclusion probabilities of the following
neighbours. If the window is small, the next unit is more likely to be
modified directly to 0 or in the next step. Then $m$ can be considered as a
leverage to control the level of spread of the sample. Again, one of the
advantages of OSOD over the classical methods is that it does not postpone the
decision on the units.

To clarify this strategy, consider

$\bm{\pi}=(0.5,0.5,0.3,0.1,0.6,0.7,0.3)^{\top}.$

Based on Result 3, there are at least three windows to make decision about the
first unit, $m=2,5,7$ (as the sum is integer for each of them). More formally
the problem can be presented as

$\bm{\pi}=(\underbrace{\underbrace{\underbrace{0.5,0.5,}_{\bm{w}_{1},\text{
with $m=2$}}\;\;\;\;\;\;\underbrace{0.3,0.1,0.6,}_{\bm{w}_{2},\text{ with
$m=3$}}}_{\bm{w}_{1},\text{ with
$m=5$}}\;\;\;\;\;\;\underbrace{0.7,0.3}_{\bm{w}_{3},\text{ with
$m=2$}}}_{\bm{w}_{1},\text{ with $m=7$}})^{\top}.$ (9)

If $m=2$, then one of the first two units will be selected in the sample. But
if $m=7$ is considered, it is possible to have both, or none of the first
units inside the final sample. Figure 1 shows the sampling process of the
method with the small windows ($m=2,\;m=3$ and $m=2$) strategy. For all the
cases, the sample is well spread in the three windows.

Considering all the population as a window (i.e. $m=7$) and apply the method,
will lead to obtain one of the sample in Table 1 (probabilities are based on
10,000 iterations). In this situation, it is possible to have a sample that is
not spread over the population at all, such as samples that are labelled in
the table by indices 6, 16, and 31.

$\scriptstyle\begin{pmatrix}0.5\\\ 0.5\\\ 0.3\\\ 0.1\\\ 0.6\\\ 0.7\\\ 0.3\end{pmatrix}$$\scriptstyle\begin{pmatrix}0\\\ 1\\\ 0.3\\\ 0.1\\\ 0.6\\\ 0.7\\\ 0.3\end{pmatrix}$$\scriptstyle\begin{pmatrix}1\\\ 0\\\ 0.3\\\ 0.1\\\ 0.6\\\ 0.7\\\ 0.3\end{pmatrix}$$\scriptstyle\begin{pmatrix}0\\\ 1\\\ 0\\\ 0.14\\\ 0.86\\\ 0.7\\\ 0.3\end{pmatrix}$$\scriptstyle\begin{pmatrix}0\\\ 1\\\ 1\\\ 0\\\ 0\\\ 0.7\\\ 0.3\end{pmatrix}$$\scriptstyle\begin{pmatrix}1\\\ 0\\\ 0\\\ 0.14\\\ 0.86\\\ 0.7\\\ 0.3\end{pmatrix}$$\scriptstyle\begin{pmatrix}1\\\ 0\\\ 1\\\ 0\\\ 0\\\ 0.7\\\ 0.3\end{pmatrix}$$\scriptstyle\begin{pmatrix}0\\\ 1\\\ 0\\\ 0\\\ 1\\\ 0.7\\\ 0.3\end{pmatrix}$$\scriptstyle\begin{pmatrix}0\\\ 1\\\ 0\\\ 1\\\ 0\\\ 0.7\\\ 0.3\end{pmatrix}$$\scriptstyle\begin{pmatrix}0\\\ 1\\\ 1\\\ 0\\\ 0\\\ 0\\\ 1\end{pmatrix}$$\scriptstyle\begin{pmatrix}0\\\ 1\\\ 1\\\ 0\\\ 0\\\ 1\\\ 0\end{pmatrix}$$\scriptstyle\begin{pmatrix}1\\\ 0\\\ 0\\\ 0\\\ 1\\\ 0.7\\\ 0.3\end{pmatrix}$$\scriptstyle\begin{pmatrix}1\\\ 0\\\ 0\\\ 1\\\ 0\\\ 0.7\\\ 0.3\end{pmatrix}$$\scriptstyle\begin{pmatrix}1\\\ 0\\\ 1\\\ 0\\\ 0\\\ 0\\\ 1\end{pmatrix}$$\scriptstyle\begin{pmatrix}1\\\ 0\\\ 1\\\ 0\\\ 0\\\ 1\\\ 0\end{pmatrix}$$\scriptstyle\begin{pmatrix}0\\\ 1\\\ 0\\\ 0\\\ 1\\\ 0\\\ 1\end{pmatrix}$$\scriptstyle\begin{pmatrix}0\\\ 1\\\ 0\\\ 0\\\ 1\\\ 1\\\ 0\end{pmatrix}$$\scriptstyle\begin{pmatrix}0\\\ 1\\\ 0\\\ 1\\\ 0\\\ 0\\\ 1\end{pmatrix}$$\scriptstyle\begin{pmatrix}0\\\ 1\\\ 0\\\ 1\\\ 0\\\ 1\\\ 0\end{pmatrix}$$\scriptstyle\begin{pmatrix}1\\\ 0\\\ 0\\\ 0\\\ 1\\\ 0\\\ 1\end{pmatrix}$$\scriptstyle\begin{pmatrix}1\\\ 0\\\ 0\\\ 0\\\ 1\\\ 1\\\ 0\end{pmatrix}$$\scriptstyle\begin{pmatrix}1\\\ 0\\\ 0\\\ 1\\\ 0\\\ 0\\\ 1\end{pmatrix}$$\scriptstyle\begin{pmatrix}1\\\ 0\\\ 0\\\ 1\\\ 0\\\ 1\\\ 0\end{pmatrix}$0.50.50.70.30.70.30.860.140.30.70.860.140.30.70.30.70.30.70.30.70.30.7$None$ $None$$None$$None$$None$$None$$None$$None$$None$$None$$None$$None$ Figure 1: Well spread strategy of OSOD on $\pi=(0.5,0.5,0.3,0.1,0.6,0.7,0.3)^{\top}$ with considering three small windows ($m_{1}=2,m_{2}=3$ and $m_{3}=2$) on the population in Equation (9). The boxed numbers indicate the design of all the possible samples. Each time at least a complete decision about one unit is taken and then just the inclusion probabilities inside the respective window are updated (and not the whole population).  Table 1: Implementing OSOD on $\bm{\pi}=(0.5,0.5,0.3,0.1,0.6,0.7,0.3)^{\top}$ with considering the whole population as the only window. The design is calculated based on 10,000 iterations. No. of Sample | The Samples | Sampling Design
---|---|---
1 | 0 | 0 | 1 | 0 | 0 | 1 | 1 | 0.0123
2 | 0 | 0 | 1 | 0 | 1 | 0 | 1 | 0.0069
3 | 0 | 0 | 1 | 0 | 1 | 1 | 0 | 0.0369
4 | 0 | 0 | 1 | 1 | 0 | 0 | 1 | 0.0078
5 | 0 | 0 | 1 | 1 | 0 | 1 | 0 | 0.0148
6 | 0 | 0 | 1 | 1 | 1 | 0 | 0 | 0.0144
7 | 0 | 1 | 0 | 0 | 0 | 1 | 1 | 0.0582
8 | 0 | 1 | 0 | 0 | 1 | 0 | 1 | 0.0308
9 | 0 | 1 | 0 | 0 | 1 | 1 | 0 | 0.1537
10 | 0 | 1 | 0 | 1 | 0 | 0 | 1 | 0.0024
11 | 0 | 1 | 0 | 1 | 0 | 1 | 0 | 0.0039
12 | 0 | 1 | 0 | 1 | 1 | 0 | 0 | 0.0032
13 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 0.0266
14 | 0 | 1 | 1 | 0 | 0 | 1 | 0 | 0.0651
15 | 0 | 1 | 1 | 0 | 1 | 0 | 0 | 0.0554
16 | 0 | 1 | 1 | 1 | 0 | 0 | 0 | 0.0027
17 | 1 | 0 | 0 | 0 | 0 | 1 | 1 | 0.0804
18 | 1 | 0 | 0 | 0 | 1 | 0 | 1 | 0.0415
19 | 1 | 0 | 0 | 0 | 1 | 1 | 0 | 0.1953
20 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | 0.0074
21 | 1 | 0 | 0 | 1 | 0 | 1 | 0 | 0.0171
22 | 1 | 0 | 0 | 1 | 1 | 0 | 0 | 0.0178
23 | 1 | 0 | 1 | 0 | 0 | 0 | 1 | 0.0085
24 | 1 | 0 | 1 | 0 | 0 | 1 | 0 | 0.0216
25 | 1 | 0 | 1 | 0 | 1 | 0 | 0 | 0.0186
26 | 1 | 0 | 1 | 1 | 0 | 0 | 0 | 0.0030
27 | 1 | 1 | 0 | 0 | 0 | 0 | 1 | 0.0153
28 | 1 | 1 | 0 | 0 | 0 | 1 | 0 | 0.0351
29 | 1 | 1 | 0 | 0 | 1 | 0 | 0 | 0.0324
30 | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 0.0048
31 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0.0061

##  6 Populations with Non-Integer Sum of Inclusion Probabilities

As pointed before, the method is perfectly working if the sum of inclusion
probabilities within the window is not an integer as long as the condition of
Result 2 is satisfied. An issue occurs if the condition is not satisfied and
the sum of the inclusion probabilities is not an integer. In such situation,
this problem can be solved by using the notion of phantom unit (Grafström et
al., 2012b, ). The ingenuity is to add a phantom unit that has inclusion
probability equal to the difference needed to reach the ceiling integer. Then,
by Result 3, the method is working and at the end of the process the phantom
unit will be removed of the sample. The strategy respects the inclusion
probabilities. As an example, consider the following inclusion probability
vector

$\displaystyle\small\begin{array}[]{ll}\bm{\pi}=(0.85,0.90,0.90,0.02,00.02,0.98,0.99,0.95,0.99,0.01,0.01,0.99,0.99,0.99)^{\top}.\end{array}$
(11)

In this case, because of extreme inclusion probabilities, the condition of
Result 2 is not satisfied on the whole population and the sum of the inclusion
probabilities is not an integer. With implementing the method, the updated
inclusion probabilities are (rounded to 2 decimal places)

$\displaystyle\bm{\pi}^{1}=\small\left\\{\begin{array}[]{llll}(0.00,\;1.00,\;1.00,\;0.20,\;0.20,\;1.00,\;1.00,\\\
\hskip 28.45274pt1.00,\;1.00,\;0.10,\;0.10,\;1.00,\;1.00,\;1.00)^{\top}&\mbox{
if }\pi^{1}_{1}=0\\\ \\\
\displaystyle(1.00,0.88,0.88,-0.01,-0.01,0.98,0.99,\\\ \hskip
28.45274pt0.94,0.99,-0.01,-0.01,0.99,0.99,0.99)^{\top}&\mbox{ if
}\pi^{1}_{1}=1,\end{array}\right.$ (17)

Because of negative values, $\bm{\pi}^{1}(1)$ is not acceptable. As the sum of
inclusion probabilities is $9.59$, the proposed solution is to add a phantom
unit with inclusion probability $\pi^{*}=0.41$. Based on Result 3 it is
possible to apply the method on

which leads to updated vectors as

$\displaystyle\bm{\pi}^{*1}=\small\left\\{\begin{array}[]{llll}(0.00,1.00,1.00,0.04,0.04,1.00,1.00,\\\
\hskip
28.45274pt1.00,1.00,0.02,0.02,1.00,1.00,1.00,\boxed{0.87})^{\top}&\mbox{ if
}\pi^{1}_{1}=0\\\ \\\ \displaystyle(1.00,0.88,0.88,0.02,0.02,0.98,0.99,\\\
\hskip
28.45274pt0.94,0.99,0.01,0.01,0.99,0.99,0.99,\boxed{0.33})^{\top}&\mbox{ if
}\pi^{1}_{1}=1.\end{array}\right.$ (23)

As the sum remains constant in the updated inclusion probabilities, the method
can be applied step by step to decide about all the units. At the last step,
the phantom unit will be ignored. With this procedure, inclusion probabilities
and also the expectation of the sample size are respected. With this strategy
it is always possible to finish the selection of the sample even if the two
necessary conditions are not fulfilled. It is important to note that if this
strategy is applied while the selection process is not at the end of its
course, the sample size is random with expectation equal to the sum of the
inclusion probabilities. For more information, the reader can find details in
(Grafström et al., 2012b, ).

##  7 Simulations

As the method has significant new advantages over other methods, we show on
the basis of a real dataset, that in terms of estimation accuracy, the method
is comparable to its competitors. The data set is available in the R package
sampling (Tillé and Matei,, 2015) under the name belgianmunicipalities. It
contains different information about the Belgian population between 2001 and
2004 by municipalities. The data set contains 589 rows corresponding to each
municipality of Belgium. Here just two variables were considered, the variable
of interest $y$ which represents the taxable income in euros in 2001, and the
auxiliary variable $x$ which provides the total population on July 1, 2004.
Figure 2 shows the variable of interest by municipalities.

Next, based on 10,000 simulations, several sampling methods were considered
and expansion estimator for each of them was calculated. OSOD method was
compared to five other sampling methods, random systematic (Madow,, 1949),
random pivotal (Deville and Tillé,, 1998), Tillé’s method (Tillé,, 1996),
Midzuno’s method (Midzuno,, 1950) and finally the conditional Poisson sampling
also called the maximum entropy method (Chen and Liu,, 1997). The reader can
find more details on each of these methods in (Tillé,, 2006).

The sample size was set to 200 municipalities and the inclusion probabilities
were unequal and set proportional to the total population variable $x$. Figure
3 shows the boxplot of the expansion estimator of the total for each sampling
design and confirms that our method is completely comparable in terms of
variance and respect the unequal inclusion probabilities. The method was
imple-mented in the R package StratifiedSampling (Jauslin et al.,, 2021) under
the name osod.

$49.5^{\circ}$N$50.0^{\circ}$N$50.5^{\circ}$N$51.0^{\circ}$N$51.5^{\circ}$N$3^{\circ}$E$4^{\circ}$E$5^{\circ}$E$6^{\circ}$ELongituteLatitudeTaxable
Income[0, $1.5\times 10^{6}$]($1.5\times 10^{6}$, $50\times 10^{6}$]($50\times
10^{6}$, $90\times 10^{6}$]($90\times 10^{6}$, $120\times 10^{6}$]($120\times
10^{6}$, $164\times 10^{6}$]($164\times 10^{6}$, $228\times
10^{6}$]($228\times 10^{6}$, $400\times 10^{6}$]($400\times 10^{6}$, $\infty$)
Figure 2: Belgian municipalities dataset. The discrete gradient scale
represents the taxable income in euros.
1.18e+111.20e+111.22e+111.24e+11OsodRandom systematicRandom
pivotalTilleMidzunoMax entropysampling designTaxable income Figure 3: Boxplot
of the Horvitz-Thompson estimator for the variable of interest taxable income
of the Belgian municipalities dataset. Different sampling design are
illustrated as well as the proposed method. We find that the proposed sampling
design is quite comparable in terms of variance estimator.

##  8 Conclusion

The OSOD method is the first method in unequal probability sampling that can
decide about the units completely based on their order. The method is flexible
on the order of the data and spreading the data over the population. The order
can be arbitrary but if the data are the result of a stream, the method can be
very useful to decide about what units should be stored as a sample over time.
Then after deciding about each unit to update the inclusion probabilities,
there is no need to use all the population and usually just a small window of
the population can be enough. Also based on the size of the windows, it is
possible to control sample rate of spreading over the population indices. The
indices can indicate time or location of the data. If the main variable is
correlated with the indices, it would be very beneficial to select the windows
as small as possible.

## Appendix

###### Proof of Result 1.

Without loss of generality and for ease of notation we only do the proof for
$t=1$. For $k=1$ it is obvious that $\textrm{E}(\pi^{1}_{1})=\pi_{1}$, and for
$k=2,3,\dots,N$,

$\displaystyle\textrm{E}(\pi^{1}_{k})=\pi^{1}_{k}(0)(1-\pi_{1})+\pi^{1}_{k}(1)\pi_{1}=\pi_{k}^{1}(0)(1-\pi_{1})+\frac{\pi_{k}-\pi^{1}_{k}(0)(1-\pi_{1})}{\pi_{1}}\pi_{1}=\pi_{k}.$

Also for the sum of the inclusion probabilities, we have

$\sum_{k\in
U}\pi_{k}^{1}(0)=0+\sum_{k=2}^{N}\pi^{1}_{k}(0)=\sum_{k=2}^{N}\min(c_{1}\pi_{k},1)=n=\sum_{k=1}^{N}\pi_{k},$

and

$\begin{array}[]{lll}\displaystyle\sum_{k\in U}\pi_{k}^{1}(1)&=&\displaystyle
1+\sum_{k=2}^{N}\pi^{1}_{k}(1)\\\ &=&\displaystyle
1+\frac{\sum_{k=2}^{N}\pi_{k}-(1-\pi_{1})\sum_{k=2}^{N}\pi^{1}_{k}(0)}{\pi_{1}}\\\
&=&\displaystyle\frac{\pi_{1}+\\{(n-\pi_{1})-n+n\pi_{1}\\}}{\pi_{1}}\\\
&=&n\\\ &=&\displaystyle\sum_{k=1}^{N}\pi_{k}.\end{array}$

∎

###### Proof of Result 2.

Without loss of generality and for ease of notation we only do the proof for
$t=1$. Since

$\displaystyle\pi^{1}_{k}(1)=\left\\{\begin{array}[]{ll}\displaystyle\frac{\pi_{k}-\min(c_{1}\pi_{k},1)(1-\pi_{1})}{\pi_{1}}\leq\frac{\pi_{k}-c_{1}\pi_{k}(1-\pi_{1})}{\pi_{1}}&\mbox{
if }\pi^{1}_{k}(1)=1\\\\[8.53581pt]
\displaystyle\frac{\pi_{k}-c_{1}\pi_{k}(1-\pi_{1})}{\pi_{1}}&\mbox{ if
}\pi_{k}^{1}(1)<1,\end{array}\right.$ (26)

for $k=2,3,\dots,N$, a necessary and sufficient condition is that $\pi_{1}\geq
1-1/c_{1}$, which gives the result. ∎

###### Proof of Result 3.

Let

$U_{A}=\\{k\in U|0<\pi_{k}^{1}(0)<1\\},\;U_{B}=\\{k\in U|\pi_{k}^{1}(0)=1\\},$
$A=\sum_{k\in U_{A}}\pi_{k}\mbox{ and }B=\sum_{k\in U_{B}}\pi_{k}.$

Then it is possible to decompose $n$ as

$\pi_{1}+A+B=n$ (27)

and

$c_{1}A+\\#U_{B}=n.$

Now, if $n$ is integer, $c_{1}A$ is an integer denoted by $d$. Thus
$A=d/c_{1}$. In this case, $\\#U_{B}=n-d$. Moreover, it is easy to see that
$B\leq\\#U_{B}$. Now, from Equation (27), we have

$\pi_{1}=n-A-B=n-\frac{d}{c_{1}}-B\geq
n-\frac{d}{c_{1}}-\\#U_{B}=n-\frac{d}{c_{1}}-(n-d)=d\left(1-\frac{1}{c_{1}}\right).$

Now if $d=1,2,\dots$ then $\pi_{1}\geq(1-1/c_{1})$ and if $d=0$ then
$\\#U_{B}=n$ or in other words, $\pi^{1}_{k}(0)=1$ for all $k=2,3,\dots,N$.
Therefore to have all $\pi^{1}_{k}(1)\geq 0$ we need to have

$\pi^{1}_{k}(1)=\frac{\pi_{k}-(1-\pi_{1})}{\pi_{1}}\geq
0,\Rightarrow\pi_{k}+\pi_{1}-1\geq 0.$ (28)

But as in such cases, $\pi^{1}_{k}(0)=\pi_{k}+\alpha_{k}\pi_{1}=1,$ for some
$0<\alpha_{k}\leq 1$, then $\pi_{k}+\pi_{1}\geq 1$ and therefore Condition
(28) is satisfied.
Then, if $n$ is integer, the condition of Result 2 is always fulfilled. ∎

## References

  * Busnel and Tillé, (2020) Busnel, Y. and Tillé, Y. (2020).  Attack-tolerant unequal probability sampling methods over sliding window for distributed streams.  In 4th International Conference on Compute and Data Analysis (ICCDA 2020), Mar 2020, San Jose, United States, pages 72–78.
  * Chao, (1982) Chao, M.-T. (1982).  A general purpose unequal probability sampling plan.  Biometrika, 69:653–656.
  * Chauvet, (2012) Chauvet, G. (2012).  On a characterization of ordered pivotal sampling.  Bernoulli, 18(4):1320–1340.
  * Chen and Liu, (1997) Chen, X.-H. and Liu, J. S. (1997).  Statistical applications of the Poisson-binomial and conditional Bernoulli distributions.  Statistica Sinica, 7:875–892.
  * Cohen et al., (2009) Cohen, E., Duffield, N., Kaplan, H., Lund, C., and Thorup, M. (2009).  Stream sampling for variance-optimal estimation of subset sums.  In Proceedings of the twentieth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 1255–1264. Society for Industrial and Applied Mathematics.
  * Deville, (1998) Deville, J.-C. (1998).  Une nouvelle (encore une!) méthode de tirage à probabilités inégales.  Technical Report 9804, Méthodologie Statistique, Insee, Paris.
  * Deville and Tillé, (1998) Deville, J.-C. and Tillé, Y. (1998).  Unequal probability sampling without replacement through a splitting method.  Biometrika, 85:89–101.
  * Gabler, (1990) Gabler, S. (1990).  Minimax Solutions in Sampling from Finite Populations.  Springer, New York.
  * (9) Grafström, A., Lundström, N. L. P., and Schelin, L. (2012a).  Spatially balanced sampling through the pivotal method.  Biometrics, 68(2):514–520.
  * (10) Grafström, A., Matei, A., Qualité, L., and Tillé, Y. (2012b).  Size constrained unequal probability sampling with a non-integer sum of inclusion probabilities.  Electronic Journal of Statistics, 6:1477–1489.
  * Hanif and Brewer, (1980) Hanif, M. and Brewer, K. R. W. (1980).  Sampling with unequal probabilities without replacement: A review.  International Statistical Review, 48:317–335.
  * Horvitz and Thompson, (1952) Horvitz, D. G. and Thompson, D. J. (1952).  A generalization of sampling without replacement from a finite universe.  Journal of the American Statistical Association, 47(260):663–685.
  * Jauslin et al., (2021) Jauslin, R., Eustache, E., Panahbehagh, B., and Tillé, Y. (2021).  StratifiedSampling: Different Methods for Stratified Sampling.  R Foundation for Statistical Computing, Vienna, Austria.  R package version 0.3.0.
  * Madow, (1949) Madow, W. G. (1949).  On the theory of systematic sampling, II.  Annals of Mathematical Statistics, 20:333–354.
  * Midzuno, (1950) Midzuno, H. (1950).  An outline of the theory of sampling systems.  Annals of the Institute of Statistical Mathematics, 1:149–156.
  * Narain, (1951) Narain, R. D. (1951).  On sampling without replacement with varying probabilities.  Journal of the Indian Society of Agricultural Statistics, 3:169–174.
  * Tillé, (1996) Tillé, Y. (1996).  An elimination procedure of unequal probability sampling without replacement.  Biometrika, 83:238–241.
  * Tillé, (2006) Tillé, Y. (2006).  Sampling Algorithms.  Springer, New York.
  * Tillé, (2019) Tillé, Y. (2019).  A general result for selecting balanced unequal probability samples from a stream.  Information Processing Letters, 152:1–6.
  * Tillé and Matei, (2015) Tillé, Y. and Matei, A. (2015).  sampling: Survey Sampling.  R package version 2.7.