Antithetic Sampling

Mohammad Jafari Jozani Université de Sherbrooke, Sherbrooke, Qc, CANADA

François Perron Université de Montréal, Montréal, Qc, CANADA

Summary. Consider the problem of estimating the mean of a population in a nonparametric framework. Assume that we can rank all of the observations but we are aloud to measure only n of them. In this paper we introduce new sampling algorithms using antithetic variables. We propose unbiased estimators which reduce the variance in comparison with the ranked set sampling estimator. Theoretical results and numerical comparisons are included.

Keywords: Antithetic sampling; Ranked set sampling.

1 Introduction

In most sampling surveys, a reasonable number of the sampling units can be fairly accurately ordered with respect to a variable of interest without actual measurement and at little cost. On the other hand, exact measurements of these units may be very tedious and/or expensive. For example, for environmental risks such as radiation (soil contamination, disease clusters) or pollution (water contamination, root disease of crops) we commonly find that exact measurements involves substantial scientific processing of materials and correspondingly high attendant cost, while the variable of interest from the experimental (sampling) units can easily be ranked. Ranked set sampling, as proposed by McIntyre (1952) in estimating the mean of the pasture yields, provides an interesting alternative to simple random sampling in these situations.

Compared to simple random sampling, ranked set sampling has been proven theoretically (e.g., Takahasi and Wakimoto, 1968) and shown empirically (see, Kaur et al., 1995 and Chen et al., 2004 for more details) to yield more precise estimator of the population mean. This is especially useful for small surveys or for situations where it is expensive or destructive to obtain data.

Address for correspondence: Département de Mathématiques et de Statistique, Université de Montréal, P.O. Box 6128, Succursale Centre-Ville, Montréal, Que., Canada H3C 3J7.

E-mail address: perronf@dms.umontreal.ca (F. Perron).

Since the pioneering articles of McIntyre (1952) and Takahasi and Wakimoto (1968) several variations of ranked set sampling method has been proposed and developed by researchers to come up with more efficient estimators of a population mean. For example, Samawi et al. (1996) introduced extreme ranked set sampling and obtained an unbiased estimator of the mean which outperforms the usual mean of a simple random sample of the same size for symmetric distributions. Muttlak (1997) suggested median ranked set sampling to increase the efficiency and to reduce ranking errors over ranked set sampling method and proved its better performance in estimating the mean of a variable of interest for some symmetric distributions. Some relevant references here, in addition to those mentioned above, are among others: Bhoj (1997) for a new parametric ranked set sampling, Li et al. (1999) for random selection in ranked set sampling, Hossain and Muttlak (1999) for paired ranked set sampling, Al-Saleh and Al-Kadiri (2000) for double ranked set sampling, Hossain and Muttlak (2001) for selected ranked set sampling and Al-Saleh and Al-Omari (2002) for multistage ranked set sampling. Al-Nasser (2007) introduced L-Ranked set sampling design as a generalization of some of the above mentioned ranked set type sampling methods and proved the optimal property of his proposed estimators for symmetric family of distributions.

A review of these articles reveals that most of them concentrate mainly on situations where specific assumptions are made concerning the parent distribution of the underlying population and the proposed estimators are beneficiated from those assumptions to provide their optimal properties such as unbiasedness or having smaller variance.

The aim of this paper is to propose better estimators than μ̂_{R}, the ranked set sampling estimator, under the assumption that we can rank the observations before measurements, in a completely non parametric basis and without making any assumption about the distribution of the underlying population. The organization of the remaining sections is as follows. In Section 2, we introduce a new estimator, μ̂_{A}, of the population mean called the antithetic estimator, based on a new proposed ranked set type sampling method. We show that our estimator is unbiased and we derive an explicit expression for its variance. It is then proved that μ̂_{A} dominates μ̂_{S}, the simple random sampling estimator, by a direct comparison of the variances. However, it is not clear how to handle the comparisons of the variances between the μ̂_{R} and μ̂_{A} based simply on a direct calculation. In Section 3, we introduce a class of random estimators. The main idea is conditioning on the order statistics. It is interesting to see that μ̂_{S} and μ̂_{R} are now viewed as random estimators. A general dominance result leading to a sufficient condition (see Theorem 3.1), which we believe is novel, for a random estimator μ̂_{1} to dominate another random estimator μ̂_{2} is established. This method works for showing that μ̂_{A} dominates μ̂_{S}. Similarly, we show that μ̂_{A} dominates μ̂_{R} when the sample size is equal to two. Finally, in Section 4, we propose a better estimator than μ̂_{R} for all sample sizes. A simulation study is carried out in Section 5.

2 The antithetic estimator

Consider X_{1}, …, X_{n^{2}} as a simple random sample of size n^{2} from a distribution with mean μ and finite variance σ^{2} with X_{(1)}, …, X_{(n^{2})} as its corresponding order statistics. Let X_{ij} = X_{n(i−1)+j} and define X_{(ij)} = X_{(n(i−1)+j)}, i, j = 1, …, n, be the n(i − 1) + j th order statistics in the sample of size n^{2}. Also, we define X_{i(j)} to be the j th order statistics based on X_{i1}, …, X_{in}, i, j = 1, …, n. We adopt the following notations in the rest of the paper:

X_{(i·)} = \frac{1}{n} ∑_{j} X_{(ij)}, X_{(·j)} = \frac{1}{n} ∑_{i} X_{(ij)}, X_{(··)} = \frac{1}{n^{2}} ∑_{ij} X_{(ij)}, i, j = 1, …, n.

μ = E[X_{1}], σ^{2} = Var[X_{1}], μ_{1(i)} = E[X_{1(i)}], σ_{1(i)}^{2} = Var[X_{1(i)}], i = 1, …, n.

μ_{(ij)} = E[X_{(ij)}], μ_{(i·)} = E[X_{(i·)}], μ_{(·j)} = E[X_{(·j)}], μ_{(··)} = E[X_{(··)}], i, j = 1, …, n.

σ_{(ij)}^{2} = Var[X_{(ij)}], σ_{(i·)}^{2} = Var[X_{(i·)}], σ_{(·j)}^{2} = Var[X_{(·j)}], σ_{(··)}^{2} = Var[X_{(··)}], i, j = 1, …, n.

For estimating the population mean μ, let μ̂_{S} be the mean of a simple random sample of size n, i.e., {X_{1(1)}, …, X_{1(n)}}, which is given by

μ̂_{S} = \frac{1}{n} ∑_{i=1}^{n} X_{1(i)}.  (1)

with E[μ̂_{S}] = μ and Var[μ̂_{S}] = \frac{σ^{2}}{n}. Also, let μ̂_{R} be the mean of a ranked set sample of size n, i.e., {X_{1(1)}, …, X_{n(n)}}, from the underlying population as

μ̂_{R} = \frac{1}{n} ∑_{i=1}^{n} X_{i(i)}.  (2)

It is shown that (e.g., Takahasi and Wakimoto, 1968; Dell, 1969; Dell and Clutter, 1972), E[μ̂_{R}] = μ with

Var[μ̂_{R}] = \frac{σ^{2}}{n} − \frac{1}{n^{2}} ∑_{i=1}^{n}(μ_{1(i)} − μ)^{2}.  (3)

Let π = {π_{1}, …, π_{n}} be a random permutation of {1, …, n} independent of the sample and select {X_{(1π_{1})}, …, X_{(nπ_{n})}} as a new ranked set type sample called the antithetic sample. For estimating the population mean, based on the antithetic sample, we propose the following estimator

μ̂_{A} = \frac{1}{n} ∑_{i=1}^{n} X_{(iπ_{i})}.  (4)

Some basic properties of μ̂_{A} are given in Theorem 2.1.

Theorem 2.1 Let μ̂_{A} be the antithetic estimator defined in (4) and μ̂_{S} be the sample mean of a simple random sample of the same size as in (1). Then, μ̂_{A} is an unbiased estimator of the population mean μ and

\frac{σ^{2}}{n^{2}} ≤ Var[μ̂_{A}] = \frac{σ^{2}}{n} + \frac{1}{n(n − 1)} [2σ^{2} − ∑_{k=1}^{n}\{σ_{(k·)}^{2} + (μ_{(k·)} − μ)^{2} + σ_{(·k)}^{2} + (μ_{(·k)} − μ)^{2}\}] ≤ \frac{σ^{2}}{n} .

Proof. Let Y_{ij} = X_{(ij)} − μ. ConsiderP, the class of all permutations on {1, …, n}. We obtain

nE[μ̂_{A} − μ] = \frac{1}{n!} ∑_{π∈P} E[∑_{i} Y_{iπ_{i}}]
= E[\frac{1}{n} ∑_{i k} Y_{ik}]
= nE[X̄ − μ]
= 0.

Furthermore, to obtain Var(μ̂_{A}), we have

n^{2}Var(μ̂_{A}) = \frac{1}{n!} ∑_{π∈P} E[\{∑_{i} Y_{iπ_{i}}\}^{2}]

= \frac{1}{n!} ∑_{π∈P} E[∑_{i} Y_{iπ_{i}}^{2} + ∑_{i≠j} Y_{iπ_{i}}Y_{jπ_{j}}]

= E[\frac{1}{n} ∑_{i k} Y_{ik}^{2} + \frac{1}{n(n − 1)} ∑_{i≠j} ∑_{k≠ℓ} Y_{ik}Y_{jℓ}]

= E[\frac{1}{n} ∑_{i k} Y_{ik}^{2} + \frac{1}{n(n − 1)}\{ ∑_{i j k ℓ} Y_{ik}Y_{jℓ} − ∑_{i j k} Y_{ik}Y_{jk} − ∑_{i k ℓ} Y_{ik}Y_{iℓ} + ∑_{i k} Y_{ik}^{2}\}]

= E[\frac{1}{n} ∑_{i k} Y_{ik}^{2} + \frac{1}{n(n − 1)}\{(∑_{i k} Y_{ik})^{2} − ∑_{k}(∑_{i} Y_{ik})^{2} − ∑_{i}(∑_{k} Y_{ik})^{2} + ∑_{i k} Y_{ik}^{2}\}]

= E[\frac{1}{n} ∑_{i k} Y_{ik}^{2} + \frac{1}{n(n − 1)}\{n^{4}Y_{··}^{2} − n^{2} ∑_{k} Y_{·k}^{2} − n^{2} ∑_{k} Y_{k·}^{2} + ∑_{i k} Y_{ik}^{2}\}]

= E[n^{2}Y_{··}^{2} + \frac{1}{(n − 1)} ∑_{i k}\{Y_{ik} − Y_{·k} − Y_{i·} + Y_{··}\}^{2}]

= σ^{2} + \frac{1}{(n − 1)}E[∑_{i j}\{Y_{ij} − Y_{i·} − Y_{·j} + Y_{··}\}^{2}]

= σ^{2} + \frac{1}{(n − 1)} ∑_{i j}[\{σ_{(ij)}^{2} − σ_{(i·)}^{2} − σ_{(·j)}^{2} + σ_{(··)}^{2}\} + \{μ_{(ij)} − μ_{(i·)} − μ_{(·j)} + μ_{(··)}\}^{2}]
= \frac{n}{(n − 1)}[(n + 1)σ^{2} − ∑_{k}\{σ_{(k·)}^{2} + (μ_{(k·)} − μ)^{2} + σ_{(·k)}^{2} + (μ_{(·k)} − μ)^{2}\}],

which leads to the result. In particular, note that

Var(μ̂_{A}) = \frac{σ^{2}}{n^{2}} + \frac{1}{n^{2}(n − 1)}E[∑_{i j}\{Y_{ij} − Y_{i·} − Y_{·j} + Y_{··}\}^{2}] ≥ \frac{σ^{2}}{n^{2}} .

In order to show that Var(μ̂_{A}) ≤ \frac{σ^{2}}{n} we will verify that

\frac{1}{n − 1} ∑_{i k}\{Y_{ik} − Y_{·k} − Y_{i·} + Y_{··}\}^{2} ≤ \frac{1}{n + 1} ∑_{i k}\{Y_{ik} − Y_{··}\}^{2},

and use the fact that

Var(∑_{i k}{Y_{ik} − Y_{··}}^{2}) = Var(∑_{i}(X_{i} − X̄)^{2}) = (n^{2} − 1)σ^{2}.

First of all, if α_{1} ≤ ⋯ ≤ α_{n} and β_{1} ≤ ⋯ ≤ β_{n} then ∑_{k}(α_{k} − ᾱ)(β_{k} − β̄) ≥ 0. This implies that

∑_{j}(Y_{ij} − Y_{i·})(Y_{kj} − Y_{k·}) ≥ 0, i, k = 1, …, n,

and

∑_{k}(Y_{ki} − Y_{k·})(Y_{kj} − Y_{k·}) ≥ 0, i, j = 1, …, n.

Notice also that a_{1} ≤ a_{2} and b_{1} ≤ b_{2} implies that (a_{1} + b_{1})^{2} + (a_{2} + b_{2})^{2} ≥ (a_{1} + b_{2})^{2} + (a_{2} + b_{1})^{2}. Similarly, if Z_{1} ≤ ⋯ ≤ Z_{n^{2}} then (Z_{1} + ⋯ + Z_{n})^{2} + ⋯ + (Z_{n(n−1)+1} + ⋯ + Z_{n(n−1)+n})^{2} ≥ (Z_{π_{1}} + ⋯ + Z_{π_{n}})^{2} + ⋯ + (Z_{π_{n(n−1)+1}} + ⋯ + Z_{π_{n(n−1)+n}})^{2} for all permutation π = {π_{1}, …, π_{n}} ∈ P. This is useful to show that

∑_{ij}(Y_{ij} − Y_{i·})^{2} = ∑_{ij} Y_{ij}^{2} − n ∑_{i} Y_{i·}^{2}
≤ ∑_{ij} Y_{ij}^{2} − n ∑_{j} Y_{·j}^{2}
= ∑_{ij}(Y_{ij} − Y_{·j})^{2}.

Finally,

\frac{1}{n − 1} ∑_{i j}(Y_{ij} − Y_{i·} − Y_{·j} + Y_{··})^{2} = \frac{1}{n − 1}[∑_{i j}(Y_{ij} − Y_{i·})^{2} + ∑_{i j}(Y_{·j} − Y_{··})^{2}]
= \frac{1}{n} ∑_{i j}(Y_{ij} − Y_{i·})^{2}

−\frac{1}{n(n − 1)}[n ∑_{i j}(Y_{·j} − Y_{··})^{2} − ∑_{i j}(Y_{ij} − Y_{i·})^{2}]

= \frac{1}{n} ∑_{i j}(Y_{ij} − Y_{i·})^{2}

−\frac{1}{n(n − 1)}[∑_{j}\{∑_{i}(Y_{ij} − Y_{i·})\}^{2} − ∑_{i j}(Y_{ij} − Y_{i·})^{2}]

= \frac{1}{n} ∑_{i j}(Y_{ij} − Y_{i·})^{2}

−\frac{1}{n(n − 1)}[∑_{j}\{∑_{i≠k}(Y_{ij} − Y_{i·})(Y_{kj} − Y_{k·})\}]

= \frac{1}{n} ∑_{i j}(Y_{ij} − Y_{i·})^{2}

−\frac{1}{n(n − 1)}[∑_{i≠k}\{∑_{j}(Y_{ij} − Y_{i·})(Y_{kj} − Y_{k·})\}]

≤ \frac{1}{n} ∑_{i j}(Y_{ij} − Y_{i·})^{2}
≤ \frac{1}{n + 1}[∑_{i j}(Y_{ij} − Y_{i·})^{2} + \frac{1}{n} ∑_{i j}(Y_{ij} − Y_{·j})^{2}]
≤ \frac{1}{n + 1}[∑_{i j}(Y_{ij} − Y_{i·})^{2} + \frac{1}{n} ∑_{i j}(Y_{ij} − Y_{·j})^{2}]

+\frac{1}{n(n + 1)} ∑_{i≠j}(∑_{k}(Y_{ki} − Y_{·i})(Y_{kj} − Y_{·j}))

= \frac{1}{n + 1}[∑_{i j}(Y_{ij} − Y_{i·})^{2} + \frac{1}{n} ∑_{i}\{∑_{j}(Y_{ij} − Y_{·j})\}^{2})]
= \frac{1}{n + 1}[∑_{i j}(Y_{ij} − Y_{i·})^{2} + ∑_{i j}(Y_{i·} − Y_{··})^{2}]
= \frac{1}{n + 1} ∑_{i j}(Y_{ij} − Y_{··})^{2}.

We conclude this section with an application of the obtained results to the special case of the uniform distribution because here the variances have very simple forms.

Example 2.1 Consider a U(0, 1) population of size n. Here, σ^{2} = \frac{1}{12},

E[X_{(i)}] = \frac{i}{n^{2} + 1} and Cov(X_{(i)}, X_{(j)}) = \frac{i(n^{2} − j + 1)}{(n^{2} + 1)^{2}(n^{2} + 2)}, 1 ≤ i ≤ j ≤ n^{2},

so

Var(μ̂_{S}) = \frac{1}{n}σ^{2}, Var(μ̂_{R}) = \frac{2}{n(n + 1)}σ^{2}, Var(μ̂_{A}) = \frac{n^{2} + 5}{(n^{2} + 1)(n^{2} + 2)}σ^{2}.

We have already shown, in Theorem 2.1, that \frac{σ^{2}}{n^{2}} is a lower bound for Var(μ̂_{A}). Also, it is shown (see Theorem 3.1) that, \frac{σ^{2}}{n^{2}} is a lower bound for Var(μ̂_{R}) as well. Therefore we can analyse how much do we lose in comparison to the lower bound \frac{σ^{2}}{n^{2}} by using μ̂_{A} or μ̂_{R} in estimating the population mean μ = \frac{1}{2}. For the U(0, 1) distribution

(n^{2}Var(μ̂_{R}) − σ^{2}) = \frac{n − 1}{n + 1} = O(1),

while

(n^{2}Var(μ̂_{A}) − σ^{2}) = 2\frac{(n^{2} − 1)}{(n^{2} + 1)(n^{2} + 2)} = O(\frac{1}{n^{2}} ).

3 A special class of random estimators

In this section we propose a class of random estimators of the mean and develop a general dominance result leading to a sufficient condition for a random estimator μ̂_{1} to dominate another random estimator μ̂_{2}. Our treatment permits consideration of a more general type of estimators than that discussed in the preceding section and includes the estimators μ̂_{S}, μ̂_{R} and μ̂_{A} as special cases. We are concern with random estimators μ̂ that can be written in the form

μ̂ = \frac{1}{n} ∑_{j=1}^{n^{2}} I_{j}X_{(j)},  (5)

where I_{1}, …, I_{n^{2}} are Bernoulli(\frac{1}{n}) random variables, independent of X_{(1)}, …, X_{(n^{2})}, satisfying the regularity condition ∑_{j=1}^{n^{2}} I_{j} = n. In the above representation of the random estimator μ̂ the inclusion probabilities of the first and second orders are given by p_{i} = E[I_{i}] = \frac{1}{n} and p_{ij} = E[I_{i}I_{j}], i, j = 1, …, n^{2}, respectively. In particular, we have Cov(I_{i}, I_{j}) = p_{ij}−\frac{1}{n^{2}} and ∑_{i=1}^{n^{2}} p_{ij} = ∑_{j=1}^{n^{2}} p_{ij} = 1, i, j = 1, …, n^{2}.

Remark 3.1 Let R_{1}, …, R_{n^{2}} be the rank statistics corresponding to the sample X_{1}, …, X_{n^{2}} and assume that the ties are broken randomly. The rank statistics are independent of X_{(1)}, …, X_{(n^{2})}. Therefore, the construction of the random variables I_{1}, …, I_{n^{2}} could involve the rank statistics.

Example 3.1 (Simple Random Sampling). The estimator μ̂_{S} as in (1), corresponding to a simple random sample of size n from the underlying population, is in the form of (5) with Bernoulli random variables given by

I_{j} = \{ \begin{matrix}1 & if j ∈ \{R_{i}: i = 1 …, n\} \\ 0 & otherwise \end{matrix} , j = 1, …, n^{2},

and p_{ij} = \frac{1}{n(n+1)}, 1 ≤ i < j ≤ n^{2}.

Example 3.2 (Ranked Set Sampling). Let {R_{i(k)} : k = 1, …, n} be the order statistics based on {R_{ik}: k = 1, …, n}, i = 1, …, n. The estimator μ̂_{R} as in (2), corresponding to a ranked set sample of size n from the underlying population, is in the form of (5) with Bernoulli random variables given by

I_{j} = \{ \begin{matrix}1 & if j ∈ \{R_{i(i)} : i = 1 …, n\} \\ 0 & otherwise \end{matrix} , j = 1, …, n^{2}.

To obtain the corresponding probability of inclusion of the second order, consider a ranked set sampling algorithm with n = 2. There are six ways to put {X_{(1)}, …, X_{(4)}} in two rows which are,

( \begin{matrix}X_{(1)} & X_{(2)} \\ X_{(3)} & X_{(4)} \end{matrix} ) ( \begin{matrix}X_{(1)} & X_{(3)} \\ X_{(2)} & X_{(4)} \end{matrix} ) ( \begin{matrix}X_{(1)} & X_{(4)} \\ X_{(2)} & X_{(3)} \end{matrix} )

( \begin{matrix}X_{(2)} & X_{(3)} \\ X_{(1)} & X_{(4)} \end{matrix} ) ( \begin{matrix}X_{(2)} & X_{(4)} \\ X_{(1)} & X_{(3)} \end{matrix} ) ( \begin{matrix}X_{(3)} & X_{(4)} \\ X_{(1)} & X_{(2)} \end{matrix} )

leading to the following selections, (X_{(1)}, X_{(4)}), (X_{(1)}, X_{(4)}), (X_{(1)}, X_{(3)}), (X_{(2)}, X_{(4)}), (X_{(2)}, X_{(3)}), (X_{(2)}, X_{(3)}). Therefore, p_{12} = 0, p_{13} = \frac{1}{6}, p_{14} = \frac{1}{3}, p_{23} = \frac{1}{3}, p_{24} = \frac{1}{6} and p_{34} = 0. For n > 2 finding the p_{ij} is a combinatorics problem which leads to the general solution

p_{ij} = \frac{1}{n^{2} − 1}[1 − \frac{A}{\binom{n^{2}−2}{n−1 n−1 n(n−2)}}], 1 ≤ i < j ≤ n^{2},

with

A = ∑_{1≤k≤ℓ≤n} \binom{i − 1}{k − 1}\binom{j − i − 1}{ℓ − k}\binom{n^{2} − j}{n − ℓ}\binom{j − ℓ − 1}{k − 1}\binom{n(n − 1) + ℓ − j}{n − k}.

Example 3.3 (Antithetic Sampling). Let π = {π_{1}, …, π_{n}} be a random permutation on {1, …, n}. The antithetic estimator μ̂_{A} as in (4) is in the form of (5) and the corresponding Bernoulli random variables are given by

I_{j} = \{ \begin{matrix}1 & if j ∈ \{n(i − 1) + π_{i}: i = 1 …, n\} \\ 0 & otherwise \end{matrix} , j = 1, …, n^{2},

with

p_{ij} = \{ \begin{matrix}0 & if ⌈\frac{i}{n}⌉ = ⌈\frac{j}{n}⌉ and i ≢ j (mod n) \\ 0 & if ⌈\frac{i}{n}⌉ ≠ ⌈\frac{j}{n}⌉ and i ≡ j (mod n) \\ \frac{1}{n(n−1)} & if ⌈\frac{i}{n}⌉ ≠ ⌈\frac{j}{n}⌉ and i ≢ j (mod n) \end{matrix} , 1 ≤ i < j ≤ n^{2}.

We pursue with a Theorem and a Corollary regarding to the properties of the random estimators μ̂ as in (5); the latter one giving simple and general sufficient condition for dominance results.

Theorem 3.1 The random estimator μ̂ as in (5) is unbiased, in estimating the population mean μ, with

Var[μ̂] = \frac{1}{n^{2}} \{σ^{2} + ∑_{i=1}^{n^{2}−1} ∑_{j=1}^{n^{2}−1} τ_{ij}E[S_{i}S_{j}]\} ≥ \frac{σ^{2}}{n^{2}} ,

where

τ_{ij} = Cov[∑_{k=1}^{i} I_{k}, ∑_{ℓ=1}^{j} I_{ℓ}] = ∑_{k=1}^{i} ∑_{ℓ=1}^{j}[p_{kℓ} − \frac{1}{n^{2}} ] and S_{i} = X_{(i+1)} − X_{(i)}, i, j = 1, …, n^{2} − 1.

Proof. The expectation of μ̂ is given by

E[μ̂] = E[E[μ̂|X_{(1)}, …, X_{(n^{2})}]]
= E[X̄]
= μ.

To obtain the variance of μ̂, we have

Var[μ̂] = Var[E[μ̂|X_{(1)}, …, X_{(n^{2})}]] + E[Var[μ̂|X_{(1)}, …, X_{(n^{2})}]]
= Var[X̄] + \frac{1}{n^{2}} E[∑_{i=1}^{n^{2}} ∑_{j=1}^{n^{2}}(p_{ij} − \frac{1}{n^{2}} )X_{(i)}X_{(j)}]
= \frac{1}{n^{2}} \{σ^{2} + E[ ∑_{i=1}^{n^{2}−1} ∑_{j=1}^{n^{2}−1} τ_{ij}S_{i}S_{j}]\}
≥ \frac{σ^{2}}{n^{2}} ,

and the inequality holds because Var[μ̂] − \frac{σ^{2}}{n^{2}} = E[Var[μ̂|X_{(1)}, …, X_{(n^{2})}]] ≥ 0. Corollary 3.2 If μ̂ and μ̂^{∗} are two random estimators, in the form of (5), such that τ_{ij} ≥ τ_{ij}^{∗}; for all i, j = 1, …, n^{2} − 1; then Var[μ̂] ≥ Var[μ̂^{∗}]. Proof. The proof is based on the fact that S_{i} ≥ 0 for i = 1, …, n^{2} − 1.

We pursue with various applications of the above results. Example 3.4 In this example we provide a new proof for superiority of μ̂_{A} over μ̂_{S} in estimating the population mean in the spirit of Corollary 3.2. By virtue of Corollary 3.2, we only need to obtain τ_{Aij} and τ_{Sij} corresponding to μ̂_{A} and μ̂_{S} and show that τ_{Aij} ≤ τ_{Sij} for all i, j = 1, …, n^{2} − 1. Let

i_{1} = ⌈\frac{i}{n}⌉, i_{2} = i − n(i_{1} − 1),

j_{1} = ⌈\frac{j}{n}⌉, j_{2} = j − n(j_{1} − 1).

We obtain

τ_{Sij} = \frac{1}{n + 1}(i ∧ j)(1 − \frac{i ∨ j}{n^{2}} ),

τ_{Aij} = \{ \begin{matrix}\frac{1}{n^{2}} (i_{2} ∧ j_{2})(n − (i_{2} ∨ j_{2})) & if & i_{1} = j_{1} \\ −\frac{1}{n^{2}(n−1)}(i_{2} ∧ j_{2})(n − (i_{2} ∨ j_{2})) & if & i_{1} < j_{1} \end{matrix} .

When 1 ≤ i_{1} < j_{1} ≤ n we have τ_{Aij} ≤ 0 ≤ τ_{Sij}. Assume that i_{1} = j_{1}. In this case

τ_{Sij} = \frac{1}{n^{2}(n + 1)}[n(i_{1} − 1) + (i_{2} ∧ j_{2})][n(n − i_{1}) + (n − (i_{2} ∨ j_{2}))].

If i_{1} > (i_{2} ∧ j_{2}) then [n(i_{1} − 1) + (i_{2} ∧ j_{2})] ≥ (n + 1)(i_{2} ∧ j_{2}), and

[n(n − i_{1}) + (n − (i_{2} ∨ j_{2}))] ≥ (n − (i_{2} ∨ j_{2})).

So τ_{Sij} ≥ τ_{Aij}. If i_{1} ≤ (i_{2} ∧ j_{2}) then [n(n − i_{1}) + (n − (i_{2} ∨ j_{2}))] ≥ (n + 1)(n − (i_{2} ∨ j_{2})), and

[n(i_{1} − 1) + (i_{2} ∧ j_{2})] ≥ (i_{2} ∧ j_{2}).

So τ_{Sij} ≥ τ_{Aij}. Therefore, τ_{Sij} ≥ τ_{Aij} for all 1 ≤ i ≤ j ≤ n^{2} − 1.

In the reminding of this section, we discuss the problem of superiority of the antithetic sampling over the ranked set sampling in estimating the population mean. First, the following Example highlights the attractive feature of μ̂_{A} in dominating μ̂_{R} for n = 2. As we see in the next section, this latter property of μ̂_{A} leads us to introduce a new variation of the ranked set sampling and correspondingly a new estimator, μ̂_{V}, which outperforms μ̂_{R}.

Example 3.5 . If n = 2 then

Var[μ̂_{R}|X_{(1)}, …, X_{(4)}] − Var[μ̂_{A}|X_{(1)}, …, X_{(4)}] = \frac{1}{24}[(X_{(2)} − X_{(1)})(X_{(4)} − X_{(3)})] ≥ 0.

In the following two Examples we consider the performance of μ̂_{A} in comparison with μ̂_{R} for n = 3 and n = 4. Let T = (τ_{ij})_{i,j=1,...,n^{2}−1}. Note that for all n ≥ 2 we have

(T_{R} − T_{A})_{ij} = (T_{R} − T_{A})_{ji} = (T_{R} − T_{A})_{n^{2}−i,n^{2}−j}, i, j = 1, …, n^{2} − 1.

Example 3.6 . For n = 3, the comparison between μ̂_{A} and μ̂_{R} leads to

100(T_{R} − T_{A}) = ( \begin{matrix}0 & 0 & 3.57 & 10.12 & 2.38 & −3.57 & 8.33 & 4.17 \\ 0 & 0 & 10.71 & 8.33 & 9.05 & −4.29 & 1.19 & 8.33 \\ 3.57 & 10.71 & 21.43 & 10.71 & 2.86 & −2.14 & −4.29 & −3.57 \\ 10.12 & 8.33 & 10.71 & 0 & 0 & 2.86 & 9.05 & 2.38 \\ 2.38 & 9.05 & 2.86 & 0 & 0 & 10.71 & 8.33 & 10.12 \\ −3.57 & −4.29 & −2.14 & 2.86 & 10.71 & 21.43 & 10.71 & 3.57 \\ 8.33 & 1.19 & −4.29 & 9.05 & 8.33 & 10.71 & 0 & 0 \\ 4.17 & 8.33 & −3.57 & 2.38 & 10.12 & 3.57 & 0 & 0 \end{matrix} ) .

In particular,

∑_{i,j∈{1,6}}(τ_{Rij} − τ_{Aij})S_{i}S_{j} → −∞ when S_{6} = 1 and S_{1} → ∞.

Therefore, at this point, we cannot claim that Var[μ̂_{A}] ≤ Var[μ̂_{R}] for n = 3. On the other hand, for Bernoulli(p) trials, S_{i}S_{j} = 0 if i ≠ j and

Var[μ̂] = \frac{pq}{n^{2}} + ∑_{k=1}^{n^{2}−1} τ_{kk}P[Binomial(n^{2}, p) = n^{2} − k + 1].

We obtain that

\frac{Var[μ̂_{R}] − Var[μ̂_{A}]}{Var[μ̂_{S}]} = \frac{3 × 0.2143}{pq}\{P[Binomial(9, p) = 3] + P[Binomial(9, p) = 6]\}
= 54(pq)^{2}(1 − 3pq).

In this application we see that (pq)^{2}(1 − 3pq) reaches its maximum when p = \frac{1}{3} and p = \frac{2}{3}.

100(T_{R} − T_{A})_{i=1,...,15,j=1,...,8} = ( \begin{matrix}0 & 0 & 1.43 & 4.07 & 7.71 & 3.85 & 0.64 & −2.05 \\ 0 & 0 & 4.29 & 10.33 & 9.56 & 10.09 & 3, 38 & −2.42 \\ 1.43 & 4.29 & 8.57 & 18.57 & 13.61 & 10.17 & 8.10 & −1.12 \\ 4.07 & 10.33 & 18.57 & 28.57 & 19.60 & 12.22 & 6.34 & 1.86 \\ 7.71 & 9.56 & 13.61 & 19.60 & 10.62 & 7.71 & 6.34 & 6.48 \\ 3.85 & 10.09 & 10.17 & 12.22 & 7.71 & 4.80 & 7.99 & 12.73 \\ 0.64 & 3.38 & 8.10 & 6.34 & 6.34 & 7.99 & 11.19 & 20.57 \\ −2.05 & −2.42 & −1.12 & 1.86 & 6.48 & 12.73 & 20.57 & 29.96 \\ 3.97 & 0.84 & −0.97 & −1.36 & 8.06 & 10.68 & 14.84 & 20.57 \\ 1.94 & 4.68 & 0.06 & −3.43 & 2.66 & 10.11 & 10.68 & 12.73 \\ 0.09 & 0.63 & 1.82 & −4.49 & −1.48 & 2.66 & 8.06 & 6.48 \\ −1.65 & −3.10 & −4.17 & −4.68 & −4.49 & −3.43 & −1.36 & 1.86 \\ 5.00 & 1.72 & −1.38 & −4.17 & 1.82 & 6.32 & −0.97 & −1.12 \\ 3.33 & 6.67 & 1.72 & −3.10 & 6.33 & 4.68 & 8.39 & −2.42 \\ 1.67 & 3.33 & 5.00 & −1.65 & 0.09 & 1.94 & 3.97 & −2.05 \end{matrix} )

Here again, τ_{Rii} ≥ τ_{Aii} for i = 1, …, 15. Therefore, Var(μ̂_{R}) ≥ Var(μ̂_{A}) for Bernoulli(p) trials.

4 A variation of ranked set sampling

In this section, following the optimal property of μ̂_{A} in dominating μ̂_{R} for n = 2, explored in Example 3.5, we introduce another variation of the ranked set sampling and correspondingly another estimator of the population mean called μ̂_{V} which outperforms μ̂_{R}. In fact, when n = 2, μ̂_{V} is the antithetic estimator. Let m_{1} = ⌊\frac{n}{2}⌋ and m_{2} = ⌈\frac{n}{2}⌉. To construct μ̂_{V} we start with a useful representation of the ranked set sampling estimator μ̂_{R}. Note that μ̂_{R} in (2) can be written as follows

nμ̂_{R} = ∑_{i=1}^{n} X_{i(i)}
= \{ \begin{matrix}∑_{k=1}^{m_{1}}(X_{k(k)} + X_{k+m_{2} (k+m_{2})}) & if n is even \\ ∑_{k=1}^{m_{1}}(X_{k(k)} + X_{k+m_{2} (k+m_{2})}) + X_{m_{2}(m_{2})} & if n is odd \end{matrix} ,

Or equivalently

μ̂_{R} = \frac{1}{n} ∑_{k=1}^{m_{2}} Z_{k},  (6)

with Z_{k} = (X_{k(k)} + X_{k+m_{2} (k+m_{2})}); k = 1, …, m_{1}. Now, we define our new proposed estimator as

μ̂_{V} = \frac{1}{n} ∑_{k=1}^{m_{2}} V_{k},  (7)

with

V_{k} = \{ \begin{matrix}(X_{k(k)} ∧ X_{k+m_{2} (k)}) + (X_{k (k+m_{2})} ∨ X_{k+m_{2} (k+m_{2})}) with probatility \frac{1}{2} \\ (X_{k(k)} ∨ X_{k+m_{2} (k)}) + (X_{k (k+m_{2})} ∧ X_{k+m_{2} (k+m_{2})}) with probatility \frac{1}{2} \end{matrix} ,

and V_{m_{2}} = X_{m_{2},(m_{2})} if n is even.

We now are ready to pursue with the following key Theorem of this section regarding to the optimal properties of μ̂_{V} such as its superiority over μ̂_{R} in estimating the population mean.

Theorem 4.1 The estimator μ̂_{V} as in (7) is an unbiased estimator of the population mean and Var(μ̂_{V}) ≤ Var(μ̂_{R}).

Proof. First of all, note that

E[μ̂_{R}] = \frac{1}{n} ∑_{k=1}^{m_{2}} E[Z_{k}], Var[μ̂_{R}] = \frac{1}{n^{2}} ∑_{k=1}^{m_{2}} Var[Z_{k}].

Similarly

E[μ̂_{V}] = \frac{1}{n} ∑_{k=1}^{m_{2}} E[V_{k}], Var[μ̂_{V}] = \frac{1}{n^{2}} ∑_{k=1}^{m_{2}} Var[V_{k}].

Let 1 ≤ ℓ_{1} < ℓ_{2} ≤ n, and define

Z = X_{ℓ_{1}(ℓ_{1})} + X_{ℓ_{2}(ℓ_{2})} V = \{ \begin{matrix}X_{ℓ_{1}(ℓ_{1})} ∧ X_{ℓ_{2}(ℓ_{1})} + X_{ℓ_{1}(ℓ_{2})} ∨ X_{ℓ_{2}(ℓ_{2})} with probatility \frac{1}{2} \\ X_{ℓ_{1}(ℓ_{1})} ∨ X_{ℓ_{2}(ℓ_{1})} + X_{ℓ_{1}(ℓ_{2})} ∧ X_{ℓ_{2}(ℓ_{2})} with probatility \frac{1}{2} \end{matrix}

Now, consider the event A_{ijkl}, 1 ≤ i < j < k < l ≤ n^{2}, given by

{R_{ℓ_{1}(ℓ_{1})}, R_{ℓ_{1}(ℓ_{2})}, R_{ℓ_{2}(ℓ_{1})}, R_{ℓ_{2}(ℓ_{2})}} = {i, j, k, l},

and assume that (i, j, k, l) is such that P[A_{ijkl}] > 0. When the event A_{ijkl} occurs it implies that

(R_{ℓ_{1}(ℓ_{1})} ∧ R_{ℓ_{2}(ℓ_{2})}, R_{ℓ_{1}(ℓ_{1})} ∨ R_{ℓ_{2}(ℓ_{2})}) ∈ {(i, l), (j, k), (i, k), (j, l)}.

Since ((R_{ℓ_{1}(ℓ_{1})}, R_{ℓ_{1}(ℓ_{2})}), (R_{ℓ_{2}(ℓ_{1})}, R_{ℓ_{2}(ℓ_{2})})) and ((R_{ℓ_{2}(ℓ_{1})}, R_{ℓ_{2}(ℓ_{2})}), (R_{ℓ_{1}(ℓ_{1})}, R_{ℓ_{1}(ℓ_{2})})) have the same distribution we obtain that P[(R_{ℓ_{1}(ℓ_{1})}∧R_{ℓ_{2}(ℓ_{2})}, R_{ℓ_{1}(ℓ_{1})}∨R_{ℓ_{2}(ℓ_{2})}) = (i, l)|A_{ijkl}] = P[(R_{ℓ_{1}(ℓ_{1})}∧R_{ℓ_{2}(ℓ_{2})}, R_{ℓ_{1}(ℓ_{1})}∨R_{ℓ_{2}(ℓ_{2})}) = (j, k)|A_{ijkl}]. Similarly P[(R_{ℓ_{1}(ℓ_{1})}∧R_{ℓ_{2}(ℓ_{2})}, R_{ℓ_{1}(ℓ_{1})}∨R_{ℓ_{2}(ℓ_{2})}) = (i, k)|A_{ijkl}] = P[(R_{ℓ_{1}(ℓ_{1})}∧R_{ℓ_{2}(ℓ_{2})}, R_{ℓ_{1}(ℓ_{1})}∨R_{ℓ_{2}(ℓ_{2})}) = (j, l)|A_{ijkl}]. Let

α_{ijkl} = P[(R_{ℓ_{1}(ℓ_{1})} ∧ R_{ℓ_{2}(ℓ_{2})}, R_{ℓ_{1}(ℓ_{1})} ∨ R_{ℓ_{2}(ℓ_{2})}) ∈ {(i, k), (j, l)}|A_{ijkl}],
μ_{ijkl} = \frac{X_{(i)} + X_{(j)} + X_{(k)} + X_{(l)}}{2} .

Conditionally on the event A_{ijkl} we obtain that

Z = \{ \begin{matrix}X_{(i)} + X_{(l)} & with probability & \frac{(1−α_{ijkl})}{2} \\ X_{(j)} + X_{(k)} & with probability & \frac{(1−α_{ijkl})}{2} \\ X_{(i)} + X_{(k)} & with probability & \frac{α_{ijkl}}{2} \\ X_{(j)} + X_{(l)} & with probability & \frac{α_{ijkl}}{2} \end{matrix} , V = \{ \begin{matrix}X_{(i)} + X_{(l)} & with probability & \frac{1}{2} \\ X_{(j)} + X_{(k)} & with probability & \frac{1}{2} \end{matrix} .

Now,

E[Z|A_{ijkl}] = μ_{ijkl} = E[V |A_{ijkl}],

establishing the unbiasedness of μ̂_{V}. Finally,

Var[Z|A_{ijkl}] = (X_{(i)} + X_{(l)} − μ_{ijkl})^{2}(1 − α_{ijkl}) + (X_{(i)} + X_{(k)} − μ_{ijkl})^{2}α_{ijkl}
= (X_{(i)} + X_{(l)} − μ_{ijkl})^{2} + \frac{α_{ijkl}}{2} (X_{(j)} − X_{(i)})(X_{(l)} − X_{(k)})
= Var[V |A_{ijkl}] + \frac{α_{ijkl}}{2} (X_{(j)} − X_{(i)})(X_{(l)} − X_{(k)})
≥ Var[V |A_{ijkl}].

5 Simulation study

In this section, we carry out a simulation study to compare numerically the performance of the proposed estimators μ̂_{A} and μ̂_{V} with that of μ̂_{R}. It has already been established, in Theorem (3.1), that \frac{σ^{2}}{n^{2}} is a lower bound for the variances of all these estimators. The aim is to reduce, as much as possible, the gap between the variance of the proposed estimators and the lower bound \frac{σ^{2}}{n^{2}} . For this reason, given a distribution P and a sample size n, the comparison between an estimator μ̂ and μ̂_{R} will be based on the measure R_{n} given by

R_{n}(P, μ̂) = \frac{Var_{P} (μ̂) − \frac{σ^{2}}{n^{2}}}{Var_{P} (μ̂_{R}) − \frac{σ^{2}}{n^{2}}} .

In our simulation study we worked on four underlying distributions: the standard normal distribution N (0, 1), the exponential distribution Exp(1), the uniform distribution U(0, 1), and the Poisson distribution Poisson(1).

Figure 1: Simulated values of R_{n}(P, μ̂) for considered distributions: (a) R_{n}(P, μ̂_{V}) (b) R_{n}(P, μ̂_{A}).

For a given sample size n, n ∈ {2, … 30}, and a given distribution P, Figure 1 displays the simulated values of R_{n}(P, μ̂_{A}) and R_{n}(P, μ̂_{V}) based on one million iterations for the considered distributions. The simulations show that R_{n}(P, μ̂_{V}) varies a little bit with n for small values of n and becomes more stable when n gets larger. Overall, the gab between Var_{P} (μ̂_{R}) and \frac{σ^{2}}{n^{2}} has been reduced by a factor of at least 10% using μ̂_{V}.

On the other hand, the simulations indicate that μ̂_{A} is always the best estimator. The theory tells us that R_{2}(P, μ̂_{A}) ≤ 1 for all P. Moreover, the simulations indicate that R_{n}(P, μ̂_{A}) decreases in n when P is fixed. This behavior has been observed in all of the other simulations that we have conducted so far. Perhaps R_{n}(P, μ̂_{A}) ≤ 1 for all P, n ≥ 1.

References

Al–Nasser, A.D. (2007) L–Ranked Set Sampling: A Generalization Procedure For Robust Visual Sampling. Communications in Statistics: Theory and Methods, 36, 33–43.

Al–Saleh, M.F. and Al–Omari, A. (2002) Multistage Ranked Set Sampling. Journal of Statistical Planning and Inference , 102, 273–286.

Al–Saleh, M.F. and Al–Kadiri, M.A. (2000) Double Ranked Set Sampling. Statistics & Probability Letters, 48, 205–212.

Bhoj, D.S. (1997) New Parametric Ranked Set Sampling. Journal of Applied Statistical Science, 6 , 275–289.

Chen, Z., Bai, Z.D., and Sinha, B.K. (2004) Ranked Set Sampling: Theory and Applications. Lecture Notes in Statistics, 176, Springer–Verlag, New York.

Dell, T.R. (1969) The Theory of Some Applications of Ranked Set Sampling. Ph.D. Thesis, University of Georgia, Athens, GA.

Dell, T.R. and Clutter, J.L. (1972) Ranked Set Sampling With Order Statistics Background. Biometrics, 28, 545–553.

Hossain, S.S. and Muttlak, H.A. (1999) Paired Ranked Set Sampling: A More Efficient Procedure. Environmetrics, 10, 195–212. Hossain, S.S. and Muttlak, H.A. (2001) Selected Ranked Set Sampling. Australian & New Zealand Journal of Statistics, 43, 311–325.

Kaur, A., Patil, G.P., Sinha, A.K., and Taillie, C. (1995) Ranked Set Sampling: An Annotated Bibliography. Environmental and Ecological Statistics, 2, 25–54.

Li., D., Sinha, B.K., and Perron, F. (1999) Random Selection in Ranked Set Sampling and Its Applications. Journal of Statistical Planning and Inference, 76, 185–201.

McIntyre, G.A. (1952) A Method for Unbiased Selective Sampling, Using Ranked Sets. Australian Journal of Agricultural Research, 3, 385–390.

Mutllak, H.A. (1997) Median Ranked Set Sampling. Journal of Applied Statistical Science, 6, 245–255.

Samawi, H., Abu–Daayeh, H.A., and Ahmed, S. (1996) Estimating the Population Mean Using Extreme Ranked Set Sampling. Biometrical Journal, 38, 577–586.

Takahasi, K. and Wakimoto, K. (1968) On Unbiased Estimates of the Population Mean Based on the Sample Stratified by Means of Ordering. Annals of the Institute of Statistical Mathematics, 20, 1–31.