  The stochastic partial differential equation (SPDE) approach is widely used
for modeling large spatial datasets. It is based on representing a Gaussian
random field $u$ on $\mathbb{R}^d$ as the solution of an elliptic SPDE $L^\beta
u = \mathcal{W}$ where $L$ is a second-order differential operator, $2\beta$
(belongs to natural number starting from 1) is a positive parameter that
controls the smoothness of $u$ and $\mathcal{W}$ is Gaussian white noise. A few
approaches have been suggested in the literature to extend the approach to
allow for any smoothness parameter satisfying $\beta>d/4$. Even though those
approaches work well for simulating SPDEs with general smoothness, they are
less suitable for Bayesian inference since they do not provide approximations
which are Gaussian Markov random fields (GMRFs) as in the original SPDE
approach. We address this issue by proposing a new method based on
approximating the covariance operator $L^{-2\beta}$ of the Gaussian field $u$
by a finite element method combined with a rational approximation of the
fractional power. This results in a numerically stable GMRF approximation which
can be combined with the integrated nested Laplace approximation (INLA) method
for fast Bayesian inference. A rigorous convergence analysis of the method is
performed and the accuracy of the method is investigated with simulated data.
Finally, we illustrate the approach and corresponding implementation in the R
package rSPDE via an application to precipitation data which is analyzed by
combining the rSPDE package with the R-INLA software for full Bayesian
inference.
