Variance of the median estimator of normal distribution [closed]

Writing $Z_i = (X_i - \mu)/\sigma$ and noting that the order statistics satisfy $X_{(i)} = \mu + \sigma Z_{(i)}$, it suffices to study the variance of $Z_{(2)}$. Now, it is well-known that the PDF of $Z_{(2)}$ is given by

$$ f_{Z_{(2)}}(x) = 6\Phi(x)(1-\Phi(x))\phi(x), $$

where $\Phi$ and $\phi$ are the CDF and PDF of the standard normal distribution, respectively. Now by noting that

$$ f_{-Z_{(2)}}(x) = f_{Z_{(2)}}(-x) = 6\Phi(-x)(1-\Phi(-x))\phi(-x) = 6(1-\Phi(x))\Phi(x)\phi(x) = f_{Z_{(2)}}(x), $$

it follows that the distribution of $Z_{(2)}$ is symmetric about $0$ and hence $ \mathbf{E}[Z_{(2)}] = 0 $. So, the variance of $Z_{(2)}$ is the same as $\mathbf{E}[Z_{(2)}^2]$. However, by utilizing the formula

$$ \int x \phi(x) \, \mathrm{d}x =- \phi(x) + \mathsf{C}, $$

which can be easily proved by u-substitution, we get

\begin{align*} \mathbf{Var}(Z_{(2)}) &= \mathbf{E}[Z_{(2)}^2] \\ &= \int_{-\infty}^{\infty} 6x^2\Phi(x)(1-\Phi(x))\phi(x) \, \mathrm{d}x \\ &= \int_{-\infty}^{\infty} 6 \bigl( x\Phi(x)(1-\Phi(x)) \bigr)'\phi(x) \, \mathrm{d}x \tag{$\because$ IbP} \\ &= \int_{-\infty}^{\infty} \bigl( \underbrace{6 \Phi(x)(1-\Phi(x)) \phi(x)}_{=f_{Z_{(2)}}(x)} + \underbrace{6x\phi(x)^2}_{\text{odd}} - 12x\phi(x)^2\Phi(x) \bigr) \, \mathrm{d}x \\ &= 1 + 0 - \int_{-\infty}^{\infty} 12x\phi(x)^2\Phi(x) \, \mathrm{d}x. \end{align*}

So it suffices to compute the last integral:

\begin{align*} &\int_{-\infty}^{\infty} 12x\phi(x)^2\Phi(x) \, \mathrm{d}x \\ &= \frac{6}{\sqrt{2}\pi^{3/2}} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} xe^{-x^2-y^2/2} \mathbf{1}_{\{ y \leq x \}} \, \mathrm{d}x \mathrm{d}y \\ &= \frac{6}{\pi^{3/2}} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} xe^{-x^2-y^2} \mathbf{1}_{\{ \sqrt{2}y \leq x \}} \, \mathrm{d}x \mathrm{d}y \tag{$y \mapsto \sqrt{2}y$} \\ &= \frac{6}{\pi^{3/2}} \biggl( \int_{\arctan(1/\sqrt{2})-\pi}^{\arctan(1/\sqrt{2})} \cos \theta \, \mathrm{d}\theta \biggr) \biggl( \int_{0}^{\infty} r^2 e^{-r^2} \, \mathrm{d}r \biggr)\\ &= \frac{6}{\pi^{3/2}} \biggl( \frac{2}{\sqrt{3}} \biggr) \biggl( \frac{\sqrt{\pi}}{4} \biggr) \\ &= \frac{\sqrt{3}}{\pi}. \end{align*}

Therefore

$$ \mathbf{Var}(Z_{(2)}) = 1 - \frac{\sqrt{3}}{\pi} \approx 0.448671. $$

The following is a histogram for the $10^8$ simulated medians of $(Z_1, Z_2, Z_3)$ together with the summary statistics:

Here is the Mathematica code used to generate the above histogram:

(* Generate 10^8 samples *)
data=Table[
    Median[RandomVariate[NormalDistribution[],3]],
    {10^8}
];
(* Draw histogram *)
Labeled[
    Histogram[
        data,
        {Min[data],Max[data],0.2},
        "PDF",
        LabelStyle->Directive["Label",14],
        Epilog->{
        Text[Style[StringTemplate["Mean = `m`\nVar = `v`\nn = `n`"][<|
            "m"->ToString[NumberForm[Mean[data],{Infinity,4}]],
            "v"->ToString[NumberForm[Variance[data],{Infinity,4}]],
            "n"->"10^8"
            |>],14,TextAlignment->Left,"Label"],{Min[data]+0.1,0.62},{-1,1}]
        }
    ],
    {Style["Frequency",14],Style["Simulated median",14]},
    {Left,Bottom},
    RotateLabel->True
]