\documentstyle[12pt,draft]{article}


%__________________LAYOUT____________________________
\oddsidemargin 0cm
\textwidth 15.3cm
\textheight 24.2cm
\topmargin -1.5cm





\def \R{I\!\!R}                                 %space
\def \N{I\!\!N}                                 %the set of natural numbers
\def \X{X\!\!\!\!\!\!X}                         %sample
\def \Y{Y\!\!\!\!\!\!Y}                         %sample
\def \Z{Z\!\!\!Z}                               %the set of whole numbers
\def \F{I\!\!F}
\def \SIGMA{{\bf\Sigma}}
\def \nn{\mbox{\boldmath$\nu$}}
\def \mm{\mbox{\boldmath$\mu$}}
\begin{document}

\author{Mikhail Nikulin, Vassilly Voinov\\
\medskip
{\Victor Segalen Bordeaux 2 University, Bordeaux, France\\
Kazakhstan Institute of Management, Economics and Strategic
Research, Almaty, Kazakhstan}}

\centerline{\large \bf Unbiased Estimation of Reliability in Stress-Strength Model}
\centerline{\large \bf and Similar Problems}
\vspace{1cm}

\leftline{\scriptsize MIKHAIL S. NIKULIN }


\leftline{\scriptsize \it  Steklov Mathematical Institute, 27 Fontanka, 191011,
Saint Petersburg, Russia  }
\bigskip

\leftline{\scriptsize VASSILY G. VOINOV}

\leftline{\scriptsize \it Kazakhstan Institute of Management, Economics and Strategic
Research, Almaty, Kazakhstan}
\vspace{0.5cm}


\noindent {\bf Abstract:} Some problems related to an application
of the unbiased estimators for sampling inspection, lifetime and
reliability testing are discussed. Continuous and discrete,
univariate and multivariate parametric probability models are
considered.

\bigskip

\noindent {\bf Keywords and phrases:} unbiased estimation, life
distributions, multivariate failure rates, multivariate polynomial
laws, sampling inspection, reliability and quality of products, stress-strength 
model

\vspace*{30pt}


\leftline{\bf 1. Introduction}

Problems of estimating probabilities like ${\bf P}({\bf X}<{\bf
Y})$, ${\bf P}({\bf X}<{\bf a})$, etc., ${\bf X}$ and ${\bf Y}$
being random vectors and ${\bf a}$ being a given threshold vector,
are often encountered in many applications of statistics. Such
problems arise most often in assessing the lifetime and
reliability (Basu (1988), Voinov and Nikulin (1993, 1996)) and testing the quality of
products (Belyaev(1975)). Among other applications it is worthy to
mention a probabilistic description of radioactively polluted
lands (Nikulin et al(1999)), sampling surveying of populations
(Nikulin, Price, Turetayev, Voinov(1999)) and problems of texture
satellite images recognition.

The problem of estimating the reliability of elements and systems
is one of the main problems in engineering. Suppose that some
stress is applied to an element or to a whole product. A product
is destroyed if at some moment of time the stress applied exceeds
the strength limit of it. The stress is a function of the
equipment on which this element or system is tested and its value
at each moment of time can be described by a random variable. The
strength limit of the product is defined by the stress at which
destruction occurs. It depends on material properties, technology,
etc. If a few randomly chosen products being in mass production
are tested then their strength limit $x$ can also be described by
a random variable $X$ independent on $Y$ describing the stress.
The strength limit $x$ and the maximum stress $y$ are realizations
of the random variables $X$ and $Y$, respectively, and therefore
the probability ${\bf P}\{Y<X\}$ that the maximum stress does not
exceed the strength limit during the testing time $t\in [0,T]$ is
considered to be a system reliability (see, e.g., Chandra and Owen
(1975)). Naturally, products are not destroyed in every test but
since this takes place, it is necessary to be able to get good
estimators of the probability ${\bf P}\{Y<X\}$. Sometimes a
supposition that components of a system fail independently of each
other is inappropriate (see e.g. Lindley and Singpurwalla (1986)).
In such cases one should consider multivariate models with
correlated components. In this case probabilities of interest look
like ${\bf P}({\bf Y}<{\bf X})$, where ${\bf Y}$ and ${\bf X}$ are
random vectors. Different techniques are used for estimating such
probabilities. Maximum likelihood estimators as well as method of
moments estimators seem to be the most widely used. Nonparametric
estimators also have interesting applications (see, e.g., Bakker
(1990), Van der Laan (1991,1994)). Nevertheless, if there exist
the complete sufficient statistics and parametric functions of
interest are estimable, one may use the minimum variance unbiased
estimators. It is important, first, because reliability tests are
rather often made for several samples obtained independently and,
based on these samples, a combined estimator (usually the
arithmetic mean one) is constructed. If the estimators used are
biased, the averaged estimator may have a bias greater than its
standard error and, hence, may cause overestimation or
underestimation of the reliability of products. Secondly, a biased
estimator can inadmissibly distort the information on the
reliability of the system consisting of several stochastically
independent subsystems, since the reliability estimator of the
whole system is equal to the product of the reliabilities of its
subsystems (Pugh (1963)).  Many papers are devoted to the problem
of reliability estimation in the case where random variables are
independent identically or nonidentically distributed (see, e.g.,
Arnold (1975), Johnson (1988), Selvavel (1989), McCool (1991)).

For an inspection of products a random sampling is usually used.
In view of this, a selection of an appropriate probability model
is very important. Continuous and discrete, univariate and
multivariate parametric models are frequently used. If
measurements are of high precision, one may use continuous models.
Sometimes one is unable to measure features of products with high
accuracy. In such situations it is much better to quantize the
measurements and use discrete probability models. To assess the
quality of products for both continuous and discrete models,
probabilistic measures like ${\bf P}({\bf X}<{\bf a})$ or ${\bf
P}({\bf X}>{\bf a})$ are usually used. In this expression ${\bf
X}$ stands for a random vector describing properties of a product
and ${\bf a}$ is a prescribed threshold vector. Unbiased
estimators of such probabilities for different univariate and
multivariate, discrete and continuous models can be found in
Voinov and Nikulin (1993,1996). That is why, in this paper we
shall consider only a multivariate discrete probability model
which relates to those probabilities and proves to be also useful
in radioecology and sampling surveying. \vspace*{30pt}

\vspace{0.5cm}

\leftline{\bf 2. Univariate continuous probability models}
\vspace{0.5cm}

Let $\X=(X_1,...,X_n)$ and $\Y=(Y_1,...,Y_m)$ be samples, where $$
X\sim f(x;\theta),\mbox{    }\theta\in \Theta,\mbox{    },Y\sim
g(y;\lambda),\mbox{    }\lambda\in \Lambda.$$ In this case $$ {\bf
P}\{Y<X\}=F(\theta,\lambda).$$ If the parameters $\theta$ and
$\lambda$ are unknown, the problem of the statistical estimation
of the function $F(\theta,\lambda)$ arises.

Let, for example, $ X_i\sim \frac{1}{\sigma}\varphi\left(\frac
{x-\mu}{\sigma}\right),\mbox{    }i=1,2,...,n,$ and $Y\sim
\varphi(y),$ i.e. the parameters of the maximum stress $Y$
distribution are known, $$
\varphi(x)=\frac{1}{\sqrt{2\pi}}e^{-{x^2}/2}.$$ In this case ${\bf
P}\{Y<X\}=\Phi(\mu/\sqrt{1+\sigma^2}).$ The problem is to
construct the best unbiased estimator $\hat{P_1}$ of the function
$\Phi(\mu/\sqrt{1+\sigma^2})$. Using the Rao-Blackwell-Kolmogorov
theorem, Mazumdar (1970) and Downton (1973) obtained this
estimator for the case $n\geq 3$ as:$$
\hat{P}_1=\frac{1}{B\left(\frac{1}{2},\frac{n-2}{2}\right)}\int\limits_{-1} ^{1}
\Phi\left(\bar{X}_n+\frac{uS_n(n-1)}{\sqrt{n}}\right)(1-u^2)^{\frac{n-4}{2} }du,$$
where $\bar{X}_n=\frac{1}{n}\sum\limits_{i=1}^{n}X_i$ and
$S_n^2=\frac{1}{n-1}\sum\limits_{i=1}^{n}(X_i-\bar{X}_n)^2$ are
statistics sufficient for $\mu$ and $\sigma^2$.

Solving the unbiasedness equation for the estimator of the
function \\${\bf P}\{Y<X\}=\Phi(\mu/\sqrt{1+\sigma^2})$, it can be
shown (Voinov (1984)) that there exists the alternative expression
for $\hat{P}_1$ being valid for $n\geq 2$ and equal to
$$\hat{P}_1=\frac{1}{2}+\frac{1}{2\sqrt{\pi}}\sum\limits_{m=0}^{\infty}
\frac{(-1)^m(\bar{X}_n/\sqrt{2})^{2m+1}}{m!(m+1/2)}\times$$
$$\times {}_2{F}_1\left(-m,-m-\frac{1}{2};
\frac{n-1}{2};\frac{(n-1)^2S_n^2}{n\bar{X}_n}\right),$$ where
${}_2F_1(\alpha,\beta;\gamma;x)$ is the Gauss hypergeometrical
function.

Let now parameters of the distribution of $Y$ be unknown and
$$X_i\sim \frac{1}{{\sigma}_x}\varphi\left(\frac
{x-{\mu}_x}{{\sigma}_x}\right),\mbox{      }Y_i\sim
\frac{1}{{\sigma}_y}\varphi\left(\frac
{y-{\mu}_y}{{\sigma}_y}\right).$$ In this case $${\bf P}\{Y<X\}=
\Phi\left(\frac{\mu_x-\mu_y}{\sqrt{\sigma_x^2+\sigma_y^2}}\right).$$
The problem of estimating the probability ${\bf P}\{Y<X\}$, based
on samples $\X=(X_1,...,X_{n_1})$ and $\Y=(Y_1,...,Y_{n_2})$ is
much more realistic than the previous one where the parameters of
the distribution of $Y$ were supposed to be known. In Voinov
(1984) it has been shown that for $n_1\geq 3$ and $n_2\geq 2$ the
uniformly minimum variance unbiased estimator (UMVUE) $\hat{P_2}$
of the function
$\Phi\left((\mu_x-\mu_y)/\sqrt{\sigma_x^2+\sigma_y^2}\right)$ is$$
\hat{P}_2=\left\{\begin{array}{lr}
  1-\int\limits_{-1}^{1}\psi(\rho_1)(1-u^2)^{\frac{n_1-4}{2}}du, & a>1, \\
  \int\limits_{-1}^{1}\psi(\rho_2)(1-u^2)^{\frac{n_1-4}{2}}du & a<-1, \\
  \frac{1}{2}+\frac{\Gamma\left(\frac{n_1-1}{2}\right)}
  {2\sqrt{\pi}\Gamma\left(\frac{n_1-2}{2}\right)}
  \int\limits_{-a}^{a}(1-u^2)^{\frac{n_1-4}{2}}du &  \\
  +\int\limits_{a}^{1}\psi(\rho_3)(1-u^2)^{\frac{n_1-4}{2}}du-
  \int\limits_{a}^{1}\psi(\rho_1)(1-u^2)^{\frac{n_1-4}{2}}du & |a|\leq 1,
\end{array}\right.$$
where $\rho_1=b(a+u)$, $\rho_2=-\rho_1$, $\rho_3=b(u-a)$,$$
a=\frac{\sqrt{n_1}(\bar{X}-\bar{Y})}{S(x)(n_1-1)},\mbox{        }
b=\frac{ \sqrt{n_2}(n_1-1)S(x)}{\sqrt{n_1}(n_2-1)S(y)},$$ $$
S^2(x)=\frac{1}{n_1-1}\sum\limits_{i=1}^{n_1}(X_i-\bar{X})^2,
\mbox{
}S^2(y)=\frac{1}{n_2-1}\sum\limits_{i=1}^{n_2}(Y_i-\bar{Y})^2,$$
$$\psi(\rho_i)=\left\{\begin{array}{lr}
  \frac{\Gamma\left(\frac{n_1-1}{2}\right)\Gamma\left(\frac{n_2-1}{2}\right)
  \rho_i(1-\rho_i)^{\frac{n_2-2}{2}}}{2\pi \Gamma\left(\frac{n_1-2}{2}\right)
  \Gamma\left(\frac{n_2}{2}\right)} &  \\
  \times_2F_{1}\left(1,\frac{n_2-1}{2};\frac{n_2}{2};1-\rho_i^2\right), & 0<\rho_i\leq{1}, \\
  0, & \rho_i>1,\mbox{    }i=1,2,3.
\end{array}\right.$$

Consider the problem of determining the reliability of a system
consisting of $N$ elements with different strength limits, which
in tests is subjected to a stress during the time $T$. The system
reliability is determined as the probability $q={\bf
P}(X<Y_1,X<Y_2,...,X<Y_N)$, where random variables $Y_1,...,Y_N$
correspond to strength limits and $X$ is a random stress. Let all
variables be stochastically independent, with $X\sim
\frac{1}{{\sigma}}\varphi\left(\frac {x-{\mu}}{{\sigma}}\right),$
and $Y_i\mbox{    }(i=1,...,N)$ having the p.d.f. $$
f(y;\alpha)=\left\{\begin{array}{lr}
  \alpha e^{-\alpha y}, & y\geq 0,\mbox{    }\alpha>0, \\
  0 & \mbox{otherwise}.
\end{array}\right.$$The parameter $\alpha$ is considered to be
known. In this case, the probability $q$ is equal (Chandra and
Owen (1975)) to $q_1+q_2$, where $$q_1=\Phi(-\mu/\sigma)\mbox{ and
}q_2=\exp\left\{-\alpha
N\mu+\frac{{\alpha}^2N^2{\sigma}^2}{2}\right\}\left[1+\Phi\left(\alpha
N\sigma-\frac{\mu}{\sigma}\right)\right].$$ If the parameter
$\sigma^2$ is known, then, based on a sample $\X=(X_1,...,X_n)$,
one can construct the UMVUE of the function $q$, being equal to
$\hat{q}(\sigma)= \hat{q}_1(\sigma)+ \hat{q}_2(\sigma)$, where
(see Chandra and Owen (1975)) $ \hat{q}_1(\sigma)=\Phi(-
\sqrt{n}\bar{X}_n/\sigma \sqrt{n-1})$, $$
\bar{q}_2(\sigma)=\exp\{-\alpha N\bar{X}_n+\alpha^2
N^2(n-1)\sigma^2/2n\}\times$$ $$\times\left[1-\Phi\left(\frac{
\sqrt{n-1}\alpha\sigma N}{ \sqrt{n}}-\frac{
\sqrt{n}\bar{X}_n}{\sigma \sqrt{n-1}}\right)\right].$$ Consider a
more typical case when the parameter $\sigma^2$ is unknown. As has
been shown by Lieberman and Resnikoff (1955) the UMVUE $\hat{q}_1$
of the function $q_1=\Phi(-\mu/\sigma)$ is equal to $$
\hat{q}_1=\left\{\begin{array}{lr}
  0, & z<0, \\
  I_z\left(\frac{n-2}{2},\frac{n-2}{2}\right), & 0\leq{z}\leq 1, \\
  1, & z>1,
\end{array}\right.$$
where $$z=\frac{1}{2}\left(1-\frac{\bar{X}_n
\sqrt{n}}{(n-1)S_{n}}\right)$$ and $I_{z}(\alpha,\beta)$ is the
incomplete $\beta$-function. To derive the UMVUE of the function
$\hat{q}_{2}(\sigma)$ with respect to the unknown parameter
$\sigma^2$, we represent $\hat{q}_{2}(\sigma)$ as
$g(\sigma)\exp\{-\alpha N\bar{X}_n\}$, where
$$g(\sigma)=e^{\frac{\alpha^{2}N^2(n-1)\sigma^2}{2n}}\left[1-
\Phi\left(\frac{ \sqrt{n-1}\alpha\sigma N}{ \sqrt{n}}- \frac{
\sqrt{n}\bar{X}_n}{\sigma \sqrt{n-1}}\right)\right]$$
$$=e^{\frac{\alpha^{2}N^2(n-1)\sigma^2}{2n}}-\frac{1}{\sigma}
\int\limits_{-\frac{ \sqrt{n}\bar{X}_n}{ \sqrt{n-1}}}^{\infty}
\varphi\left(\frac{y}{\sigma}\right)\exp\left(-\frac{\sqrt{n-1}\alpha
N}{\sqrt{n}}y\right)dy.$$ According to the known result of Neyman
and Scott (1960) the UMVUE $T_1$ of the function
$\exp\left(\frac{\alpha^{2}N^2(n-1)\sigma^2}{2n}\right)$ is the
statistic
$$T_1=\frac{\Gamma\left(\frac{n-1}{2}\right)}{\left[\frac{T(n-1)^2
\alpha^2
N^2}{4n}\right]^{\frac{n-3}{4}}}I_{\frac{n-3}{2}}\left(\frac{\alpha
N(n-1)S_n}{ \sqrt{n}}\right).$$ Solving the equation of
unbiasedness $ET_2=\frac{1}{\sigma}\varphi(\frac{y}{\sigma})$
w.r.t. $T_2$, it can be shown (Voinov and Nikulin (1993)) that the
UMVUE $\hat{q}$ of $q$ is $\hat{q}_1+T_3$, where
$$T_3=\left\{\begin{array}{lr}
  T_1e^{-\alpha N\bar{X}_n}, &\bar{X}_n\leq-\frac{S_{n}(n-1)}{\sqrt{n}}  \\
  T_1e^{-\alpha N\bar{X}_n}-\frac{\Gamma\left(\frac{n-1}{2}\right)
  {e^{-\alpha N\bar{X}_n}}}{\sqrt{\pi}\Gamma\left(\frac{n-2}{2}\right)}
  \int\limits_{-\frac{ \sqrt{n}\bar{X}_n}{S_{n}(n-1)}}^{1}(1-x^2)^{\frac{n-4}{2}} &  \\
  \times\exp\left\{-\frac{\alpha NS_n(n-1)x}{ \sqrt{n}}\right\}dx, & |\bar{X}_n|<
  \frac{S_{n}(n-1)}{\sqrt{n}}, \\
  0, &\bar{X}_n\geq\frac{S_{n}(n-1)}{ \sqrt{n}}.
\end{array}\right.$$

Let ${\X}=(X_{1},X_{2},...,X_{m})^{T}$ and
$\Y=(Y_{1},Y_{2},...,Y_{n})^{T}$ be two independent random samples
from families with probability density functions
$f_{1}(x;\theta_{1})$ and $f_{2}(y;\theta_{2})$ respectively,
where
\begin{equation}
    f_{1}(x;\theta_{1})=q_{1}(\theta_{1})h_{1}(x),  \hspace{.5in}
    a<x<\theta_{1}<b,
\label{a121}
\end{equation}
and
\begin{equation}
    f_{2}(y;\theta_{2})=q_{2}(\theta_{2})h_{2}(y), \hspace{.5in}
    a<y<\theta_{2}<b.
\label{a122}
\end{equation}
It is supposed that $-\infty<a<b<\infty$ are known constants;
$h_{1}(x), h_{2}(y)$ are positive continuous functions and
$q_{1}(\theta_{1}), q_{2}(\theta_{2})$ are differentiable, where
$$\frac{1}{q_{i}(\theta)}=\int\limits_{a}^{\theta} h_{i}(x)dx,
\hspace{.5in} i=1,2.$$ Let $X_{(1)}<X_{(2)}<...<X_{(m)}$ and $
Y_{(1)}<Y_{(2)}<...<Y_{(n)}$ be sets of order statistics
from~(\ref{a121}) and~(\ref{a122}). The complete sufficient
statistics for parameters $\theta_{1}$ and $\theta_{2}$ are known
to be $X_{(m)}$ and $X_{(n)}$ respectively. Rohatgi (1989) solving
the unbiasedness equation showed that the UMVUE $\varphi
(X_{(m)},Y_{(n)})$ of an estimable function
$g(\theta_{1},\theta_{2})$ is
\begin{equation}
\begin{array}{lc}
   \varphi
   (X_{(m)},Y_{(n)})=g(X_{(m)},Y_{(n)})+\frac{g_{1}(X_{(m)},Y_{(n)})}
   {mq_{1}(X_{(m)})h_{1}(X_{(m)})}+
   \\ \\
   +\frac{g_{2}(X_{(m)},Y_{(n)})}{nq_{2}(Y_{(n)})h_{2}(Y_{(n)})}+
   \frac{g_{12}(X_{(m)},Y_{(n)})}{mnq_{1}(X_{(m)})q_{2}(Y_{(n)})h_{1}
   (X_{(m)})h_{2}(Y_{(n)})},
\end{array}
\label{a123}
\end{equation}
where $$g_{1}(x,y)=\frac{\partial g(x,y)}{\partial
x},\hspace{.2in} g_{2}(x,y)=\frac{\partial g(x,y)}{\partial
y},\hspace{.2in} g_{12}(x,y)=\frac{\partial^{2} g(x,y)}{\partial
x\partial y}.$$ Suppose the samples $\X$ and $\Y$ are Type $II$
censored. Another words one observes only
$X_{(1)},X_{(2)},...,X_{(r_{1})},\hspace{.1in} 1\leq r_{1}\leq m,
$ and $Y_{(1)},Y_{(2)},...,Y_{(r_{2})},\hspace{.1in} 1\leq
r_{2}\leq n,$ respectively. In this case the UMVUE $\varphi
(X_{(r_{1})},Y_{(r_{2})})$ of an estimable function
$g(\theta_{1},\theta_{2})$ is (Selvavel (1989))
\begin{equation}
\begin{array}{lc}
   \varphi(X_{(r_{1})},Y_{(r_{2})})=\frac{(r_{1}-1)!(r_{2}-1)!}{m!n!}
   q_{1}^{r_{1}-1}(X_{(r_{1})})q_{2}^{r_{2}-1}(Y_{(r_{2})})\times
   \\ \\
   \times {\cal
   D}_{11}^{m-r_{1}+1}\left\{\frac{1}{q_{1}^{m}(X_{(r_{1})})}{\cal
   D}_{01}^{n-m-r_{2}+r_{1}}\left(\frac{g(X_{(r_{1})},Y_{(r_{2})})}{q_{2}^{n
   }(Y_{(r_{2})})}\right)\right\},
\end{array}
\label{a124}
\end{equation}
if $n-m-r_{2}+r_{1}\geq 0,$ and
\begin{equation}
\begin{array}{lc}
   \varphi(X_{(r_{1})},Y_{(r_{2})})=\frac{(r_{1}-1)!(r_{2}-1)!}{m!n!}
   q_{1}^{r_{1}-1}(X_{(r_{1})})q_{2}^{r_{2}-1}(Y_{(r_{2})})\times
   \\ \\
   \times {\cal
   D}_{11}^{m-r_{2}+1}\left\{\frac{1}{q_{2}^{n}(Y_{(r_{2})})}{\cal
   D}_{10}^{m-n-r_{1}+r_{2}}\left(\frac{g(X_{(r_{1})},Y_{(r_{2})})}{q_{1}^{m}
   (X_{(r_{1})})}\right)\right\},
\end{array}
\label{a125}
\end{equation}
if $n-m-r_{2}+r_{1}<0,$ where $${\cal D}_{ij}={\cal D}_{ij}\omega
(\theta_{1},\theta_{2})=\frac{1}{h_{1}^{i}(\theta_{1})h_{2}^{j}(\theta_{2})}
\frac{\partial^{i+j}\omega(\theta_{1},\theta_{2})}{\partial^{i}\theta_{1}
\partial^{j}\theta_{2}}.$$
Since
\begin{equation}
  {\bf P}\{X<Y\}=\left\{ \begin{array}{lr}
        q_{1}(\theta_{1})q_{2}(\theta_{2})\int\limits_{a}^{\theta_{2}}
        \left(\int\limits_{x}^{\theta_{2}}h_{2}(y)dy\right)h_{1}(x)dx,&
        \theta_{1}>\theta_{2}, \\
         q_{1}(\theta_{1})q_{2}(\theta_{2})\int\limits_{a}^{\theta_{1}}
         \left(\int\limits_{x}^{\theta_{2}}h_{2}(y)dy\right)h_{1}(x)dx,&
       \theta_{1}<\theta_{2},
        \end{array}
        \right.
\label{a126}
\end{equation}
from~(\ref{a124}),~(\ref{a125}) and~(\ref{a126}) it follows that
the UMVUE $\psi (X_{r_{1}},Y_{r_{2}})$ of ${\bf P}\{X<Y\}$ if
exists is (Selvavel (1989))
\begin{equation}
\begin{array}{lr}
     \psi
     (X_{r_{1}},Y_{r_{2}})=\frac{(r_{1}-1)!(r_{2}-1)!}{m!n!}q_{1}^{r_{1}-1}
     (X_{(r_{1})})q_{2}^{r_{2}-1}(X_{(r_{2})})\times
     \\ \\
     \times {\cal
     D}_{11}^{m-r_{1}+1}\left[\frac{1}{q_{1}^{m}(X_{(r_{1})})}{\cal
D}_{01}^{n-m-r_{2}+r_{1}}\left(\frac{r(X_{(r_{1})},Y_{(r_{2})})}{q_{2}^{n}
     (Y_{r_{2}})}\right)\right],
\end{array}
\label{a127}
\end{equation}
if $X_{(r_{1})}>Y_{(r_{2})}$, and
\begin{equation}
\begin{array}{lcr}
  \psi
     (X_{r_{1}},Y_{r_{2}})=\frac{(r_{1}-1)!(r_{2}-1)!}{m!n(m+r_{2}-r_{1}-1)}
     q_{1}^{r_{1}-1}(X_{(r_{1})})q_{2}^{r_{2}-1}(Y_{(r_{2})})\times
     \\ \\
  \times {\cal
D}_{11}^{m-r_{1}+1}\left[\frac{1}{q_{1}^{m}(X_{(r_{1})})q_{2}^{m+r_{2}-r_{1}}
     (Y_{(r_{2})})}\right.\times
     \\ \\
\left.\times\left\{r(X_{(r_{1})},Y_{(r_{2})})+\frac{(n-m-r_{2}+r_{1})}{(m+r _{2}-r_{1})}
     \right\}\right],
\end{array}
\label{a128}
\end{equation}
if $X_{(r_{1})}<Y_{(r_{2})}$.


\vspace{0.5cm}



\leftline{\bf 3. Multivariate continuous probability models}
\vspace{0.5cm}

Let random variables $X$ and $Y$ follow the absolutely continuous
bivariate exponential distribution of Block and Basu (1974). This
distribution is widely used for describing both physical and
biological systems.

The joint survival functions of this distribution is
\begin{equation}
\begin{array}{lr}
  \bar{F}(x,y)=\frac{\lambda}{\lambda_{1}+\lambda_{2}}\exp
  \{-\lambda_{1}x-\lambda_{2}y-\lambda_{12}\max(x,y)\}-
   \\ \\
  -\frac{\lambda_{12}}{\lambda_{1}+\lambda_{2}}\exp
  \{-\lambda\max(x,y)\}, \hspace{.3in} x,y>0,
\end{array}
\label{a129}
\end{equation}
where $\lambda=\lambda_{1}+\lambda_{2}+\lambda_{12}$. From now on
we shall follow the paper of Klein and Basu (1985), see also
Basu(1988).

Let marginals of~(\ref{a129}) be identical, i.e.
$\lambda_{1}=\lambda_{2}=\alpha$. Denote $\lambda_{12}=\beta$. Let
${\X}=({\bf X}_{1},...,{\bf X}_{n})$ be a sample from~(\ref{a129})
with ${\bf X}_{i}=(X_{i},Y_{i})^{T}, i=1,2,...,n.$ The complete
sufficient for parameters $\lambda$ and $\beta$ statistics are
\\

 $T=\sum\limits_{i=1}^{n}\min(X_{i},Y_{i})$ and
$V=\sum\limits_{i=1}^{n}\max(X_{i},Y_{i})-\sum\limits_{i=1}^{n}\min(X_{i},Y _{i}).$
 \\
The conditional density function of $X_{1}$ and $Y_{1}$ given
$T=t$ and $V=v$ can be derived as
\begin{equation}
   p(x_{1},y_{1}\mid t,v)=\left\{ \begin{array}{llll}
       \frac{(n-1)^{2}(t-x_{1})^{n-2}(v-y_{1}+x_{1})^{n-2}}{2t^{n-1}v^{n-1}},
       \\
       \hspace{.2in} 0<x_{1}<t, \hspace{.1in} x_{1}<y_{1}<v+x_{1};
       \\ \\
\frac{(n-1)^{2}(t-y_{1})^{n-2}(v-x_{1}+y_{1})^{n-2}}{2t^{n-1}v^{n-1}},
        \\
       \hspace{.2in} 0<y_{1}<t, \hspace{.1in} y_{1}<x_{1}<v+y_{1}.
       \end{array}
   \right.
\label{a1211}
\end{equation}

Based on an unique observation ${\bf X}_{1}=(X_{1},Y_{1})^{T}$ the
statistic
\begin{equation}
     \varphi (X_{1},Y_{1})=\left\{ \begin{array}{lc}
         1, & X_{1}>x, \hspace{.1in} Y_{1}>y,
         \\
         0, & \mbox{otherwise},
        \end{array}
   \right.
\label{a1212}
\end{equation}
is an unbiased estimator of $\bar{F}(x,y)={\bf P}\{X>x,Y>y\}.$
Hence, the conditional expectation of $ \varphi (X_{1},Y_{1})$
given $T$ and $Y$ gives the UMVUE $$\psi (T,V)={\bf E}\{\varphi
(X_{1},Y_{1})\mid T,V\}$$ of $\bar{F}(x,y).$ There are three
different cases.

${\bf i)}\hspace{.2in} T>x=y>0.$
$$\psi(T,V)=\left(1-\frac{x}{T}\right)^{n-1}.$$

${\bf ii)}\hspace{.2in} x<y<T.$

$$\psi(T,V)=(1-\frac{y}{T})^{n-1}+\frac{n-1}{2V^{n-1}T^{n-1}}\times$$

$$
\times\sum\limits_{k=0}^{n-1}(-1)^{k}\frac{{n-1\choose{k}}(V-y+T)^{n-k-1}}
 {n+k-1}[(T-x)^{n+k-1}-(T-y)^{n+k-1}].$$

${\bf iii)}\hspace{.2in} T>x>y>0.$ $$
\psi(T,V)=(1-\frac{x}{T})^{n-1}+\frac{n-1}{2V^{n-1}T^{n-1}}\times$$
$$\times\sum\limits_{k=0}^{n-1}(-1)^{k}\frac{{n-1\choose{k}}(V-x+T)^{n-k-1}}
 {n+k-1}[(T-y)^{n+k-1}-(T-x)^{n+k-1}].$$

Let ${\bf X}=(X_{1},...,X_{p})^{T}$ and ${\bf
Y}=(Y_{1},...,Y_{q})^{T}$ be two random vectors such that the
distribution of ${\bf Z}=({\bf X}^{T},{\bf Y}^{T})^{T}$ is
multivariate normal $N_{p+q}(\mm,{\SIGMA}),$ where $
  {\mm}=({\mm}_{x}^{T},{\mm}_{y}^{T})^{T}=(({\mu}_{x_{1}},...,
  {\mu}_{x_{p}})^{T},({\mu}_{y_{1}},...,{\mu}_{y_{q}})^{T})^{T}
$ and the covariance matrix $$
 {\SIGMA}=\left( \begin{array}{cc}
       {\SIGMA_{xx}} & {\SIGMA_{xx}}
       \\
       {\SIGMA_{yx}} & {\SIGMA_{yy}}
       \end{array}
       \right)
$$ is $(p+q)\times(p+q)$ positive definite matrix.

Let ${\bf a}=(a_{1},...,a_{p})^{T}$ and ${\bf
b}=(b_{1},...,b_{q})^{T}$ be two known vectors. The problem is to
estimate $ R={\bf P}\{{\bf a}^{T}{\bf X}>{\bf b}^{T}{\bf Y}\}. $

This problem arises, e.g.,(see Gupta and Gupta (1990)) in systems
where the energy is supplied to systems by $p$ sources and is
consumed through $q$ sources. It is also useful in biochemistry of
brain.

The distribution function of a random variable $U={\bf b}^{T}{\bf
y}-{\bf a}^{T}{\bf X}$ is $\Phi(\frac{u-\mu}{\sigma}),$ where
$\mu={\bf b}^{T}{\mm}_{y}-{\bf a}^{T}{\mm}_{x}$ and
$\sigma^{2}={\bf a}^{T}\SIGMA_{xx}{\bf a}-2{\bf
a}^{T}\SIGMA_{xy}{\bf b}+{\bf b}^{T}\SIGMA_{yy}{\bf b}.$ From this
it follows that $$
 R={\bf P}\{{\bf a}^{T}{\bf X}>{\bf b}^{T}{\bf Y}\}={\bf
 P}\{U<0\}=\Phi(-\frac{\mu}{\sigma}).$$

Let ${\Z}=({\bf Z}_{1},...,{\bf Z}_{n})$ be a random sample of
size $n$ from $${\bf Z}=({\bf X}^{T},{\bf Y}^{T})^{T} \sim
N_{p+q}({\mm},{\SIGMA}),$$ where ${\bf Z}_{i}=({\bf
X}_{i}^{T},{\bf Y}_{i}^{T})^{T} , i=1,2,...,n.$

The statistics $U_{i}={\bf b}^{T}{\bf Y}_{i}-{\bf a}^{T}{\bf
X}_{i}, \hspace{.2in} i=1,2,...,n,$ form a random sample from
$N(\mu,\sigma^{2})$. Denoting
\begin{eqnarray*}
  \bar{U}=\frac{1}{n}\sum\limits_{i=1}^{n}U_{i}, \hspace{.1in}
  S^{2}=\frac{1}{n}\sum\limits_{i=1}^{n}(U_{i}-\bar{U})^{2},\hspace{.1in}
  \bar{\bf
  X}=\frac{1}{n}\sum\limits_{i=1}^{n}{\bf X}_{i},
  \\ \\
  {\bf S_{xx}}=\frac{1}{n}\sum\limits_{i=1}^{n}({\bf X}_{i}-\bar{\bf
  X})({\bf X}_{i}-\bar{\bf X})^{T}, \hspace{.1in} {\bf S_{xy}}=\frac{1}{n}
  \sum\limits_{i=1}^{n}({\bf X}_{i}-\bar{\bf
  X})({\bf Y}_{i}-\bar{\bf Y})^{T},
  \\ \\
  {\bf S_{yx}}=\frac{1}{n}\sum\limits_{i=1}^{n}({\bf Y}_{i}-\bar{\bf
  Y})({\bf X}_{i}-\bar{\bf X})^{T},\hspace{.1in} {\bf S_{yy}}=
  \frac{1}{n}\sum\limits_{i=1}^{n}({\bf Y}_{i}-\bar{\bf
  Y})({\bf Y}_{i}-\bar{\bf Y})^{T},
  \\ \\
   c=\frac{\bar{U}}{\sqrt{n-1}S}=\frac{\sqrt{n}({\bf b}^{T}
   {\bf \bar{Y}}-{\bf a}^{T}{\bf \bar{X}})}{\sqrt{n-1}[{\bf a}^{T}
   {\bf S_{xx}a}-2{\bf a}^{T}{\bf S_{xy}b}+{\bf b}^{T}{\bf
   S_{yy}b}]^{1/2}},
\end{eqnarray*}
the UMVUE $\hat{R}$ of $R={\bf P}\{{\bf a}^{T}{\bf X}>{\bf
b}^{T}{\bf Y}\}$ can be written as (Gupta and Gupta (1990)) $$
 \hat{R}=\left\{ \begin{array}{lr}
     0,& c\geq 1,\\
     \frac{1}{B\left(\frac{1}{2},\frac{n-2}{2}\right)}\int\limits_{c}^{1}
     (1-u^{2})^{\frac{n-4}{2}}du,&
     |c|<1,\\
     1,& c\leq-1.
     \end{array}
     \right.
$$

\vspace{0.5cm}



\leftline{Multivariate discrete probability models}
\vspace{0.5cm}

Consider a class of discrete distributions useful in assessing the
joint survival function of a system with discrete time, sampling
inspection, radioecology, etc. The class is induced by the
following urn scheme. Let an urn contain balls marked by vectors
${\bf a}={(a_1,...,a_m)}^T$, $a_j\in \{\tilde{a}_{j1},
\tilde{a}_{j2},\ldots,\tilde{a}_{jk_j}\}$, $ j=1,...,m$,
$\tilde{a}_{ij}$ -- being arbitrary integers. Denote by $$ {\bf
a}_j=(a_{j1},a_{j2},\ldots,a_{jm})^T, j=1,2,\ldots,K, $$ where
$K=k_1\cdot k_2\cdots k_m, $  all possible values of ${\bf a}$
defined by possible values of their components $a_j$. Let $p_{{\bf
a}_j}$, $j=1,2,\ldots,K$, be the probability for a ball bearing a
vector ${\bf a}_j$ be drawn with replacement and
$\sum\limits_{j=1}^K p_{{\bf a}_j}=1.$

Let a random vector ${\bf X}={(X_1,\ldots, X_m)}^T$ take the value
${\bf r}={(r_1,\ldots, r_m)}^T$ if sums of numbers of $j$th
components of vectors on $n$ balls drawn with replacement are

\[\sum\limits_{i=1}^K a_{ij}l_{{\bf a}_i}=r_j,\]
\[n\min_{1\leq i\leq k_{j}}\{\tilde{a}_{ji}\}
\leq r_j \leq n\max_{1\leq i\leq k_j}\{\tilde{a}_{ji}\},\mbox{   }
j=1,2,...,m,\] where $l_{{\bf a}_i},i=1,\ldots,K$, denotes the
number of balls in the sample bearing the vector ${\bf a}_i$ and
$l_{{\bf a}_i}$th are nonnegative integers such that $
\sum\limits_{i=1}^K l_{{\bf a}_i}=n. $

Evidently that

\begin{equation}
{\bf P}({\bf X}={\bf r},{\bf
p})=\sum\frac{n!}{\prod\limits_{i=1}^K l_{{\bf
a}_i}!}\prod\limits_{i=1}^K p_{{\bf a}_i}^{l_{{\bf a}_i}},
\label{a1213}
\end{equation}
where ${\bf p}=(p_{{\bf a}_1},\ldots,p_{{\bf a}_K})^T$ is the
vector of parameters and the summation is performed over all sets
of nonnegative solutions $l_{{\bf a}_i},i=1,2,\ldots,K$,  of the
system of linear diophantine equations

\begin{equation}
\left\{  \begin{array}{ll} \sum\limits_{i=1}^K a_{ij}l_{{\bf
a}_i}=r_j,  &         j=1,2, ..., m,\\ \sum\limits_{i=1}l_{{\bf
a}_i}=n.
        \end{array}
        \right.
\label{a1214}
\end{equation}
It is easily shown that~(\ref{a1213}) is a proper probability
distribution because $\sum_{\bf r}{\bf P}({\bf X}={\bf r},{\bf
p})=1.$ This distribution can be considered as a generalization of
a probability distribution introduced in Panaretos and Xekalaki
(1986).

To evaluate probabilities~(\ref{a1213}) one needs to find all
solutions of the system of linear equations~(\ref{a1214}) in
nonnegative integers or what is the same to construct all vector
partitions of ${\bf r}$ satisfying the evident relation
$\sum\limits_{i=1}^K l_{{\bf a}_i}=n$. Many algorithms solving a
system of linear Diophantine equations are known. One may use,
e.g., the Smith's technique. Unfortunately, this technique gives
only the partial solution in integers. An algorithm giving the
partial solution in integers in polynomial time is also known, but
the problem of enumerating all solutions in nonnegative integers,
which is $NP$-complete, still exists. Fortunately, using unbiased
estimators for probabilities~(\ref{a1213}) it becomes possible to
evaluate them exactly even for sample sizes of the order of $n=16$
(Nikulin et al (1999).

 Since $ \sum\limits_{i=1}^K p_{{\bf a}_i}=1, $
probabilities~(\ref{a1213}) depend on $K-1$ parameters. We may
consider ${\bf P}({\bf X}={\bf r},{\bf p})$ as the sum of
multinomial probabilities

\begin{equation}
{\bf P}({\bf Y}={\bf l},{\bf
p})=\frac{n!}{\prod\limits_{i=1}^Kl_{{\bf a}_i}!}
\prod\limits_{i=1}^K p_{{\bf a}_i}^{l_{{\bf a}_i}}, \label{a1215}
\end{equation}
where ${\bf l}=(l_{{\bf a}_1},\ldots,l_{{\bf a}_k})^T$. It is well
known that the complete sufficient for parameters of
~(\ref{a1215}) statistic is the subvector ${\bf Z}_{K-1}={(Z_1,
\ldots, Z_{K-1})}^T$ of the vector ${\bf Z}={({\bf Z}_{K-1},
Z_{K})}^T, $ where $Z_{K}=nM-\sum\limits_{i=1}^{K-1}Z_i,$ $M$
being the number of samples, each of which contains $n$ balls, and
$ Z_j=\sum\limits_{i=1}^{M}X_{ij},\quad j=1,2,\ldots, K-1,$ where
$X_{ij}$ is the number of balls bearing ${\bf a}_j$ in the $i$th
sample of $n$ balls.

Using a standard technique (see, e.g., Voinov and Nikulin (1996))
it is easily proved that the UMVUE ${\cal {P}}_0={\bf P}({\bf
X}={\bf r}|{\bf Z})$  of the probability ${\bf P}({\bf X}={\bf
r},{\bf p})$  is

\begin{equation}
{\cal {P}}_0=\sum A({\bf Z},{\bf l}), \label{a1216}
\end{equation}
where
\begin{equation}
A({\bf Z},{\bf l}) = \left\{  \begin{array}{lll}
\frac{\prod\limits_{i=1}^{K}{Z_i\choose {l_{{\bf a}_i}}}}
{{nM\choose{n}}}, & Z_i\geq l_{{\bf a}_i}, & M\geq 1,  \\ 0, &
\mbox { otherwise,}
   \end{array}
        \right.
\label{a1217}
\end{equation}
and the summation in~(\ref{a1216}) is performed over all sets of
nonnegative solutions $l_{{\bf a}_i}$ of~(\ref{a1214}).

 The UMVUE ${\cal P}_1$ of ${\bf Var}{\cal P}_0$  is

\begin{equation}
{\cal P}_1 = {\cal P}_0^2 - {\cal P}_2, \label{a1218}
\end{equation}
where
\begin{equation}
{\cal P}_2 = \sum B({\bf Z}, {\bf l})+2\sum\limits_{i>j} C({\bf
Z}, {\bf l}), \label{a1219}
\end{equation}

\begin{equation}
B({\bf Z},{\bf l}) = \left\{  \begin{array}{lll} \frac{
\prod\limits_{k=1}^{K} {2l_{{\bf a}_k}\choose{l_{{\bf a}_k}}}
{Z_k\choose{2l_{{\bf a}_k}}} } { {2n\choose{n}}{nM\choose{2n}} },
& Z_k\geq 2l_{{\bf a}_k}, & M\geq 2, \\ 0,  &  \mbox{ otherwise}.
   \end{array}
        \right.
\label{a1220}
\end{equation}
and
\begin{equation}
C({\bf Z},{\bf l}) = \left\{  \begin{array}{lll} \frac{
\prod\limits_{k=1}^{K} {l_{{\bf a}_k}^i+l_{{\bf
a}_k}^j\choose{l_{{\bf a}_k}^i}} {Z_k\choose{l_{{\bf
a}_k}^i+l_{{\bf a}_k}^j}} } { {2n\choose{n}}{nM\choose{2n}} }, &
Z_k\geq l_{{\bf a}_k}^i+l_{{\bf a}_k}^j , & M\geq 2, \\ 0,  &
\mbox{ otherwise,}
   \end{array}
        \right.
\label{a1221}
\end{equation}
$\{l_{{\bf a}_1}^{i}, \ldots, l_{{\bf a}_K}^{i} \}$,
$i=1,2,\ldots,q$,  is the $i$th set of arbitrarily ordered sets of
solutions of the system of equations~(\ref{a1214}), $1\leq i,j
\leq q,$ where $q$ is the total number of sets of solutions for
${\bf r}$ given. The first sum in~(\ref{a1219}) is performed over
all $q$ sets and the second one over sets with $i>j.$

 Formulas~(\ref{a1216}) and~(\ref{a1219}) are easily proved by direct
evaluation of their expectations with respect to the probability
distribution of the complete sufficient for parameters $p_{{\bf
a}_i}$ statistic ${\bf Z}_{K-1}$.  It is worthy to note that
UMVUEs ${\cal {P}}_0$ are probabilities themselves, i.e.,
$\sum_{\bf r}{\cal {P}}_0=1$. Another words, they are range
preserving in the sense of Hoeffding (1983).

\vspace{0.5cm}


\leftline{\bf References}
\vspace{0.5cm}

\begin{enumerate}
\item Basu A.P. (1988). Multivariate exponential distributions and
their applications in reliability, {\it P.R.Krishnaiah and C.R.Rao
(eds.), Handbook of Statistics}, {\bf 7}, 467--476.

\item Belyaev Yu.K. (1975).{\it Probabilistic methods of sampling
inspection},\\ Nauka, Moscow (in Russian).

\item Block H.W., Basu A.P. (1974). A continuous bivariate
exponential extension, {\it JASA}, {\bf 69}, 1031--1037.

\item Chandra S. and Owen D.B. (1975). On estimating the
reliability of a component subject to several different stresses
(strengths), {\it Naval Res. Logist. Quart.}, {\bf 22}, 31--39.

\item Downton F. (1973). The estimation of ${\bf P}\{Y<X\}$ in the
normal case, {\it Technometrics}, {\bf 15}, 551--558.

\item Gupta R.D., Gupta R.C. (1990). Estimation of $Pr({\bf
a}'{\bf x}>{\bf b}'{\bf y})$ in the multivariate normal case, {\it
Statistics}, {\bf 21}, 91--97.

\item Hoeffding W. (1983). Range-preserving unbiased estimators in
the multinomial case, {\it Inst.Statist.Mimeo Ser.}, {\bf 1515},
1--10.

\item Klein J.P., Basu A.P. (1985). Estimating reliability for
bivariate exponential distributions, {\it Sankhy\'{a}}, {\bf 47B},
346--353.

\item Lieberman G.J., Resnikoff G.J. (1955). Sampling plans for
inspection by variables, {\it JASA}, {\bf 50}, 457--516.

\item Lindley D.V., Singpurwalla N.D. (1986). Multivariate
distributions for the life lengths of components of a system
sharing a common environment, {\it J. Appl. Prob.}, {\bf 23},
418--431.

\item Mazumdar M. (1983). Some estimates of reliability using
interference theory, {\it Naval Res. Logist. Quart.}, {\bf 17},
159--165.

\item Neyman J., Scott E.L. (1960). Correction for bias introduced
by a transformation of variables, {\it Ann. Math. Stat.}, {\bf
31}, 645--656.

\item Nikulin M.S., Novak M.M., Smirnov T.I., Voinov V.G. (1999).
A probabilistic description of radioactive contamination, (to be
published in {\it Applied Radiation and Isotops}).

\item Nikulin M.S., Price P.M., Turetayev D., Voinov V.G. (1999).
An application of a multivariate discrete probability distribution
for sampling survey analysis of supply chains, {\it Technical
Report N9905, University Victor Segalen Bordeaux 2, Bordeaux}.

\item Panaretos J., Xekalaki E. (1986). On generalized binomial
and multinomial distributions and their relation to generalized
Poisson distributions, {\it Ann.Inst.Statist.Math.}, {\bf 38A},
223--231.

\item Pugh E.L. (1963). The best estimate of reliability in the
exponential case, {\it Oper. Res.}, {\bf 11}, 57--61.

\item Rohatgi V.K. (1989). Unbiased estimation of parametric
functions in sampling from one-truncation parameter families, {\it
Austr. J. Stat.}.

\item Selvavel K. (1989). Unbiased estimation in sampling from
two-trun\-cation parameter families when both samples are type II
censored, {\it Commun.Stat.:Theory and Meth.}, {\bf 18},
3519--3531.

\item Unnikrishnan Nair N., Asha G. (1997). Some classes of
multivariate life distributions in discrete time, {\it J. of
Multivar. Anal.}, {\bf 62}, 181--189.

\item Voinov V.G. (1984). On unbiased estimation of ${\bf
P}\{Y<X\}$ in the normal case, {\it Zap. Nauchn. Sem. LOMI}, {\bf
136}, 5--12(in russian).

\item Voinov V.G., Nikulin M.S. (1993). {\it Unbiased estimators
and their applications. Vol.1:Univariate case},
Dordrecht/Boston/London: Kluwer.

\item Voinov V.G., Nikulin M.S. (1996). {\it Unbiased estimators
and their applications. Vol.2:Multivariate case},
Dordrecht/Boston/London: Klu\-wer.

\end{enumerate}

\end{document}