How do I generate a generalised error random variable? The generalised error distribution (sometimes called the generalised normal distribution) is a generalisation of the normal distribution allowing variation of the kurtosis away from mesokurtosis.  Like the normal distribution, the generalised error distribution is symmetric, and therefore unskewed, but it allows variation of the kurtosis over a wide range of values.  Given the shape parameter $\beta >0$, the standardised form of the generalised error distribution has the density kernel:
$$\text{GenErr}(z|\beta) \propto \exp \Bigg( -\bigg( \frac{\Gamma(3/\beta)}{\Gamma(1/\beta)} \bigg)^{\beta/2} \cdot |z|^\beta \Bigg)
\quad \quad \quad \text{for } z \in \mathbb{R}.$$
The kurtosis of the distribution is a continuous monotonically decreasing function of the parameter $\beta$, and it can take values above the lower bound $\lim_{\beta \rightarrow \infty} \mathbb{Kurt}(Z) = 3/\sqrt{5} \approx 1.341641$ (which is not far above unit kurtosis, which is the lowest kurtosis of any distribution).  In the special case where $\beta = 2$ the generalised error distribution reduces to the standard normal distribution.

My Question: What is the simplest method of generating this random variable?
 A: I do not know the simplest method, but I can put forward one method (which is not particularly simple) to get the ball rolling.  Since the generalised error distribution is symmetric around its mean of zero, this means that if we take a standard normal random variable $Z \sim \text{N}(0,1)$ we should be able to obtain a generalised error random variable as:
$$Z_* = g(Z) \sim \text{GenErr}(\beta),$$
where the transformation $g$ is a monotonically increasing function that is "symmetric around the y-axis" (and differentiable everywhere except possibly at zero) ---i.e., which satisfies $g'(z)>0$ for all $z \neq 0$ and $g(-z) = -g(z)$ for all $z \in \mathbb{R}$.  Since the function $g$ is symmetric about the y-axis, and since this implies that $g(0)=0$, it is sufficient to find the function for values $z>0$.  Thus, without loss of explanatory power, we will look at the transformation only over these values.  Let $F_*$ and $Q_*$ denote the CDF and quantile function of the generalised error distribution.  Using the monotonic transformation of interest, we have:
$$\begin{equation} \begin{aligned}
\Phi(z)
= \mathbb{P}(0 < Z \leqslant z) 
= \mathbb{P}(0 < Z_* \leqslant g(z)) 
= F_*(g(z)|\beta). \\[6pt]
\end{aligned} \end{equation}$$
We therefore have:
$$g(z) = Q_*(\Phi(z)|\beta)
\quad \quad \quad
\text{for all } z > 0.$$
There are probably simpler methods to generate this random variable
