Approximate Posterior Predictive Quantiles with Numerical Methods I have a posterior function which is easy to approximate using numerical methods (the posterior has only 2 parameters, and is approximately Gaussian because of the large sample). However, I need to evaluate the quantile function of the posterior predictive at every sample point, and the posterior predictive is definitely not Gaussian (the variables are drawn from a generalized Pareto distribution). Is there a fast way to get this numerically (i.e. without resorting to random sampling)? I'm going to have to repeat this operation possibly millions of times whenever I call a function, so it's important to do it quickly.

A 2-dimensional Gaussian quadrature, as described in this paper, may well do the job for you.
I'll briefly describe the procedure for the case where your Gaussian posterior distribution has a low correlation between the two parameters, then describe generalizing the procedure to the alternative case.
To concretize, let's assume you don't numerically integrate beyond $\pm 3$ standard deviations from the posterior means ($3$ being an arbitrary number.)  The low correlation means we can integrate over the square $(-3, -3)$ to $(3,3)$ using a grid of points without wasting a lot of points on extremely low probability regions.  Gaussian quadrature works on the interval $(-1,1)$, and remapping one to the other is of course trivial.  An n-point 1D Gaussian quadrature will have (weight, abscissa) values as shown on this page, but of course numerical integration routines will have them internally as well.  To extend to 2 dimensions, we simply take the values on the grid, e.g., a 5x5 Gaussian quadrature would have a point at $(0.9061798..., 0.9061798...)$, another at $(0.9061798..., 0.5384693...)$ and so forth, with weights $w_{ij}$ equal to the product of the one-dimensional weights  $w_i$ and $w_j$.
We can also calculate the posterior probabilities $p(\sigma_i, \alpha_j)$ at each point on the grid, and multiply by the weights to get probability-adjusted weights $w'_{ij}$. The final calculation of the posterior predictive distribution is just a weighted sum:
$$p(x) = \sum_{i,j}p_{gP}(x; \sigma_i, \alpha_j)w'_{ij}$$
The weights don't change with different values of $x$, so this computation is about as fast as it can get.
To check for convergence, one typically increases the number of points in the grid and sees whether the values are about the same.  More efficient use of points can be made via Gauss-Kronrod integration; G-K points for, say, $n=10$ reuse the points for $n=5$, so computations aren't wasted.
Gaussian quadrature of degree $n$ is exact for polynomials of degree $2n-1$; in the multidimensional case, this only applies to polynomials of each of the variables separately, i.e., $n=3$ does not imply that, say, $x_1^3x_2^2$ will be integrated exactly.
Generalizing to correlated variables is described in the paper as well, but the basic approach is about what you would guess; find a transform that renders the transformed variables approximately independent, then integrate over the transformed variables instead.  This can be done in the Gaussian case via spectral decomposition:

*

*Decompose the covariance matrix $\Sigma = S\Lambda S^T$, where $\Lambda$ is a diagonal matrix.

*Form $A = S\sqrt\Lambda$

*Integrate over the two independent Gaussian variates $z = (z_1, z_2)$, transforming them to $(\sigma_i, \alpha_j)$ by premultiplying by $A$.

