\documentstyle[12pt,draft]{article}            
\oddsidemargin 0cm
\textwidth 14.8cm
\textheight 24.2cm
\topmargin -1.5cm                                                                                                                                               
               


\def\eqnum#1{\eqno(#1)}
\def \suml{\sum\limits}
%\def \un{\epsfbox{un.eps}}
%\def \Un{\mathop{{\lower3pt\hbox{\epsfxsize=10pt\epsfbox{un.eps}}}}}
%\def \petitun{\mathop{{\lower3pt\hbox{\epsfxsize=4pt\epsfbox{un.eps}}}}}
\def\PP{\mbox{\boldmath$P$}}




%17.02.99


\begin{document}
\medskip
\centerline{{\Large \bf On parametric inference for step-stress models}}  

\vspace{0.5cm}

\centerline{\bf Vilijandas B. Bagdonavi\v{c}ius}
\centerline{University of Vilnius, Lithuania }
\vspace{0.5cm}

\centerline{\bf Mikhail S. Nikulin}
\centerline{Steklov Mathematical Institute, St.Petersburg, Russia $\&$}
\centerline{ University Victor Segalen Bordeaux-2, France}
\vspace{1cm}


\noindent {\small {\bf Key Words} -  Accelerated life testing; Additive accumulation of damages; 
Cox model; Step-stress. }
\vspace{1cm}


\noindent {\small {\bf Summary \& Conclusions} -  Applications of the additive accumulation of  damages (AAD) and 
the proportional hazards (PH) models in accelerated life testing with step-stresses are discussed. A new 
model including both of them is proposed. It is more natural then PH model and wider then the AAD model. 
Maximum likelihood function construction is discussed.}
\vspace{0.5cm}


\centerline{1. INTRODUCTION}
\vspace{0.5cm}

\noindent {\it Acronyms}\\
\begin{tabular}{rl} 
MLE & maximum likelihood estimator(s)\\ 
ALT & accelerated life testing \\ 
PH  & proportional hazards\\
AAD & additive accumulative of damages\\
cdf & cumulative distribution function(s)\\
pdf & probability density function(s)\\
ahr & accumulated hazard rate\\
GPH & generalized PH\\ 
GLPH& generalized linear PH\\
TFR & tampered failure rate
\end{tabular}
\vspace{0.5cm}

\noindent {\it Notations}\\

\begin{tabular}{rl} 
$x_0$                & design stress\\
$k$                  & number of test stresses\\
$j$                  & index for test-stresses : $j=1,...,k$ unless otherwise specified\\ 
$x_j$                & ordered test-stresses\\
$x(\cdot)$           & time varying stress, particularly a step-stress\\
$n$                  & number of units placed on test\\
$T_{(i)}$            & ordered failure times\\
$t_i$                & stress-change times, $i=1,...,k-1$\\ 
$t_f$                & censoring time at $x_k$\\
$^{\prime}$          & first derivative\\
\end{tabular}

\begin{tabular}{rl} 
$r,\, \varphi$       & functions of a stress\\
$r_j$                & $r(x_j)$\\
$\Delta t_j$         & $t_j-t_{j-1}$\\
$s_i$                & $\sum^i_{j=1}r_j\Delta t_j$\\
$\rho$               & $r(x_1)/r(x_0)$\\
$F_{x(\cdot)}$       & Cdf of times-to-failure under the stress $x(\cdot)$\\
$F_0$                & baseline Cdf\\
$f_{x(\cdot)}$       & pdf of times-to-failure under the stress $x(\cdot)$\\
$\lambda_{x(\cdot)}$ & hazard rate under the stress $x(\cdot)$\\
$\Lambda_{x(\cdot)}$ & $\int_0^t \lambda_{x(\cdot)}(u)du$ -ahr under  $x(\cdot)$\\
$\lambda_0$          & baseline hazard rate\\
$\Lambda_0$          & $\int_0^t \lambda_0(u)du$ - baseline ahr\\
$sp_{x(\cdot)}$       & the right end  of the support of distribution under  $x(\cdot)$\\
$\beta,\, \delta,\,\gamma,\,\alpha_0$ & unknown parameters\\
$k(t)$               & $min\{j:t_j\geq t\}$ \\
$\theta_i$             & $r_0\alpha_0/r_i$\\
$\hat{\beta},\, \hat{\delta},\,\hat{\gamma},\,\hat{\alpha}_0$ & ML estimators of $\beta,\, \delta,\,\gamma,\,\alpha_0$ \\
$T_{x(\cdot)}$       & the time-to-failure under the stress $x(\cdot)$\\
$R$                  & $\Lambda_{x(\cdot)}(T_{x(\cdot)})$-the exponential resource
\end{tabular}
\vspace{1cm}


\par In ALT units are tested at higher-then-usual levels of stress
to induce early failures. The results are extrapolated to estimate the life distribution at 
the design stresss using models which relate the lifetime to the stress.
\par In step-stress testing a unit is placed on test at an initial low stress and if it does
not fail in a predetermined time $t_1$, the stress is increased and so on.
\par Denote by $x(\cdot)$ a step-stress of the form
$$
x(\tau) = \left\{
\begin{array}{ll}
x_1,&t_0\leq \tau < t_1,\\
x_2,&t_1\leq \tau < t_2,\\ 
\cdots& \cdots\\
x_k,&t_{k-1}\leq \tau < t_k,
\end{array}\right. \eqnum{1}
$$ 
where $t_0 = 0$, $t_k = \infty$.
\par The two most known {\it general} models of ALT with time varying stresses are: the AAD model 
(Bagdonavi\v{c}ius [1]) and the PH model (Cox [2]).
\par The {\it AAD model}: 
$$
F_{x(\cdot)}(t)= F_0(\int^t_0 r[x(\tau)]d\tau). \eqnum{2}
$$
\par The {\it PH model}:
$$
\lambda_{x(\cdot)}(t)= r[x(t)]\lambda_0(t). \eqnum{3}
$$
\par The AAD and PH models are rather restrictive. In the case of AAD model  the stress changes locally 
only the  scale. If the PH model is true, $x_0$ is a usual stress, $x_1$ is an accelerated with 
respect to $x_0$ stress, $x_0,x_1 \in {\cal E}$, i.e. $S_{x_0}(t)\geq S_{x_1}(t)$ for all $t \geq 0$,   
  and 
$$
x(t)=\left\{
\begin{array}{ll}
x_1,& 0 \leq t < t_1,\\
x_0,& t \geq t_1,
\end{array} \right.
$$
is a step-stress, then $\alpha_{x(\cdot)}(t)=\alpha_{x_0}(t)$ for all $ t_1 > 0$, $t > t_1$. So if one 
group of items is tested under the usual stress and the 
second group - under an accelerated stress $x_1$ untill the moment $t_1$ and after this moment both groups are observed under the same 
usual stress $x_0$, the failure rate after the moment $t_1$ is the same for both groups. If items 
are aging, it is not  natural.
\vspace{0.5cm}

\centerline{2. FORMULATION OF AAD AND PH MODELS FOR CONSTANT}
\centerline{ AND STEP-STRESSES}

\vspace{0.5cm}


In terms of Cdf the PH model (3) is written:
$$
F_{x(\cdot)}(t)= 1 - exp\{-\int^t_0 r[x(\tau)]d\Lambda_0(\tau)\}. \eqnum{4}
$$

\par In the particular case of constant in time stress $x_j$ and the step-stress (1) the AAD
model (2) implies 
$$
F_{x_j}(t) = F_0(r_jt),
$$
and for $t\in [t_{i-1}, t_i)$, $(i=1,2,...,k)$
$$
F_{x(\cdot)}(t) =F_0\left(\sum^{i-1}_{j=1}r_j\Delta t_j+r_i(t-t_{i-1})\right).
$$

In terms of $F_{x_j}$,
$$
F_{x(\cdot)}(t)=F_{x_i}(s_{i-1}+t-t_{i-1}), \eqnum{5}
$$
where $s_i$ can be determined recurrently from the equations
$$
F_{x_{i+1}}(s_i)=F_{x_i}(s_{i-1}+t_i-t_{i-1}). \eqnum{6}
$$
For the step-stress (1) and $k=2$ the PH model (3) implies
$$
\lambda_{x(\cdot)}(\tau)=\left\{\begin{array}{ll}
\lambda_{x_0}(\tau), & 0 \leq \tau \leq t, \\
\rho \lambda_{x_0}(\tau), & \tau > t. 
\end{array}
\right.
$$
\par In this particular case of step-stresses the Cox model is sometimes called the TFR 
model (Bhattacharyya and Stoejoeti (1989)). If $k\geq 2$, the formula (4) implies 
$F_{x_j}(t)=1-e^{-r_j\Lambda_0(t)}$ and for $t\in [t_{i-1},t_i)$
$$
F_{x(\cdot)}(t)=1-\exp\{-\sum^{i-1}_{j=1}r_j(\Lambda_0(t_j)-\Lambda_0(t_{j-1}))-
r_i(\Lambda_0(t)-\Lambda_0(t_{i-1}))\}. \eqnum{7}
$$
In the case when the distribution under the design stress $x_0$ is Weibull :
$$
F_{x_0}(t)=1-e^{-t^\delta/ \alpha_0}, \quad \Lambda_0(t)=\frac{t^\delta}{r_0\alpha_0},
$$
we have
$$
F_{x_i}(t)= 1-e^{-t^\delta/ \theta_i},
$$
 and for $t\in [t_{i-1},t_i)$
$$
F_{x(\cdot)}(t)= 1-\exp\left\{-\sum^{i-1}_{j=1}\frac{t_j^\delta-t_{j-1}^\delta}{\theta_j}-
\frac{t^\delta-t_{i-1}^\delta}{\theta_i}\right\}.
$$
This particular case of the PH model was considered by Khamis and Higgins [4].
\vspace{0.5cm}

\centerline{ 3. NEW MODEL}
\vspace{0.5cm}

Consider generalisation of both the AAD and PH models. For this purpose we formulate 
these models in other terms.  For any $x(\cdot)$ the random variable $R=\Lambda_{x(\cdot)}(T_{x(\cdot)})$ has the standard 
exponential distribution with the survival function $S_R(t)=e^{-t}$, $t \geq 0$ and is called the 
exponential resource, the number $\Lambda_{x(\cdot)}(t)$ being the  exponentional resource used untill 
the moment $t$ under the stress $x(\cdot)$. So the failure rate $\lambda_{x(\cdot)}(t)$ is the rate 
of exponential resource usage.
\par Under the PH model we have 
$$
\Lambda^{\prime}_{x(\cdot)}(t)=r[x(t)]\lambda_0(t),
$$
i.e. the rate of resource usage is proportional to some baseline rate and the constant of 
proportionality is a function of a stress.  
\par GPH model is true on ${\cal E}$ if for all $x(\cdot) \in {\cal E}$ 
$$
\lambda_{x(\cdot)}(t)=r\{x(t)\}q\{\Lambda_{x(\cdot)}(t)\}\lambda_{0}(t). 
$$
The rate of resource usage at the moment $t$ depends not only on a stress applied at this moment but also on the resource used untill $t$. 
\par The particular cases of the GPH model are the PH model ($q(u)\equiv 1$) and the AAD model 
($\lambda_0(t) \equiv \lambda_0=const$). 

\par The function $r$ in the GPH model (as in the PH model)  can be parametrized in the following way :
$$
r[x(\tau)]=e^{\beta^T \varphi(x(\tau))},
$$
where $\varphi$ is some known function of the stress. For example, if $\varphi(x)=x, \, \ln{x}, \, 
1/x$, we have generalizations of the 
 log-linear,   power rule,  Arrhenius models, respectively, see, for example, 
Nelson (1990). We'll write $x(\cdot)=\varphi(x(\cdot))$ in what follows.
\par Consider the case when the rate of exponential
 resource usage with respect to the baseline rate is proportional to a monotone function of the used 
resource. One of the most natural  parametrizations of the function $q$ in this case would be 
$$
q(u)=e^{\gamma u}, \quad \gamma \in {\bf R}.
$$
So we consider
\par GLPH model :
$$
\lambda_{x(\cdot)}(t)=e^{\beta x(t)+\gamma \Lambda_{x(\cdot)}(t)}\lambda_{0}(t). \eqnum{8}
$$
This model implies
$$
F_{x(\cdot)}(t)=\left\{\begin{array}{ll}
1-\{1-\gamma\int^t_0e^{\beta x(u)}d\Lambda_0(u)\}^{1/\gamma},&\gamma \neq 0\\
\exp\{-\int^t_0e^{\beta x(u)}d\Lambda_0(u)\},&\gamma = 0
\end{array}\right.
$$
If $\gamma \leq 0$, the support of distribution is $[0,\infty)$, if $\gamma > 0$, it is  
$[0,sp_{x(\cdot)})$, where $sp_{x(\cdot)}$ verifies the equation 
$$
\int^{sp_{x(\cdot)}}_0e^{\beta x(u)}d\Lambda_0(u)=1.
$$ 
In the particular cases of step-stress (1) we have the formula (7) when $\gamma =0$ ; 
if $\gamma \neq 0$ then  for all $t \in [t_{i-1},t_i)$
$$
F_{x(\cdot)}(t)=1-\{1-\gamma\sum^{i-1}_{j=1}e^{\beta x_j}(\Lambda_0(t_j)-\Lambda_0(t_{j-1}))-
\gamma e^{\beta x_i}(\Lambda_0(t)-\Lambda_0(t_{i-1}))\}^{1/\gamma} \eqnum{9}
$$
and for the constant stress $x_j$
$$
F_{x_j}(t)=1-\{1-\gamma e^{\beta x_j}\Lambda_0(t)\}^{1/\gamma}.
$$
In the case when the distribution of time-to-failure under the design stress $x_0$ is Weibull, i.e. 
$$
F_{x_0}(t)=1-e^{-t^\delta/\alpha_0}, \quad \lambda_0(t)=\frac{1}{\gamma}e^{-\beta x_0}\left(1-
e^{-\frac{\gamma}{\alpha_0}+\delta}\right),
$$
then under the GLPH model (8) and the step-stress (1)
$$
F_{x(\cdot)}(t)=1-\{1+\sum^{i-1}_{j=1}e^{\beta (x_j-x_0)}[e^{-\gamma t_j^\delta/\alpha_0}-
e^{-\gamma t_{j-1}^\delta/\alpha_0}]+
$$
$$
e^{\beta (x_i-x_0)}[e^{-\gamma t^\delta/\alpha_0}-
e^{-\gamma t_{i-1}^\delta/\alpha_0}]\}^{1/\gamma}. \eqnum{10}
$$
Note that if the rate of resource usage depends on used resource ( i.e. $\gamma \neq 0$), the distribution 
of the time-to-failure under constant in time stresses $x_j \neq x_0$ is not Weibull :
$$
F_{x_j}(t)=1-\{1-e^{\beta (x_j-x_0)}[1-e^{-\gamma t^\delta/\alpha_0}]\}^{1/\gamma}.
$$
Under the GLPH model distributions of times-to-failure under constant in time stresses are from the same family of 
distributions if we take one of the following two families of distributions : \\
\noindent The family $K_1$ :
$$
F_{x_0}(t)=1-(1+\alpha_0 t^\delta)^{1/\gamma}, \quad t \geq 0\quad (\alpha_0 > 0, \, \delta > 0,\,
\gamma < 0).
$$
If $0<\delta \leq  1$,  the  failure rate $\lambda_{x_0}(t)$ is {\it decreasing}. If $\delta > 1$, then the 
failure rate is increasing until the moment $t=\left(\frac{\delta -1}{\alpha_0}\right)^{1/\delta}$ and 
 decreasing after this moment i.e. has the $\bigcap $-{\it shape}.\\
\noindent The family $K_2$ :
$$
F_{x_0}(t)=1-(1+\alpha_0 t^\delta)^{1/\gamma}, \quad t \in [0,\left(-\frac{1}{\alpha_0}\right)^{1/\delta})
\quad (\alpha_0 < 0, \, \delta > 0,\, \gamma > 0).
$$
If $0<\delta < 1$, the hazard  rate is decreasing in the interval $[0, \left(\frac{\delta -1}{\alpha_0}\right)^{1/\delta})$ and increasing 
in the interval  $(\left(\frac{\delta -1}{\alpha_0}\right)^{1/\delta}, \left(-\frac{1}{\alpha_0}\right)^{1/\delta})$. 
The hazard  rate has $\bigcup$-{\it shape} . So this family is very appealing. If $\delta > 1$, the hazard rate is 
{\it increasing} on $(0,\left(-\frac{1}{\alpha_0}\right)^{1/\delta})$.
\par Suppose that the time-to-failure distribution is from one of the above mentioned families. 
If $ F_{x_0}(t) \in K_i$ (i=1,2) and the GLPH model is satisfied then
$$
\Lambda_0(t)=-\frac{\alpha_0}{\gamma}e^{-\beta x_0}t^{\delta}
$$
and 
$$
F_{x_j}(t)=1-(1+\alpha_jt^{\delta})^{1/\gamma},
$$
where $\alpha_j=\alpha_0e^{\beta (x_j-x_0)}$. If $F_{x_0} \in K_2$, the support of the distribution is 
finite : $[0,\left(\frac{1}{-\alpha_j}\right)^{\frac{1}{\delta}})$. 
 If $x(\cdot)$ is a step stress of the form (1) then
$$
F_{x(\cdot)}(t)=1-\{1+\alpha_0 \varphi(t,\beta,\delta)\}^{1/\gamma}, \quad t < sp_{x(\cdot)}, \eqnum{11}
$$
where 
$$
\varphi(t,\beta,\delta)=\sum^{k(t)-1}_{j=1}e^{\beta(x_j-x_0)}(t_j^\delta -t_{j-1}^\delta)+e^{\beta(x_{k(t)}-x_0)}(t^\delta - 
t_{k(t)-1}^\delta).
$$
If $\alpha_0 > 0$, $\gamma < 0$, then $sp_{x(\cdot)} = +\infty$. 
If $\alpha_0 < 0$, $\gamma > 0$, the right limit of the support $(0, sp_{x(\cdot)})$ 
 verifies the equation $\varphi(sp_{x(\cdot)},\beta,\delta)=-1/\alpha_0$, i.e.
$$
sp_{x(\cdot)}=\{t_{k(t)-1}^\delta-\frac{1}{\alpha_0}e^{\beta(x_{k(t)}-x_0)}-\sum^{k(t)-1}_{j=1}e^{\beta(x_j-x_0)}(t_j^\delta -t_{j-1}^\delta).
$$
The probability density is 
$$
f_{x(\cdot)}(t)=\frac{\delta}{\gamma}e^{\beta(x_{k(t)}-x_0)}t^{\delta -1}\{1+\alpha_0
\varphi(t,\beta,\delta)\}^{1/\gamma -1}, \quad t < sp_{x(\cdot)},
$$
and the failure rate 
$$
\lambda_{x(\cdot)}(t)=\frac{\delta}{\gamma}e^{\beta(x_{k(t)}-x_0)}\frac{t^{\delta-1}}
{(1+\alpha_0 \varphi(t,\beta,\delta))}
$$
has the same shape as in the case of constant in time stresses (decreasing or 
$\bigcap$-shape for the family $K_1$, increasing or $\bigcup$-shape for the family $K_2.$
\vspace{0.5cm}

\centerline{4. ESTIMATION}
\vspace{0.5cm}


\par Suppose that all units are initially placed on test at $x_1$ and run until $t_1$ when the stress 
is changed to $x_2>x_1$. At $x_2$, testing continues until $t_2$ when stress is changed to $x_3>x_2$, 
and so on, until $x_k$. Under $x_k$, testing continues until all remaining units fail or until $t_f$ 
whichever comes first.
\par Denote by $D$ the  number of failures during the experiment, $T_{(1)} \leq ... \leq T_{(D)}$ the 
observed moments of failures. In the case of a finite support the following reparametrization can be 
introduced: $\sigma=sp_{x(\cdot)},\,\beta,\,\gamma,\,\delta$. 
Then $\alpha_0=-1/\varphi(\sigma,\beta,\delta)$. 
The likelihood function: \\
\noindent  for the family $K_1$ 
$$
L(\alpha_0, \beta,\gamma,\delta)=\left(\frac{\delta}{\gamma}\right)^De^{\beta\sum^D_{i=1}
(x_{k(T_{(i)})}-x_0)}\times 
$$
$$
\prod^D_{i=1}T_{(i)}^{\delta -1}\{1+\alpha_0 \varphi(T_{(i)};
\beta,\delta)\}^{\frac{1}{\gamma}-1}\{1+\alpha_0 \varphi(t_f;
\beta,\delta)\}^{\frac{n-D}{\gamma}}
$$
and for the family $K_2$ 
$$
L(\sigma, \beta,\gamma,\delta)=\left(\frac{\delta}{\gamma}\right)^De^{\beta\sum^D_{i=1}
(x_{k(T_{(i)})}-x_0)}\times
$$
$$
\prod^n_{i=1}T_{(i)}^{\delta -1}\left(1-\frac{ \varphi(T_{(i)},
\beta,\delta))}{\varphi(\sigma,\beta,\delta)}\right)^{\frac{1}{\gamma}-1}
\left(1+\frac{ \varphi(t_f,\beta,\delta))}{\varphi(\sigma,\beta,\delta)}\right)^{\frac{n-D}{\gamma}} 
{\bf 1}\{T_{(D)} < \sigma\}.
$$
\par The Cdf $F_{x_0}(t)$ is estimated by 
$$
 \hat{F}_{x_0}(t)=1-(1+\hat{\alpha}_0t^{\hat{\delta}})^{\frac{1}{\hat{\gamma}}}.
$$
In the case of the class $K_1$ approximate confidence intervals for $F_{x_0}(t)$ and other reliability
 characteristics under the design stress $x_0$ are obtained by standard ML methods using the estimated Fisher 
information matrix. In the case of the class $K_2$ when the support is finite, rigorous asymptotic theory
 is to be established.
If the time-to-failure under the design stress $x_0$ has the Weibull (or another) distribution, then estimation under 
the GLPH model can be done similarly taking the formula (10) (or (9) with respective $\Lambda_0$) for 
$F_{x(\cdot)}$ instead of (11) in all formulas.
\vspace{1cm}



\centerline{ \sc References}
\def\ref{\smallskip\hang\noindent}
\vspace{0.5cm}

\ref [1] Bagdonavi\v{c}ius, V., (1978).  Testing the hyphothesis of the additive accumulation of damages. {\it Probab. Theory and its Appl.},  
{\bf 23} (2) , 403-408.
  
\ref [2] Bhattacharyya, G.K., Stoejoeti, Z., (1989). A tampered failure rate model for step-stress accelerated life model, 
{\it Communication in Statistics. - Theory Meth.},{\bf 18} (5), 1627-1643.

\ref [3] Cox, D.R., (1972). Regression models and life tables, J.R.Statist. Soc., 
B, {\bf 34}, 187-220.

\ref [4] Khamis, I.H. and Higgins, J.J., (1998). A new model for step-stress testing. {\it IEEE Transactions on Reliability}, 
{\bf 47}, 2, 131-134.

\vspace{1cm}

\centerline{AUTHORS}
\vspace{0.5cm}

\noindent  Prof. Vilijandas Bagdonavi\v{c}ius ; Dept. of Mathematics Statistics; Vilnius University; 
Vilnius 2049, Lithuania.\\
\noindent  e-mail :      Vilijandas.Bagdonavicius@maf.vu.lt\\
\noindent  e-mail :      Vilijandas.Bagdonavicius@mi2s.u-bordeaux2.fr
\par {\bf Vilijandas Bagdonavi\v{c}ius} is a Professor in the Department of Mathematical Statistics of 
Vilnius University; Invited Professor of University Victor Segalen Bordeaux 2 in 1994-1999 (half time). 
He received his Ph.D in Mathematics and Physics in 1979 from Vilnius University. Author of many articles on mathematical statistics, 
survival analysis and reliability; co-author with M.Nikulin of "Semiparametric models in accelerated life testing", v.98, (1995) and 
"Additive and Multiplicative Semiparametric Models in Accelerated life testing and Survival Analysis", v.108, (1998), published in 
"Queen's Papers in Pure and Applied Mathematics", Kingston, Ontario, Canada.
\vspace{0.5cm}

\noindent Prof. Mikhail Nikulin ; UFR MI2S ; University Victor Segalen Bordeaux 2; 146 Rue de Leo Saignat 33076 Bordeaux, France
 \& The Laboratory of Statistical Methods of the Steklov Mathematical Institute at St. Petersburg, Russia.\\
\noindent  e-mail : M.S.Nikouline@mi2s.u-bordeaux2.fr
\par {\bf Mikhail Nikulin} earned his Ph.D in the Theory of Probability and Mathematical Statistics (1973) from The Steklov Mathematical 
Institute in Moscow. Member of the St.Petersburg's Mathematical Society, the elected member of the International Statistics Institute, 
a board member of the newly book series "Statistical Distributions and Models with Applications", Gordon and 
Breach Science Publishers. The co-author with V.Voinov of " Unbiased Estimations and Their Applications", V.1 : Univariate Case, (1993); 
V.2 : Multivariate Case (1996), Kluwer Academic Publishers, and the co-author with P.E.Greenwood of "A Guide to Chi-squared Testing", John 
Wiley and Sons.
 
 \end{document}