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

Question

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?

Answer 1

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.

Answer 2

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.

Answer 3

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}