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?
 A: 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


*

*Erf denotes the: Error function

*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.
A: 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.
A: 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}
