# Are there order statistics for a Gaussian variable raised to a power?

Let $X$ be a random variable with a standard normal distribution. Let $Y = |X|^{2p}$. I am trying to find the distribution for $Y_{(n)}$, i.e., the largest value of $Y$ out of $n$ samples.

I have derived the pdf to be: $$f_{Y_{(n)}} = n \left(\frac{1}{p\sqrt{2\pi}} y^{\frac{1}{2p} - 1} \exp\left(-\frac{1}{2}y^{1/p} \right)\right) \left(\int_0^y \frac{1}{p\sqrt{2\pi}} t^{\frac{1}{2p} - 1} \exp\left(-\frac{1}{2}t^{1/p}\right) \, dt \right)^{n-1}$$

But Mathematica says $EY_{(n)}$ is infinite. Intuitively, I feel that it should be some finite value in terms of p and n. Any ideas?

• What can you tell us about $\rho$? Integer? Positive? Why do you need the 2 in $|X|^{2p}$? – wolfies Mar 1 '16 at 5:06
• But Mathematica says $EY_{(n)}$ is infinite. ///////// Why does the existence of $EY_{(n)}$ concern you? To be precise, why is the existence of the moments relevant to your question ... which is to find the pdf of the sample maximum of $Y$? – wolfies Mar 1 '16 at 5:11
• Let p be a real number greater than 1. The 2 is extraneous, I'll admit. I would like to bound the expected value of this function in terms of p and n, in order to understand the behavior as I increase these parameters. Also, if the expected value is infinite, then the pdf may be wrong, as it doesn't make sense for the max of some Gaussian rvs to be infinite. – Edward Yu Mar 1 '16 at 5:15
• $E(Y_{(n)})\le nE(Y)<\infty$ and for $p>1$ the first inequality is asymptotically sharp. You can easily compute $E(Y)$ in terms of $\Gamma$ function by $u=-x^2/2$ substitution. For $p=1$ you get $C\log n$ instead of $n$ and for $p=1/2$ you get $C \sqrt{\log n}$ instead. – A.S. Mar 1 '16 at 10:09
• Do you happen to have a proof or reference for the inequality $E(Y_{(n)})\le nE(Y)<\infty$ ? – Edward Yu Mar 1 '16 at 17:25

Let $Z \sim N(0,1)$. Then $X = |Z|$ has a half-Normal distribution with pdf $f(x)$:

Let ${X_1, \dots, X_n}$ denote a random sample of size $n$ drawn on $X$, and let $X_{(n)}$ denote the sample maximum. The pdf of $X_{(n)}$, say $g(x)$ is easy to derive with a computer algebra system:

with domain of support $X>0$.

Let $Y = X^\alpha$ where $\alpha >1$. Since $X > 0$, it follows that $Y_{(n)} = {X_{(n)}}^\alpha$.

The pdf of $Y_{(n)}$, say $h(y)$, is then given by the transformation:

defined on the positive real line. All done.

The following diagram plots the pdf of $Y_{(n)}$ when $\alpha = 2$, for different sample sizes $n$:

Monte Carlo check

Here is a quick Monte Carlo check. In the following:

• the blue squiggly curve denotes the simulated empirical pdf of the sample maximum $Y_{(n)}$

• the dashed red curve plots the exact theoretical pdf of $Y_{(n)}$ derived above

when $\alpha = 2$ and $n = 10$:

Looks good.

Notes

1. Erf denotes the: Error function

2. The OrderStat and Transform functions used above are from the mathStatica suite for Mathematica. As disclosure, I should add that I am one of the authors.

Comment:

Trying to visualize what might be going wrong, I simulated this in R for 100,000 samples of size $n = 5$ with $p = 1.5.$ The simulated distribution of $Y_{(5)}$ (left plot) is extremely right-skewed even for small $p$, which may be causing some trouble in Mathematica. However, several runs gave nearly the same values for $E(Y)$ and $SD(Y)$ each time, suggesting that the true values are finite. A histogram (right) shows the distribution of $\log(Y_{(5)})$.

 B = 10^5; p = 1.5; n = 5
y = abs(rnorm(B*n))^(2*p)
DTA = matrix(y, nrow=B)  # B x n: each row a sample of n
y.5 = apply(DTA, 1, max) # vector of B sample maxima
mean(y.5);  sd(y.5)
## 5.437638  # approx E(Y.5)
## 6.134618  # approx SD(Y.5)


Ref: Fragmentary as it is, the Wikipedia page on 'folded normal distribution' may be of some help.

• Hmm, I think you are right that it is a computer error within Mathematica. I tried n = 10 and p = 10 and the integral failed to converge. – Edward Yu Mar 1 '16 at 5:18
• Not sufficiently familiar with Mathematica to help you de-bug. Maybe experiment with huge numbers $M$ for upper limit on the integral. Then there is no question of divergence. – BruceET Mar 1 '16 at 5:24

In this answer I will try to derive an analytic formula for you. I think my answer is the same as yours. Anyway, you can use it to compare your result and check what might have gone wrong.

Let $\Phi$ denote the cumulative distribution function of the standard normal distribution. That is $$\Phi(x) = \int_{-\infty}^{x}\frac{1}{\sqrt{2\pi}}e^{-t^2/2}dt$$ In the following let the $X_i$ be standard normally distributed and $Y_i = |X_i|^{2p}$. Then $$P(Y_{(n)}\le y) = P(Y_1\le y;\cdots;Y_n\le y) = \prod_{i=1}^{n}P(Y_i\le y) =: F_Y(y)^n$$ and thus $$f_{Y_{(n)}}(y) = \frac{d}{dy}F_Y(y)^n = nF_Y(y)^{n-1}f_Y(y).$$ We conclude that if we can determine $F_Y(y)$, then we have the pdf of $Y_{(n)}$. $$F_Y(y) = P(Y\le y) = P(|X|^{2p}\le y) = P(-y^{1/2p}\le X \le y^{1/2p}) = 2\Phi(y^{1/2p})-1$$ From this we see that $$f_Y(y) = \frac{d}{dy}2\Phi(y^{1/2p})-1 = 2f_X(y^{1/2p})\frac{1}{2p}y^{-\frac{2p-1}{2p}}.$$ Finally we get that \begin{align} f_{Y_{(n)}}(y) &= n\left(2\Phi(y^{1/2p})-1\right)^{n-1}2f_X(y^{1/2p})\frac{1}{2p}y^{-\frac{2p-1}{2p}}\\ &= n\left(2\int_{-\infty}^{y^{1/2p}}\frac{1}{\sqrt{2\pi}}e^{-t^2/2}dt-1\right)^{n-1}\left(\frac{1}{\sqrt{2\pi}}e^{-y^{1/p}/2}\right)\frac{1}{p}y^{-\frac{2p-1}{2p}} \end{align}