1. Introduction
In many applications where two random variables represent the lifetime of different systems, it is of interest to study their aging properties. This will allow to determine which system is performing better with respect to some given property: the aging rate, lifetime expectancy, skewness of lifetimes, etc. For this purpose, stochastic ordering between random variables provides a convenient way to describe such comparisons. These orderings may be defined through relations between distributions, survival or failure rate functions of the relevant random variables. The monographs by Shaked and Shantikumar [Reference Shaked and Shantikumar18] or Marshall and Olkin [Reference Marshall and Olkin16], give a good account of various stochastic orders and their applications.
We will be interested in the star-shaped order, introduced by Barlow and Proschan [Reference Barlow and Proschan4] and defined by a monotonicity property on a suitable transformation on the distribution functions, as expressed by Definition 1 below. It can be easily seen that the definition is equivalent to allowing at most one crossing between the distribution functions of scaled lifetimes, as referred in Proposition C.11 of Marshall and Olkin [Reference Marshall and Olkin16] or 4.B.2 in Shaked and Shantikumar [Reference Shaked and Shantikumar18]. It follows from this characterization that the star-shaped order may, thus, be interpreted as a comparison of the lifetime aging rate for systems that started functioning simultaneously. From the practical point of view, since the distribution functions of the lifetime variables under comparison often do not have an explicit formula, it may be technically difficult to verify the star-shaped ordering. Thus, it becomes relevant to establish equivalent conditions for which the monotonicity of the referred function holds. Using the sign technique referred to in Marshall and Olkin [Reference Marshall and Olkin16] or Shaked and Shantikumar [Reference Shaked and Shantikumar18], Arab and Oliveira [Reference Arab and Oliveira1] analyzed the ordering relationships within the Gamma and the Weibull families of distributions, later extended to the comparison of lifetimes of parallel systems in Arab et al. [Reference Arab, Oliveira and Hadjikyriakou2]. However, when the underlying distributions depend on a large number of parameters, this sign analysis becomes rather hard to control and often does not allow for a conclusion. An alternative approach may be based on a criterion proposed by Saunders and Moran [Reference Saunders and Moran17] when the distribution functions depend on a single real parameter. This criterion turned out to be useful to exhibit order relations within parametric families of distributions (see [Reference Khaledi and Kochar9] or [Reference Kochar and Xu11,Reference Kochar and Xu12] among many others). As what concerns the lifetimes of more complex systems, the Saunders and Moran's criterion was used by Kochar and Xu [Reference Kochar and Xu10] to obtain a characterization for parallel systems each one formed by two types of components with exponentially distributed lifetimes. More recently, Arab et al. [Reference Arab, Oliveira and Hadjikyriakou2] proved that the lifetimes of parallel systems with homogeneous and independent exponential components get smaller (or age faster) with respect to the star-shaped order, as the number of components increases. We note that, for the case of series and parallel systems with more than two heterogeneous and, especially, nonexponentially distributed components, not much work seems to have been done regarding the star-shaped comparability. Since these are complex models, depending on more than one parameter, the Saunders and Moran's [Reference Saunders and Moran17] result cannot be used. Comparability with respect to star-shaped ordering for multiparameter models has recently been addressed by Belzunce et al. [Reference Belzunce, Pinar, Ruiz and Sordo5] and Arriaza et al. [Reference Arriaza, Belzunce and Martinez-Riquelme3], who proved sufficient conditions allowing to order within a few parametric families of distributions, including the generalized Gamma distribution and types I and II Beta generalized distribution. We propose a different approach, exploring extensions of the Saunders and Moran's criterion to order models depending on several parameters.
This paper is structured as follows. In Section 2, we present the extension of the Saunders and Moran's [Reference Saunders and Moran17] criterion to families of distributions depending on multidimensional parameters. In Section 3, we discuss a few applications of the obtained criterion to complex systems with heterogeneous components, describing conditions on the parameters so that the star-shaped comparability holds when the lifetimes of the components satisfy suitable proportionality assumptions, including the popular proportional hazard rate (PHR) and proportional reversed hazard rate (PRHR) models.
2. A criterion for the star-shaped order
Let 
$\mathcal {F}$ denote the family of distribution functions vanishing at 
$0$ with support contained in 
$[0, +\infty )$. Let 
$X$ be a nonnegative random variable with distribution function 
$F_X\in \mathcal {F}$, density function 
$f_X$, and survival function 
$\bar {F}_X$. We start by defining the star-shaped order relation, following Shaked and Shantikumar [Reference Shaked and Shantikumar18].
Definition 1 Let 
$X$ and 
$Y$ be two nonnegative random variables with distribution functions 
$F_X,F_Y\in \mathcal {F}$, respectively. The random variable 
$X$ (or its distribution 
$F_X$) is said to be smaller than 
$Y$ (or its distribution 
$F_Y$) in the star-shaped order, denoted by 
$X\leq _{*} Y$ (or 
$F_X\leq _\ast F_Y$), if 
$({1}/{x})F_Y^{-1}(F_X(x))$ is increasing with 
$x>0$ (or equivalently, 
${F_Y^{-1}(u)}/{F_X^{-1}(u)}$ is increasing with 
$u\in (0,1)$).
Remark 2 It is easy to verify that the star-shaped order is scale invariant, implying that in case of families of distributions that have a scale parameter, we are able to choose the parameter in a convenient way.
The decision about the star-shaped order often relies on sign variation techniques, as follows from (4.B.2) from Shaked and Shantikumar [Reference Shaked and Shantikumar18]. Expectedly, the sign variation analysis raises technical difficulties, especially when dealing with distributions involving a large number of parameters, such as parallel systems, series systems, or order statistics. Saunders and Moran [Reference Saunders and Moran17] proved a more tractable condition for the star-shaped order to hold, providing a full characterization of such relation within a family of distributions.
Theorem 3 ([Reference Saunders and Moran17])
Let 
$\{ F_a\colon a\in I\subseteq \mathbb {R}\}$ be a family of distributions such that 
$F_a\in \mathcal {F}$ with density 
$f_a$ which does not vanish on any subinterval of its support. Then 
${F_a^{-1}(\alpha )}/{F_a^{-1}(\beta )}$ decreases (resp., increases) with respect to 
$a\in J\subseteq I$, for each fixed 
$\alpha >\beta$, if and only if 
$D(a,x)=({1}/{xf_a(x)})({\partial F_a}/{\partial a})(x)$ increases (resp., decreases) with 
$x$ on the support of 
$F_a$, for every fixed 
$a\in J\subseteq I$.
The simple criterion for the star-shaped relationships that follows is derived as a direct consequence of the monotonicity result presented above.
Theorem 4 Let 
$\{ F_a\colon a\in I\subseteq \mathbb {R}\}$ be a family of distributions as in Theorem 3. Then 
$F_a\leq _\ast F_b$, for every 
$a\leq b$ such that 
$a,b\in J\subseteq I$ if and only if 
$D(a,x)=({1}/{xf_a(x)})({\partial F_a}/{\partial a})(x)$ is decreasing with 
$x$ on the support of 
$F_a$, for every 
$a\in J\subseteq I$.
Proof. Take 
$a,b\in J\subseteq I$, such that 
$a\leq b$. For 
$\alpha \geq \beta$, we have that
$$F_a\leq_\ast F_b\Leftrightarrow\frac{F_b^{{-}1}(\beta)}{F_a^{{-}1}(\beta)}\leq \frac{F_b^{{-}1}(\alpha)}{F_a^{{-}1}(\alpha)} \Leftrightarrow \frac{F_a^{{-}1}(\alpha)}{F_a^{{-}1}(\beta)}\leq \frac{F_b^{{-}1}(\alpha)}{F_b^{{-}1}(\beta)},$$which is equivalent to 
${F_a^{-1}(\alpha )}/{F_a^{-1}(\beta )}$ being increasing with respect to 
$a$. Taking into account Theorem 3, the conclusion follows.
Note that Theorem 4 states a necessary and sufficient condition for the star-shaped order to hold between distributions 
$F_a$ for every 
$a\in J$. However, in general, distributions may depend on more than one parameter, as happens for parallel or series systems with heterogeneous components and, in general, coherent systems. Hence, it is natural to seek for extensions of Theorem 4 to families of distributions indexed by higher dimensional parameters.
To state our result, we need to introduce some notations. Let 
$\mu \in I\subseteq \mathbb {R}^{n}$ and 
$v\in \mathbb {R}^{n}$, and consider 
$\mu +tv$, 
$t\in \mathbb {R}$, the line that goes through 
$\mu$ and has vector 
$v$. We will denote by 
$L_{(\mu ,v)}=\{ \lambda _t\in I\subseteq \mathbb {R}^{n}\colon \lambda _t=\mu +tv, t\in \mathbb {R}\}$. Moreover, given a family of distributions 
$F_{\lambda }$, 
$\nabla F_{\lambda }(x)$ stands for the gradient of 
$F_{\lambda }(x)$ with respect to the parameter 
$\lambda$ and by 
$\langle v,\nabla F_{\lambda }(x)\rangle$, we denote the inner product between 
$v$ and 
$\nabla F_{\lambda }(x)$.
Theorem 5 Let 
$\{ F_{\lambda }\colon \lambda \in I\subseteq \mathbb {R}^{n}\}$ be a family of distributions such that 
$F_{\lambda }\in \mathcal {F}$ and has density function 
$f_{\lambda }$ which does not vanish on any subinterval of its support. Let 
$\mu \in I$, 
$v=(v_1,v_2,\ldots ,v_n)\in \mathbb {R}^{n}$ and 
$J\subseteq I$. Then 
$F_{\lambda _t}\leq _*F_{\lambda _{t'}}$, for every 
$\lambda _t,\lambda _{t'}\in L_{(\mu ,v)}\cap J$ and 
$t\leq t'$, if and only if 
$R(x)={\langle v,\nabla F_{\lambda }(x)\rangle }/{xf_{\lambda }(x)}$ is decreasing with 
$x$ on the support of 
$F_{\lambda }$, for every 
$\lambda \in L_{(\mu ,v)}\cap J$.
Proof. We want to prove that 
$G_{t}\leq _\ast G_{t'}$, for 
$t\leq t'$, where 
$G_t(x)=F_{\lambda _t}(x)$, for every 
$x$ on the support of 
$F_\lambda$. By Theorem 4, this is equivalent to 
$({1}/{xg_t(x)})({\partial G_t}/{\partial t})(x)$ being decreasing with 
$x$, where 
$g_t(x)=G'_t(x)=f_{\lambda }(x)$. Therefore, we may conclude that 
$F_{\lambda _t}\leq _\ast F_{\lambda _{t'}}$ if and only if 
${\langle v,\nabla F_{\lambda }(x)\rangle }/{xf_{\lambda }(x)}$ is decreasing with 
$x$, for every 
$\lambda \in L_{(\mu ,v)}\cap J$.
Remark 6 Note that one could think, of comparing distributions whose parameters belong to some general parametric curve, instead of straight lines, which would lead to an obvious extension of Theorem 5.
In the case of families of distributions with two-dimensional parameters, the following version, using the slope of the line 
$L_{(\mu ,v)}$, is convenient.
Proposition 7 Let 
$\{ F_{\lambda }\colon \lambda \in I\subseteq \mathbb {R}^{2}\}$ be a family of distributions such that 
$F_{\lambda }\in \mathcal {F}$ and has density function 
$f_{\lambda }$ which does not vanish on any subinterval of its support. Let 
$\mu \in I$, 
$v=(v_1,v_2)\in \mathbb {R}^{2}$ and 
$J\subseteq I$. If 
$v_1<0$ (resp., 
$v_1>0$), then 
$F_{\lambda _t}\leq _* F_{\lambda _{t'}}$, for every 
$\lambda _t,\lambda _{t'}\in L_{(\mu ,v)}\cap J$ where 
$t\leq t'$, if and only if 
$Q(x)=({1}/{xf(x)})(({\partial F_{\lambda }}/{\partial \lambda _1})(x)+k({\partial F_{\lambda }}/{\partial \lambda _2})(x))$ is increasing (resp., decreasing) with 
$x$ on the support of 
$F_{\lambda }$, for every 
$\lambda \in L_{(\mu ,v)}\cap J$, where 
$k={v_2}/{v_1}$.
Proof. According to Theorem 5, we have that, 
$F_{\lambda _t}\leq _*F_{\lambda _{t'}}$, for every 
$\lambda _t,\lambda _{t'}\in L_{(\mu ,v)}\cap J$, where 
$t\leq t'$, if and only if 
$R(x)=({1}/{xf_{\lambda }(x)})(({\partial F_{\lambda }}/{\partial \lambda _1})(x)v_1+({\partial F_{\lambda }}/{\partial \lambda _2})(x)v_2)$ is decreasing with 
$x>0$, for every 
$\lambda \in L_{(\mu ,v)}\cap J$. Factorizing 
$R(x)$ by 
$v_1$ and taking into account the sign of 
$v_1$, the conclusion follows.
3. Applications
We now apply the results proved in the previous section to prove comparability, with respect to the star-shaped order, for some models that are popular in reliability theory. Throughout this section, 
$X_1,\ldots ,X_n$ will represent the lifetimes of the components of a complex system. The lifetime of a parallel system is 
$X_{(n)}=\max (X_1,\ldots ,X_n)$, while the lifetime of a series system is given by 
$X_{(1)}=\min (X_1,\ldots ,X_n)$.
3.1. Parallel systems with dependent components
First, we provide a condition for the star-shaped order to hold between parallel systems, for which their lifetime components are dependent and identically distributed. We say that the joint distribution of 
$(X_1,\ldots ,X_n)$ follows an 
$n$-dimensional FGM (Farlie–Gumbel–Morgenstern, cf. [Reference Kotz, Balakrishan and Johnson14]) distribution if
\begin{equation} F_{(X_1,\ldots,X_n)}(x_1,\ldots,x_n)=\prod_{i=1}^{n}F(x_i)\left( 1+\sum_{1\leq j< k\leq n} a_{jk}\bar{F}(x_j)\bar{F}(x_k)\right), \end{equation}where 
$|\sum _{1\leq j< k\leq n} a_{jk}|\leq 1$. Then, the distribution function of 
$X_{(n)}$ is given by
\begin{equation} F_c(x)=F^{n}(x)(1+c\bar{F}^{2}(x)), \end{equation}where 
$c=\sum _{1\leq j< k\leq n} a_{jk}\in [-1,1]$. Note that the constant 
$c$ describes the strength of dependence among the random variables, while its sign reveals the direction of the dependence, that is, if 
$c>0$ 
$(c<0)$, the components are positively (negatively) dependent.
Proposition 8 Let 
$\{ F_c \colon c\in [-1,1]\}$ be a family of distributions defined as in (2). Then 
$F_a\leq _\ast F_b$, whenever 
$-1\leq a< b\leq {n}/{(n+2)}$, if 
$Q_1(x)={F(x)\bar {F}(x)}/{xf(x)}$ is decreasing with 
$x$ on the support of 
$F_c$.
Proof. According to Theorem 4, we need to prove that 
$D(c,x)=({1}/{xf_c(x)})({\partial F_c}/{\partial c})(x)$ is decreasing with 
$x>0$, for 
$-1\leq c\leq {n}/{(n+2)}$. After simplifications, we have 
$D(c,x)=Q_1(x)h(x)$, where 
$h(x)={\bar {F}(x)}/{(n(1+c)-2c(n+1)F(x)+c(n+2)F^{2}(x))}$. Since 
$({\partial F_c}/{\partial c})(x)\geq 0$, we have that 
$D(c,x)\geq 0$. Now, taking into account that obviously 
$Q_1(x)\geq 0$, it follows that 
$h(x)\geq 0$, for 
$x>0$. If 
$c=0$, then 
$h$ is decreasing, and the conclusion follows. Finally, assume that 
$|c| \leq 1$ and 
$c\neq 0$. Given that 
$F$ is increasing and nonnegative, the monotonicity of 
$h$ is the same as the monotonicity of the companion function
$$\tilde{h}(x)=\frac{1-x}{n(1+c)-2c(n+1)x+c(n+2)x^{2}},$$for 
$x\in (0,1)$. Differentiating, it is easily seen that 
$\tilde {h}'$ has the same sign as 
$N(x)=-n(1+c)+2c(n+1)-2c(n+2)x+c(n+2)x^{2}$. When 
$c>0$, 
$N(x)\leq N(0)=c(n+2)-n\leq 0$, while when 
$c<0$, 
$N(x)\leq N(1)\leq 0$. Thus, 
$\tilde {h}'(x)\leq 0$, implying that 
$h$ is decreasing. Taking into account that 
$Q_1$ is a positive decreasing function, the conclusion follows.
Remark 9 Moreover, Proposition 8 implies that, with respect to systems with independent components, negatively dependent components results in faster aging of parallel systems, while positive dependence mean slower aging rates.
Remark 10 The decreasingness assumption on 
$Q_1$ is a technical assumption in Proposition 8. We provide some additional information showing that it is satisfied by Gamma distributions with shape parameter smaller than 
$1$. Note that verifying that 
$Q_1$ is decreasing is equivalent to proving that, for every 
$c>0$, 
$H(x)=\bar {F}(x)F(x)-cxf(x)$ changes sign at most once, and if the sign change occurs, it is in the order “
$+,-$”, as 
$x$ goes from 0 to 
$+\infty$, where 
$F$ is the distribution function of a Gamma distribution with shape parameter 
$\alpha <1$ and scale parameter 
$\lambda =1$. To conclude about the sign variation of 
$H$, we need to study its derivative.
$$H'(x)=\frac{x^{\alpha-1}e^{{-}x}}{\Gamma(\alpha)}({-}F(x)+\bar{F}(x)-c\alpha+cx)= \frac{x^{\alpha-1}e^{{-}x}}{\Gamma(\alpha)} V(x).$$ As 
$x>0$, the sign of 
$H'$ is the same as the sign of 
$V$. Differentiating again, we obtain 
$V'(x)=-2({x^{\alpha -1}e^{-x}}/{\Gamma (\alpha )})+c$, and 
$V''(x)=2({x^{\alpha -2}e^{-x}}/{\Gamma (\alpha )})K(x)$, where 
$K(x)=-\alpha +1+x$. Thus, 
$K$ is increasing. Furthermore, 
$K(0)=-\alpha +1$ and 
$\lim _{x\to +\infty }K(x)=+\infty$. Since 
$\alpha < 1$, it follows that 
$K(x)>0$, for 
$x>0$ therefore, 
$V''(x)>0$, for 
$x>0$ thus, 
$V'$ is increasing. Moreover, 
$\lim _{x\to 0^{+}} V'(x)=-\infty$ and 
$\lim _{x\to +\infty }V'(x)=c$; hence, 
$V'$ has sign variation “
$-,+$”, for 
$x>0$. Consequently, 
$V$ has monotonicity “
$\searrow \nearrow$”. Now, as 
$\lim _{x\to 0^{+}} V(x)=1-c\alpha$, and 
$\lim _{x\to +\infty } V(x)=+\infty$, the sign variation of 
$V$ may be “
$+,-,+$”, “
$-,+$”, or “
$+$”. In the first case, 
$H$ has monotonicity “
$\nearrow \searrow \nearrow$”, so, given that 
$\lim _{x\to 0^{+}}H(x)=0$ and 
$\lim _{x\to +\infty }H(x)=0$, 
$H$ has sign variation “
$+,-$”. If 
$V$ has sign variation “
$-,+$”, 
$H$ will have monotonicity “
$\searrow \nearrow$”, implying that 
$H(x)< 0$, for 
$x>0$. Finally, the case where 
$V(x)\geq 0$ cannot happen, since 
$\lim _{x\to 0^{+}}H(x)=0$ and 
$\lim _{x\to +\infty }H(x)=0$. Thus, 
$H$ changes sign at most once in the order “
$+,-$”, and therefore 
$Q_1$ is increasing.
Remark 11 For other families of distributions, such as Pareto (with support 
$(0,+\infty )$) or Weibull distributions, it is easy to verify that 
$Q_1$ is also decreasing, regardless of the value of the shape parameter, proceeding analogously as in Remark 10.
3.2. Complex systems based on PHR and PRHR models
We now prove some ordering relationships for two models that have received extensive usage when modeling lifetime or survival time data: the PHR model, introduced by Cox [Reference Cox6], and the PRHR model, introduced by Gupta et al. [Reference Gupta, Gupta and Gupta7]. The PRHR model was used, for example, by Tsodikov et al. [Reference Tsodikov, Asselain and Yakovlev19] to describe a stochastic model of spontaneous carcinogenesis, where the progression time of the tumor was modeled by a PRHR model, while Lane et al. [Reference Lane, Looney and Wansley15] modeled bank failure through a PHR model. We recall the definition of these models: a PHR (resp., PRHR) model with baseline distribution 
$F$ has a distribution function satisfying 
$\bar {F}_a(x)=\bar {F}^{a}(x)$ (resp., 
$F_a(x)=F^{a}(x)$), for 
$a>0$. Hence, the PHR and PRHR models introduce a family of distributions depending on one parameter. We first characterize the star-shaped ordering for each of these models as a straightforward consequence of the Saunders and Moran's criterion, Theorem 4.
Proposition 12 Let 
$F\in \mathcal {F}$ with a density that does not vanish on any subinterval of its support, be some baseline distribution. Then 
$F_a\leq _\ast F_b$, for every 
$0< a\leq b$ if and only if 
$\bar {g}(x)={\ln (\bar {F}(x))\bar {F}(x)}/{xf(x)}$ is increasing with 
$x$ on the support of 
$F$, in the case of the PHR model, or 
$g(x)={\ln (F(x))F(x)}/{xf(x)}$ is decreasing with 
$x$ on the support of 
$F$, in the case of the PRHR model.
Proof. Consider the case of PHR model. According to Theorem 4, we have 
$D(a,x)=({1}/{xf_a(x)})({\partial F_a}/{\partial a})(x)= -({1}/{a})\bar {g}(x)$. Hence, since 
$a>0$, the conclusion follows. The PRHR case follows analogously.
Remark 13 The sample maxima and minima from independent, identically distributed random variables are typical examples of PHR and PRHR models.
Remark 14 It can be verified that the function 
$g$ and 
$\bar {g}$ considered in the previous proposition are monotone for several families of distributions popular in reliability or aging models, such as the Gamma, Weibull, or Pareto (with support 
$(0,+\infty )$). For Weibull and Pareto distributions, the monotonicity of both functions does not depend on the value of the shape parameter. However, the same does not happen for Gamma distributions. Although 
$g$ is increasing, for every shape parameter 
$\alpha >0$, the function 
$\bar {g}$ is increasing if 
$\alpha <1$ and decreasing if 
$\alpha >1$. Verifying this may be accomplished as in Remark 10.
3.3. Complex systems with heterogeneous components
Throughout this subsection, we characterize the star-shaped ordering of a few different types of heterogeneous systems, looking both at parallel and series systems. A first model looks at parallel systems with components whose lifetimes distributions are subject to different scale changes. Let 
$F\in \mathcal {F}$ be a distribution function, with density function that does not vanish on any subinterval of its support and consider random variables 
$X_i$ with distribution function 
$F_i(x)=F(\lambda _i x)$, where 
$\lambda _i>0$, for 
$i=1,\ldots ,n$. The distribution function of 
$X_{(n)}$ is
\begin{equation} F_{\lambda}(x)=\prod_{i=1}^{n} F(\lambda_i x),\quad \lambda=(\lambda_1,\ldots,\lambda_n). \end{equation}First, recall below a definition concerning a special class of functions.
Definition 15 ([Reference Karlin8])
A function 
$f\colon \mathbb {R}^{2}\to \mathcal {R}$ is said to be totally positive of order 
$2$ (TP2) if for every real numbers 
$x_1< x_2$ and 
$y_1< y_2$, 
$f(x_1,y_2)f(x_2,y_1)\leq f(x_1,y_1)f(x_2,y_2)$.
Proposition 16 Let 
$\{ F_{\lambda }\colon \lambda \in (0,+\infty )^{n}, \lambda _1<\lambda _2<\cdots <\lambda _n \}$ be a family of distributions defined as in (3). Let 
$\mu , v\in (0,+\infty )^{n}$. If 
$G(a,x)={F(ax)}/{f(ax)}$ is TP2 for 
$a>0$, then 
$F_{\lambda _t}\leq _\ast F_{\lambda _{t'}}$, for every 
$\lambda _t, \lambda _{t'}\in L_{(\mu ,v)}\cap J$ where 
$t\leq t'$ and 
$J= \{ \lambda \in (0,+\infty )^{n}\colon v_i\lambda _j-\lambda _iv_j\leq 0,i< j, i,j=1,\ldots ,n\}$.
Proof. According to Theorem 5, we need to prove that 
$R(x)=({1}/{xf_{\lambda }(x)})\sum _{i=1}^{n} v_i({\partial F_{\lambda }}/{\partial \lambda _i})(x)$ is decreasing with 
$x$, where 
$f_{\lambda }$ is the density function of 
$F_{\lambda }$. We have that
$$\frac{\partial F_{\lambda}}{\partial \lambda}(x)=xf(\lambda_i x)\prod_{\substack{j=1\\ i\neq j}}^{n} F(\lambda_j x) \quad \text{and}\quad f_{\lambda}(x)=\sum_{i=1}^{n} \lambda_if(\lambda_i x)\prod_{\substack{j=1\\ i\neq j}}^{n} F(\lambda_j x).$$ Thus, 
$R(x)={\sum _{i=1}^{n} v_i P_i(x)}/{\sum _{i=1}^{n} \lambda _i P_i(x)}$, where 
$P_i(x)=f(\lambda _i x)\prod _{j=1, i\neq j}^{n} F(\lambda _j x)$. Differentiating 
$R$, we get that the sign of 
$R'$ is the same as the sign of
$$K(x)=\sum_{i=1}^{n-1}\sum_{j=i+1}^{n} (v_i\lambda_j-\lambda_iv_j)(P_i'(x)P_j(x)-P_i(x)P_j'(x)).$$ Since, for 
$i< j$ and 
$i,j=1,\ldots ,n$, we have that 
$v_i\lambda _j-\lambda _i v_j\leq 0$, we need to prove that 
$P_i'(x)P_j(x)-P_i(x)P_j'(x)\geq 0$, for 
$i< j$. The function 
$P_i'(x)P_j(x)-P_i(x)P_j'(x)$ is the numerator of the derivative of 
$L(x)={P_i(x)}/{P_j(x)}={f(\lambda _i x)F(\lambda _j x)}/{f(\lambda _j x)F(\lambda _i x)}$. Given that, for 
$i< j$, 
$\lambda _i<\lambda _j$ and 
$G(a,x)$ is TP2, for every 
$a>0$, it follows that 
$L$ is increasing. Hence, the proof is concluded.
Remark 17 The Gamma distribution is one of the common families in reliability or aging models that verify the TP2 property. In fact, proving that 
$G(a,x)$ is TP2 is the same as proving that 
$K(x)={f(ax)F(bx)}/{f(bx)F(ax)}=({a}/{b})^{\alpha -1} ({e^{bx}F(bx)}/{e^{ax}F(ax)})$ is increasing with 
$x$, for 
$0< a< b$. Again, this is verified in an analogous way as in Remark 10. Other families of distributions that satisfy the TP2 property assumed in Proposition 16 are, for example, Weibull and Pareto (with support 
$(0,+\infty )$).
The following corollary complements the ordering result proved in Kochar and Xu [Reference Kochar and Xu10], where only two types of components with exponential lifetimes were allowed in each parallel system.
Corollary 18 Let 
$\{ F_{\lambda }\colon \lambda \in (0,+\infty )^{n}, \lambda _1<\lambda _2<\cdots <\lambda _n \}$ be a family of distributions defined as in (3), with 
$F(x)=1-e^{-x}$. Let 
$\mu , v\in (0,+\infty )^{n}$. Then 
$F_{\lambda _t}\leq _\ast F_{\lambda _{t'}}$, for every 
$\lambda _t, \lambda _{t'}\in L_{(\mu ,v)}\cap J$ where 
$t\leq t'$ and 
$J= \{ \lambda \in (0,+\infty )^{n}\colon v_i\lambda _j-\lambda _iv_j\leq 0, i,j=1,\ldots ,n, i< j\}$.
Proof. Taking into account Proposition 16, we only need to prove that 
$G(a,x)={F(ax)}/{f(ax)}$ is TP2, for 
$a>0$. But this is equivalent to proving that, for 
$a< b$,
$$K(x)=\frac{f(ax)F(bx)}{f(bx)F(ax)}=\frac{e^{bx}-1}{e^{ax}-1}$$is increasing with 
$x>0$, which is easily seen to be true.
We now have a look into complex systems based on components whose lifetimes follow a PRHR model. Assume that the lifetimes 
$X_i$ have distribution function 
$F_i(x)=F^{\lambda _i}(x)$, for some baseline function 
$F\in \mathcal {F}$ with density 
$f$ that does not vanish on any interval of the support of 
$F$, where 
$\lambda _i>0$, for every 
$i=1,\ldots ,n$. Then, the distribution function of 
$X_{(1)}$ is
\begin{equation} F_{\lambda}(x)=1-\prod_{i=1}^{n} (1-F^{\lambda_i}(x)), \lambda=(\lambda_1,\ldots,\lambda_n). \end{equation}Proposition 19 Let 
$\{ F_{\lambda }\colon \lambda \in (0,+\infty )^{n}, \lambda _1<\lambda _2<\cdots <\lambda _n \}$ be a family of distributions defined as in (4). Let 
$\mu , v\in (0,+\infty )^{n}$. If 
$g(x)={\ln (F(x))F(x)}/{xf(x)}$ is decreasing with 
$x$ on the support of 
$F$, then 
$F_{\lambda _t}\leq _\ast F_{\lambda _{t'}}$, for every 
$\lambda _t, \lambda _{t'}\in L_{(\mu ,v)}\cap J$ where 
$t\leq t'$ and 
$J= \{ \lambda \in (0,+\infty )^{n}\colon v_i\lambda _j-\lambda _iv_j\leq 0,i< j, i,j=1,\ldots ,n\}$.
Proof. Taking into account Theorem 5, we need to prove that 
$R(x)=({1}/{xf_{\lambda }(x)})\sum _{i=1}^{n} v_i ({\partial F_{\lambda }}/{\partial \lambda _i})(x)$ is decreasing for every 
$\lambda \in J$. We have that
$$\frac{\partial F_{\lambda}}{\partial \lambda_i}(x)=\ln(F(x))F^{\lambda_i}(x)\prod_{\substack{j=1\\ i\neq j}}^{n} (1-F^{\lambda_i}(x)),$$and
$$f_{\lambda}(x)=f(x)\sum_{i=1}^{n} \lambda_iF^{\lambda_i-1}(x)\prod_{\substack{j=1\\ i\neq j}}^{n} (1-F^{\lambda_i}(x)).$$ Hence, 
$R(x)=g(x)h(x)$, where 
$h(x)={\sum _{i=1}^{n} v_i P_i(x)}/{\sum _{i=1}^{n} \lambda _i P_i(x)}$, with 
$P_i(x)=F^{\lambda _i-1}(x)\prod _{j=1, i\neq j}^{n} (1-F^{\lambda _j}(x))$. It is easily seen that the first term in 
$R'(x)=g'(x)h(x)+g(x)h'(x)$ is negative. Observe that, since 
$v_i>0$, for every 
$i=1,\ldots ,n$, it follows that 
$h(x)\geq 0$. Hence, to prove that 
$R$ is decreasing, it is enough to establish that 
$h$ is increasing, given that 
$g$ is a negative decreasing function. The sign of 
$h'$ is easily seen to be the same as the sign of 
$\sum _{i=1}^{n-1}\sum _{j=i+1}^{n}(v_i\lambda _j-\lambda _iv_j)(P_i'(x)P_j(x)-P_i(x)P_j'(x))$. Given that 
$v_i\lambda _j-\lambda _iv_j\leq 0$, it remains to prove that 
$P_i'(x)P_j(x)-P_i(x)P_j'(x)\leq 0$, for every 
$i,j=1,\ldots ,n$, 
$i< j$. Observe that 
$P_i'(x)P_j(x)-P_i(x)P_j'(x)$ is the numerator of the derivative of 
$K(x)={P_i(x)}/{P_j(x)}=F^{\lambda _i-\lambda _j}(x)({(1-F^{\lambda _j}(x))}/{(1-F^{\lambda _i}(x))})$. Hence, we need to prove that 
$K$ is decreasing. Given that 
$F$ is increasing and nonnegative, the monotonicity of 
$K$ is the same as the monotonicity of the companion function
$$\tilde{K}(x)=x^{\lambda_i-\lambda_j}\frac{1-x^{\lambda_j}}{1-x^{\lambda_i}} =\frac{x^{-\lambda_j}-1}{x^{-\lambda_i}-1},\quad \text{for } x\in(0,1),$$which is easily seen to be decreasing on 
$(0,1)$.
Assuming the components follow a PHR model, we may derive a similar result about 
$X_{(n)}$. Consider that 
$X_i$ has survival function 
$\bar {F}_i(x)=\bar {F}^{\lambda _i}(x)$, where 
$\lambda _i>0$, for every 
$i=1,\ldots ,n$ and 
$F$ is a baseline distribution as above. Then, the distribution function of 
$X_{(n)}$ is, with 
$\lambda =(\lambda _1,\ldots ,\lambda _n)$,
\begin{equation} F_{\lambda}(x)=\prod_{i=1}^{n} (1-\bar{F}^{\lambda_i}(x)). \end{equation}Proposition 20 Let 
$\{ F_{\lambda }\colon \lambda \in (0,+\infty )^{n}, \lambda _1<\lambda _2<\cdots <\lambda _n \}$ be a family of distributions defined as in (5). Let 
$\mu , v\in (0,+\infty )^{n}$. If 
$\bar {g}(x)={\ln (\bar {F}(x))\bar {F}(x)}/{xf(x)}$ is increasing with 
$x$ on the support of 
$F$, then 
$F_{\lambda _t}\leq _\ast F_{\lambda _{t'}}$, for every 
$\lambda _t, \lambda _{t'}\in L_{(\mu ,v)}\cap J$ where 
$t\leq t'$ and 
$J= \{ \lambda \in (0,+\infty )^{n}\colon v_i\lambda _j-\lambda _iv_j\leq 0, i< j, i,j=1,\ldots ,n\}$.
Proof. The proof is analogous to that of Proposition 19, taking into account that
\begin{align*} & \frac{\partial F_{\lambda}}{\partial \lambda_i}(x)={-}\ln(\bar{F}(x))\bar{F}^{\lambda_i}(x)\prod_{\substack{j=1\\ i\neq j}}^{n} (1-\bar{F}^{\lambda_i}(x)),\\ & f_{\lambda}(x)=f(x)\sum_{i=1}^{n} \lambda_i\bar{F}^{\lambda_i-1}(x)\prod_{\substack{j=1\\ i\neq j}}^{n} (1-\bar{F}^{\lambda_i}(x)), \end{align*}where 
$f$ is the density function of 
$F$, and 
$R(x)=-\bar {g}(x)h(x)$, where 
$h(x)={\sum _{i=1}^{n} v_i P_i(x)}/{\sum _{i=1}^{n} \lambda _i P_i(x)}$, with 
$P_i(x)=\bar {F}^{\lambda _i-1}(x)\prod _{j=1, i\neq j}^{n} (1-\bar {F}^{\lambda _j}(x))$. Thus, we now need to prove that 
$h'(x)\leq 0$, in order to conclude that 
$R'(x)\leq 0$. Since 
$v_i\lambda _j-\lambda _iv_j\leq 0$, for 
$i,j=1,\ldots ,n$, 
$i\neq j$, this follows after proving that 
$P_i'(x)P_j(x)-P_i(x)P_j'(x)\geq 0$, for 
$i< j$, which is easily achieved analogously as in the proof of Proposition 19.
The previous result is an extension of Theorem 3.3 in Kochar and Xu [Reference Kochar and Xu13], where the authors considered one of the systems to be formed by homogeneous components.
3.4. Parallel systems with homogeneous distributions
Arab et al. [Reference Arab, Oliveira and Hadjikyriakou2] proved in their Corollary 7.2, that parallel homogeneous systems with components that have exponential lifetimes age faster as the number of components increases. We may prove this also holds when the components have exponentiated Weibull lifetimes 
$X_1,\ldots ,X_n$, whose distribution function is given by 
$F(x) = (1-e^{-(\lambda x)^{\beta }})^{\alpha }$, for 
$x>0$, where 
$\alpha , \beta > 0$ are shape parameters and 
$\lambda >0$ is a scale parameter. The distribution function of 
$X_{(n)}$ is given by
\begin{equation} F_X(x)=(1-e^{-(\lambda x)^{\beta}})^{a}, \end{equation}where 
$a=\alpha n$, for 
$n\geq 1$. Taking into account Remark 2, we may, without loss of generality, consider 
$\lambda =1$.
Proposition 23 Let 
$\{ F_{(\beta ,a)}\colon a>0, \beta >0\}$ be a family of distributions defined as in (6). Let 
$(\beta ',a'),(\tilde {\beta },\tilde {a})\in (0,+\infty )^{2}$, such that 
$\beta '\geq \tilde {\beta }$, and 
$v=(\tilde {\beta },\tilde {a})-(\beta ',a')$. If 
$\beta '=\tilde {\beta }$, then 
$F_{a'}\leq _\ast F_{\tilde {a}}$, for every 
$a'\geq \tilde {a}$. If 
$\beta '>\tilde {\beta }$, then 
$F_{\lambda _{t}}\leq _\ast F_{\lambda _{t'}}$, for every 
$\lambda _t,\lambda _{t'}\in L_{((\beta ',a'),v)}$ where 
$t\leq t'$, if and only if 
$\lambda _t,\lambda _{t'}\in J=\{ (\beta ,a)\in (0,+\infty )^{2}\colon {(\tilde {a}-a')}/{(\tilde {\beta }-\beta ')}\geq -{a}/{\beta }\}$.
Proof. If 
$\beta '=\tilde {\beta }$, then the family of distributions 
$F_{(\beta ,a)}$ follows a PRHR model. Thus, taking into account Proposition 12, since
$$\frac{\ln(F(x))F(x)}{xf(x)}=\frac{(e^{x}-1)\ln(1-e^{{-}x})}{x}$$is increasing with 
$x>0$, the conclusion follows. Consider now 
$\beta '>\tilde {\beta }$. By Proposition 7, we need to prove that
$$Q(x)=\frac{1}{xf_{(\beta,a)}(x)}\left(\frac{\partial F_{(\beta,a)}}{\partial \beta}(x)+k\frac{\partial F_{(\beta,a)}}{\partial a}(x)\right)$$is increasing with 
$x>0$, for every 
$(\beta ,a)\in L_{((\beta ',a'),v)}$, if and only if 
$k={(\tilde {a}-a')}/{(\tilde {\beta }-\beta ')}\geq -{a}/{\beta }$. We begin by studying the case where 
$k\geq 0$. Taking into account Lemma 8 in Arab and Oliveira [Reference Arab and Oliveira1], we need to prove that, for every 
$c\in \mathbb {R}$, 
$Q(x)-c$ changes sign at most once, as 
$x$ goes from 
$0$ to 
$+\infty$, and if the sign change occurs it is in the order “
$-,+$”. Note 
$Q(x)-c$ and 
$H(x)=(1-e^{-x^{\beta }})^{a-1}P(x)$, where
$$P(x)=(ae^{{-}x^{\beta}}x^{\beta}\ln(x)+k(1-e^{{-}x^{\beta}})\ln(1-e^{{-}x^{\beta}})-ca\beta x^{\beta}e^{{-}x^{\beta}})$$have, for each 
$x>0$, the same sign. Hence, it is enough to characterize the sign of 
$P$. We look at the sign of 
$P'$, that coincides, for 
$x>0$, with the one of 
$V(x)=(-ax^{\beta }\ln (x)+a\ln (x)+{a}/{\beta }+k\ln (1-e^{-x^{\beta }})+k-ca\beta +ca\beta x^{\beta })$. Differentiating 
$V$, we have 
$V'(x)=x^{\beta -1}K(x)$, where 
$K(x)=-a\beta \ln (x)-a+{a}/{x^{\beta }}+{k\beta }/{(e^{x^{\beta }}-1)}+ca\beta ^{2}$. Thus, 
$K'(x)=-{a\beta }/{x}-{a\beta }/{x^{\beta +1}}-{k\beta ^{2} x^{\beta -1}e^{x^{\beta }}}/{(e^{x^{\beta }}-1)^{2}}$. Since 
$k\geq 0$, we have that 
$K'(x)\leq 0$, implying that 
$K$ is decreasing. Given that, 
$\lim _{x\to 0^{+}}K(x)=+\infty$ and 
$\lim _{x\to +\infty }K(x)=-\infty$, it follows that the sign variation of 
$K$, which is the same as sign variation of 
$V'$, is “
$+,-$”. Therefore, 
$V$ has monotonicity “
$\nearrow \searrow$”. Moreover, 
$\lim _{x\to 0^{+}}V(x)=\lim _{x\to +\infty }V(x)=-\infty$, which implies that 
$V$ has sign variation “
$-,+,-$” or “
$-$”. In the first case, we have that 
$P$ has monotonicity “
$\searrow \nearrow \searrow$”. Since 
$\lim _{x\to 0^{+}}P(x)=\lim _{x\to +\infty }P(x)=0$, the sign variation of 
$P$, which coincides with the sign variation of 
$H$, is “
$-,+$”. In the second case, where 
$V(x)\leq 0$, then 
$P$ would be decreasing. But this is impossible, given the behavior of 
$P$ near 
$0$ and at 
$+\infty$. Suppose now that 
$-({-a}/{\beta })\leq k={(\tilde {a}-a')}/{(\tilde {\beta }-\beta ')}\leq 0$. Differentiating 
$Q$, we obtain
$$Q'(x)=\frac{((\beta kx^{\beta}-\beta k)e^{x^{\beta}}+\beta k)\ln(1-e^{{-}x^{\beta}})+(\beta k+a)x^{\beta}}{x^{\beta+1}a\beta}.$$ Some elementary calculus arguments show that 
$Q'(x)\geq 0$, implying that 
$Q$ is increasing with 
$x>0$. If 
$k\leq -{a}/{\beta }$, 
$Q$ cannot be a monotone function, since 
$\lim _{x\to 0^{+}}Q(x)=\lim _{x\to +\infty }Q(x)=+\infty$. Therefore, the proof is concluded.
Remark 24 According to Proposition 23, given 
$(\tilde {\beta },\tilde {a})$ and 
$(\beta ',a')$, we do not only have that 
$F_{(\beta ,a)}\leq _\ast F_{(\tilde {\beta },\tilde {a})}$, for 
$(\beta ,a)\in J\cap L_{((\beta ',a'),v)}$, but we also have that the distributions depending on the parameters in the set 
$J\cap L_{((\beta ',a'),v)}$ are ordered, with respect to the star-shaped order. That is, if 
$k$ is the slope of the line going through 
$(\tilde {\beta },\tilde {a})$ and 
$(\beta ',a')$, then every point 
$(\beta ,a)$ in this line, and above the line 
$-k\beta '$ (condition given by the set 
$J$) defines distributions comparable with each other and with 
$F_{(\tilde {\beta },\tilde {a})}$. However, it does not allow us to decide about star-shaped comparability between 
$F_{(\beta ,a)}$ and 
$F_{(\tilde {\beta },\tilde {a})}$, when 
$(\beta ,a)\notin J$. Nevertheless, if we keep changing the value of 
$k$ (and, therefore, the position of the point 
$(\beta ',a')$), it follows from Proposition 23 that the set of points 
$(\beta ,a)$ for which we have 
$F_{(\beta ,a)}\leq _\ast F_{(\tilde {\beta },\tilde {a})}$, is given by the set 
$\{ (\beta ,a)\in (0,+\infty )^{2}\colon \beta \geq \tilde {\beta }$ and 
$a\geq {\tilde {a}\beta }/{(\tilde {\beta }-2\beta )} \}$.
4. Conclusion
Saunders and Moran [Reference Saunders and Moran17] introduced a criterion to obtain the star-shaped ordering between distributions within a family depending on a single parameter. Although this criterion was used by many authors to exhibit star-shaped comparability within parametric families of distributions, Saunders and Moran's [Reference Saunders and Moran17] result cannot be used for the comparison of complex systems, such as parallel or series systems with heterogeneous components, whose distributions depend on more than one parameter. In this paper, we extend the characterization given by Saunders and Moran [Reference Saunders and Moran17] for the star-shaped order for families of distributions that involve multiple parameters. Based on this new criterion, sufficient conditions are obtained for comparing, with respect to the star-shaped order, parallel and series systems with heterogeneous components. Moreover, the new characterization is used to establish star-shaped comparability between parallel systems with homogeneous components, when the multidimensional parameter move along lines that depend on the scale parameters of the individual lifetime distributions.
Acknowledgments
The authors I.A., P.E.O., and B.S. are partially supported by the Centre for Mathematics of the University of Coimbra – UIDB/00324/2020, funded by the Portuguese Government through FCT/MCTES. B.S. was also supported by FCT, through the grant PD/BD/150459/2019, co-financed by the European Social Fund.
The authors would like to thank the anonymous reviewers for their detailed remarks and extensive references, that helped improving on previous versions of the paper.
