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.
$\begingroup$
$\endgroup$
11
asked Dec 29, 2021 at 16:36
-
1$\begingroup$ The generalized Pareto distribution has three parameters, but your posterior at first glance would appear to only be on one of them (otherwise you'd have more than two parameters.) Is this correct? If so, how are you dealing with the other two parameters, and which ones are which? $\endgroup$– jbowmanCommented Dec 29, 2021 at 19:05
-
1$\begingroup$ Could you describe the actual model you are using? The generic answer would be to sample, but there may be problem specific solutions. $\endgroup$– TimCommented Dec 29, 2021 at 20:11
-
$\begingroup$ If you have $N$ ordered samples $x_i$, you can use a transformation with the Gaussian cdf $\Phi$: Fit the function $f(p)$ to the datapoints $(\frac{i}{N+1}, \Phi(x_i))$, and then approximate the quantile function $Q(p)$ as $\Phi^{-1}(f(p))$. For instance, you might fit to functions of the form $f(p)=p+ap(1-p)$ with $-1<a<1$, since any such $f$ leads to a quantile function with $Q(0)=-\infty$, $Q(1)=+\infty$, and $Q(p)<Q(p')$ when $p<p'$. $\endgroup$– user225256Commented Dec 29, 2021 at 20:54
-
$\begingroup$ @jbowman my posterior has 2 parameters, the scale and shape; the location parameter is known to be 0. $\endgroup$– Closed Limelike CurvesCommented Dec 29, 2021 at 22:08
-
3$\begingroup$ I should think a 2-dimensional Gaussian quadrature, e.g., jaeckel.org/ANoteOnMultivariateGaussHermiteQuadrature.pdf would work quite well. You could pre-compute almost everything, and the final calculation of the posterior predictive dist'n value at $x$ would just be a weighted average of some number, say 50 or so, generalized Pareto distribution values at $x$. $\endgroup$– jbowmanCommented Dec 30, 2021 at 2:26
|
Show 6 more comments
1 Answer
$\begingroup$
$\endgroup$
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$.
