Let $\boldsymbol{X} = (X\_1,X\_2)^{\rm T}\sim \mathcal{N}\_2(\boldsymbol{\mu}, \mathrm{\Sigma})$, where
\begin{eqnarray\*}
\boldsymbol{\mu} = (\mu\_1, \mu\_2)^{\rm T}& = &(\sqrt{\xi\_1\xi\_2/(\xi\_1+\xi\_2)}, 0)^{\rm T}\\
\mathrm{\Sigma} & = &\begin{pmatrix} 1 & -\rho\\
-\rho & 1\end{pmatrix}\\
\rho & = & \sqrt{\xi\_1\xi\_3/(\xi\_1+\xi\_2)(\xi\_2+\xi\_3)}.
\end{eqnarray\*}
It is given that $\xi\_1\leq\xi\_2\leq\xi\_3$, $\xi\_i\geq 0$ and $\sum\_{i=1}^3\xi\_i = 1$. I have a function
\begin{equation\*}
 \pi(\boldsymbol{\mu};\boldsymbol{\xi}) = 1-\mathbb{P}(\boldsymbol{X}\leq \boldsymbol{t}) = 1-\mathbb{P}(X\_1\leq t \cap X\_2\leq t),
\end{equation\*}
where $t>1$. Under the constraints $\xi\_1\leq\xi\_2\leq\xi\_3$, $\xi\_i\geq 0$ and $\sum\_{i=1}^3\xi\_i = 1$, numerically we are getting the maximum of $\pi(\boldsymbol{\mu};\boldsymbol{\xi})$ where the non-zero mean $\mu\_1$ is maximized. This occurs when $\xi\_1 = \xi\_2 = \xi\_3 = 1/3$.
Can we use some kind of stochastic ordering arguments to prove the result theoretically?



A problem related to bivariate normal stochastic order
$\newcommand{\der}{\mathrm{der}}\newcommand{\tdert}{\mathrm{dert}}\newcommand{\erf}{\operatorname{erf}}\newcommand{\eqs}{\overset{\text{sign}}=}\newcommand{\tder}{\widetilde\der}$The statement about the maximum is true, but the proof has hardly anything to do with stochastic ordering arguments. Rather, the problem is reduced to certain problems of real algebraic geometry -- that is, to solving (rather complicated) systems of algebraic inequalities over $\mathbb R$.


Let $x:=\xi\_1$ and $y:=\xi\_2$, so that $\xi\_3=1-x-y$, $0\le x\le y\le1-x-y$,
$$\mu\_1=M(x,y):=\sqrt{\frac{x y}{x+y}},$$
$$\rho=-R(x,y),\quad R(x,y):=\sqrt{\frac{x(1-x-y)}{(x+y)(1-x)}}.$$
Let also
$$p\_t(x,y):=1-P(X\_1\le t,X\_2\le t)$$
and
\begin{equation\*}
 G:=\{(x,y)\colon 0\le x\le y\le1-x-y\}=\{(x,y)\colon 0\le x\le1/3,x\le y\le(1-x)/2\}. \tag{-1}
\end{equation\*}
We want to show that
\begin{equation\*}
 p\_t(x,y)\le p\_t(1/3,1/3) \tag{0}
\end{equation\*}
for all $(x,y)\in G$ and all real $t\ge1$.


Note that
\begin{equation\*}
 p\_t(x,y):=1-P(X\_1\le t,X\_2\le t)=P(M(x,y),-R(x,y)), \tag{1}
\end{equation\*}
where
\begin{equation\*}
 P(m,r):=1-\int\_{-\infty}^t du \int\_{-\infty}^t dv\, f\_{m,r}(u,v) 
\end{equation\*}
and $f\_{m,r}$ is the density function of the bivariate normal distribution with means $m,0$, variances $1,1$, and correlation $r$.


The key is Plackett's observation ([formula (3)](https://www.jstor.org/stable/2332716)) that
\begin{equation\*}
 D\_r f\_{m,r}(u,v)=D\_v D\_u f\_{m,r}(u,v),
\end{equation\*}
where $D\_w$ denote the partial derivative with respect to a variable $w$. It follows that
\begin{equation\*}
\begin{aligned}
 D\_r P(m,r)&:=-\int\_{-\infty}^t du \int\_{-\infty}^t dv \, D\_r f\_{m,r}(u,v) \\ 
 &=-\int\_{-\infty}^t du \int\_{-\infty}^t dv \, D\_v D\_u f\_{m,r}(u,v) \\ 
 &=-f\_{m,r}(t,t)<0. 
\end{aligned} 
\tag{2}
\end{equation\*}
Next,
\begin{equation\*}
\begin{aligned}
P(m,r)&=1-\int\_{-\infty}^t du \int\_{-\infty}^t dv\, f\_{0,r}(u-m,v) \\ 
&=1-\int\_{-\infty}^{t-m} dw \int\_{-\infty}^t dv\, f\_{0,r}(w,v) 
\end{aligned} 
\end{equation\*}
and hence
\begin{equation\*}
\begin{aligned}
 D\_m P(m,r)=\int\_{-\infty}^t dv\, f\_{0,r}(t-m,v)=\int\_{-\infty}^t dv\, f\_{m,r}(t,v)>0. 
\end{aligned} 
\tag{3}
\end{equation\*}


Also, for $(x,y)\in G\setminus\{(0,0)\}$,
\begin{equation\*}
 (D\_x M(x,y),D\_y M(x,y))=\frac{(y^2,x^2)}{(x + y)^2},
\end{equation\*}
\begin{equation\*}
 (D\_x R(x,y),D\_y R(x,y))=\frac{((1-2 x-y)y,-(1-x) x)}{(1-x)^2 (x + y)^2},
\end{equation\*}
so that $D\_x M(x,y)\ge0$ and $D\_x R(x,y)\ge0$.


So, by (1), a chain rule of differentiation, and the inequalities in (2) and (3),
\begin{equation\*}
 D\_x p\_t(x,y)=D\_m P(M(x,y),-R(x,y))D\_x M(x,y)-D\_r P(M(x,y),-R(x,y))D\_x R(x,y)\ge0. 
\end{equation\*}


So, the maximum of $p\_t(x,y)$ over $(x,y)\in G$ occurs at one of right boundaries of $G$, that is, where either $y=x$ and $y=(1-x)/2$.


Next, for $x\in[0,1/3]$ (cf. (-1)),
\begin{equation\*}
 \frac d{dx}\,M(x,(1-x)/2)=\frac{1-2x-x^2}{(1+x)^2}>0,
\end{equation\*}
\begin{equation\*}
 \frac d{dx}\,R(x,(1-x)/2)=\frac1{(1+x)^2}>0. 
\end{equation\*}
So, by (1) and the inequalities in (2) and (3), $p\_t(x,(1-x)/2)$ is increasing in $x\in[0,1/3]$.


It remains to show that
\begin{equation\*}
 \der(t,x)\overset{\text{(?)}}\ge0 \text{ for $x\in[0,1/3]$ and $t\ge1$,} \tag{4}
\end{equation\*}
where
\begin{equation\*}
\begin{aligned}
 \der(t,x)&:=8 \sqrt{\pi } \sqrt{x}\, e^{\left(\sqrt{2} \sqrt{x}-2 t\right)^2/8}\, D\_x p\_t(x,x) \\ 
 &=1+\erf\left(v(t,x)/2\right)-\frac{2 
 }{\sqrt{\pi } (1-x)}\,\sqrt{\frac{x}{1-2 x}}\,e^{u(t,x)/4}, \\ 
 u(t,x)&:=t^2 \left(-4 \sqrt{2} \sqrt{\frac{1-2x}{1-x}} (1-x)+8 x-6\right) \\ 
 &+2 t \left(-2
 \sqrt{2} x^{3/2}+2 (1-x) \sqrt{\frac{1-2x}{1-x}} \sqrt{x}+\sqrt{2}
 \sqrt{x}\right)+x (2 x-1), \\ 
 v(t,x)&:=\sqrt{1-x} \left(\sqrt{\frac{1-2x}{1-x}} \left(\sqrt{2} t-\sqrt{x}\right)+2
 t\right). 
\end{aligned} 
\end{equation\*}


By a chain rule of differentiation,
\begin{equation\*}
 D\_x p\_t(x,x)=D\_m P(M(x),-R(x))M'(x)-D\_r P(M(x),-R(x))R'(x), 
\end{equation\*}
where
\begin{equation\*}
 M(x):=M(x,x),\quad R(x):=R(x,x).
\end{equation\*}
Using (2) and (3), we get
\begin{equation\*}
 \der(t,x)=A\_1(x)e^{A\_2(t,x)}+1+\erf(A\_3(t,x));\tag{4.5}
\end{equation\*}
here and in what follows, $A\_1,A\_2,\dots$ are certain algebraic functions and $\erf$ is the error function.
So,
\begin{equation\*}
 D\_t\der(t,x)
 =A\_4(t,x)e^{A\_2(t,x)}+A\_5(t,x)e^{-A\_3(t,x)^2/2}. \tag{5}
\end{equation\*}
thus getting rid of the function $\erf$. Moreover,
\begin{equation\*}
 A\_5(t,x)>0
\end{equation\*}
(for $t$ and $x$ as in (4)).


Letting now
\begin{equation\*}
 \tdert(t,x):=\frac{D\_t\der(t,x)}{A\_5(t,x)e^{-A\_3(t,x)^2/2}},
\end{equation\*}
we get
\begin{equation\*}
 D\_t\der(t,x)\eqs\tdert(t,x)=A\_6(t,x)e^{A\_7(t,x)}+1, \tag{6}
\end{equation\*}
where $a\eqs b$ means that $a$ and $b$ are of the same sign; this reduces the two exponential expressions $e^{A\_2(t,x)}$ and $e^{-A\_3(t,x)^2/2}$ in (5) to one such expression.
Now we have
\begin{multline\*}
 D\_t\tdert(t,x)=e^{A\_7(t,x)}(1-x) 
 \sqrt{\frac{(1-2x)(1-x)}{x}} \Big(2+\sqrt2\,\sqrt{\frac{1-2x}{1-x}}\,\Big) \\ 
\times \big(3-4 x+2 \sqrt{2} \sqrt{(1-2 x) (1-x)}\,\big)>0
\end{multline\*}
(for $x\in(0,1/3]$).  

Also,  

\begin{equation\*}
 \tdert(0,x)
 =\frac{(1-2 x) \left(2-2 x+\sqrt{2} \sqrt{(1-2 x) (1-x)}\right)}{\left(2+\sqrt{2}
 \sqrt{2+\dfrac{1-2 x}{1-x}}\right) (1-x)^2}>0.
\end{equation\*}
So, $\tdert(t,x)>0$ for all $t\ge0$ and all $x\in(0,1/3]$.


So, by (6), it is enough to show that
\begin{equation\*}
 \der\_1(x):=\der(1,x)\overset{\text{(?)}}\ge0 \tag{7}
\end{equation\*}
for all $x\in(0,1/3]$.


Using (4.5) and some algebra, we get
\begin{equation\*}
 \der\_1(x)=1-\frac{2/\sqrt\pi}{1-x}\,\sqrt{\frac x{1-2x}}\,e^{p(x)/4}+\erf(q(x)),
\end{equation\*}
where
\begin{equation\*}
 p(x):=-6+2 \sqrt{2} \sqrt{x}+7 x-4 \sqrt{2} x^{3/2}+2 x^2+\sqrt{\frac{1-2 x}{1-x}} \left(-4 \sqrt{2}+4
 \sqrt{x}+4 \sqrt{2} x-4 x^{3/2}\right),
\end{equation\*}
\begin{equation\*}
 q(x):=\frac{1}{2} \sqrt{1-2 x} \left(\sqrt{2}-\sqrt{x}\right)+\sqrt{1-x}.
\end{equation\*}


Let
\begin{equation\*}
\der\_{11}(x):=\der\_1'(x)2\sqrt\pi\,
\frac{(1-2 x) (1-x)^2 \sqrt{(1-2x)x} }
{p\_1(x)e^{p(x)/4}},
\end{equation\*}
where
\begin{equation\*}
\begin{aligned}
 p\_1(x)&:=(1-2x)(1-x) x p'(x)+2+2x-8 x^2 \\ 
 &=2 + 2 x - 8 x^2 + 
 x (1 - 2 x) ((7 + \sqrt2/\sqrt{x} - 6 \sqrt2 \sqrt{x} + 4 x) (1 - 
 x) +p\_2(x)),
\end{aligned}
\end{equation\*}
\begin{equation\*}
 p\_2(x):=2 (\sqrt2 - \sqrt{x})/
 s(x) + (4 \sqrt2 + 2/\sqrt{x} - 6 \sqrt{x}) (1 - x) s(x),
\end{equation\*}
\begin{equation\*}
 s(x):=\sqrt{\frac{1-2 x}{1-x}}. 
\end{equation\*}
It is elementary to show that $p\_2(x)\ge5$ (for $x\in(0,1/3]$), and hence
\begin{equation\*}
\begin{aligned}
 p\_1(x)&\ge2 + 2 x - 8 x^2 + 
 x (1 - 2 x) ((7 + \sqrt2/\sqrt{x} - 6 \sqrt2 \sqrt{x} + 4 x) (1 - 
 x) +5)>0. 
\end{aligned}
\end{equation\*}
So,
\begin{equation\*}
 \der\_{11}(x)\eqs\der\_1'(x). 
\end{equation\*}
Next, let
\begin{equation\*}
 \der\_{111}(x):=\der\_{11}'(x)
\sqrt{(1-2x)x} \, e^{p(x)/4+q(x)^2} \, p\_1(x)^2 =\tder\_{111}(\sqrt x), 
\end{equation\*}
where
\begin{equation\*}
\begin{aligned}
 \tder\_{111}(z)&:=
 -3 \sqrt{2}+50 z+56 \sqrt{2} z^2-282 z^3-61 \sqrt{2} z^4 \\ 
 &+392 z^5-284 \sqrt{2} z^6+248 z^7+540
 \sqrt{2} z^8 \\ 
 &-832 z^9-160 \sqrt{2} z^{10}+384 z^{11}-64 \sqrt{2} z^{12} \\ 
 &-2 \sqrt{\frac{1-2z^2}{1-z^2}} \\ 
 &\times\big(1-6 \sqrt{2} z-32 z^2+48 \sqrt{2} z^3 +33 z^4 \\ 
 & -106 \sqrt{2} z^5 +136 z^6+48\sqrt{2} z^7-230 z^8 \\ 
 &+80 \sqrt{2} z^9+48 z^{10}-64 \sqrt{2} z^{11}+32 z^{12}\big). 
\end{aligned} 
\end{equation\*}
The function $\tder\_{111}$ switches its sign only once on the interval $(0,1/\sqrt3]$, from $-$ to $+$, and hence $\der\_{11}'$ has the same sign pattern on the interval $(0,1/3]$. So, for some $x\_\*\in(0,1/3)$, the function $\der\_{11}$ is decreasing on $(0,x\_\*]$ and increasing on $[x\_\*,1/3]$. Also, $\der\_{11}(0+)=-3/2<0$ and $\der\_{11}(1/3)=-1.11\ldots<0$. So, $\der\_{11}<0$ on $(0,1/3]$. So, $\der\_1$ is decreasing on $(0,1/3]$, to $\der\_1(0)=1.31\ldots>0$.


Thus, (7) is proved. $\quad\Box$


