| Title: | Spatio-Temporal Modeling of Large Data Using a Spectral SPDE Approach |
|---|---|
| Description: | Functionality for spatio-temporal modeling of large data sets is provided. A Gaussian process in space and time is defined through a stochastic partial differential equation (SPDE). The SPDE is solved in the spectral space, and after discretizing in time and space, a linear Gaussian state space model is obtained. When doing inference, the main computational difficulty consists in evaluating the likelihood and in sampling from the full conditional of the spectral coefficients, or equivalently, the latent space-time process. In comparison to the traditional approach of using a spatio-temporal covariance function, the spectral SPDE approach is computationally advantageous. See Sigrist, Kuensch, and Stahel (2015) <doi:10.1111/rssb.12061> for more information on the methodology. This package aims at providing tools for two different modeling approaches. First, the SPDE based spatio-temporal model can be used as a component in a customized hierarchical Bayesian model (HBM). The functions of the package then provide parameterizations of the process part of the model as well as computationally efficient algorithms needed for doing inference with the HBM. Alternatively, the adaptive MCMC algorithm implemented in the package can be used as an algorithm for doing inference without any additional modeling. The MCMC algorithm supports data that follow a Gaussian or a censored distribution with point mass at zero. Covariates can be included in the model through a regression term. |
| Authors: | Fabio Sigrist, Hans R. Kuensch, Werner A. Stahel |
| Maintainer: | Fabio Sigrist <[email protected]> |
| License: | GPL-2 |
| Version: | 1.7.5 |
| Built: | 2024-11-25 06:48:42 UTC |
| Source: | CRAN |
This is an R package for spatio-temporal modeling of large data sets. It provides tools for modeling of Gaussian processes in space and time defined through a stochastic partial differential equation (SPDE). The SPDE is solved in the spectral space, and after discretizing in time and space, a linear Gaussian state space model is obtained. When doing inference, the main computational difficulty consists in evaluating the likelihood and in sampling from the full conditional of the spectral coefficients, or equivalently, the latent space-time process. In comparison to the traditional approach of using a spatio-temporal covariance function, the spectral SPDE approach is computationally advantageous. This package aims at providing tools for two different modeling approaches. First, the SPDE based spatio-temporal model can be used as a component in a customized hierarchical Bayesian model (HBM). The functions of the package then provide parametrizations of the process part of the model as well as computationally efficient algorithms needed for doing inference with the HBM. Alternatively, the adaptive MCMC algorithm implemented in the package can be used as an algorithm for doing inference without any additional modeling. The MCMC algorithm supports data that follow a Gaussian or a censored distribution with point mass at zero. Covariates can be included in the model through a regression term.
Fabio Sigrist, Hans R. Kuensch, Werner A. Stahel
Maintainer: Fabio Sigrist <[email protected]>
Fabio Sigrist, Hans R. K\"unsch, and Werner A. Stahel, "Stochastic Partial Differential Equation Based Modeling of Large Space-Time Data Sets", Journal of the Royal Statistical Society: Series B, Volume 77, Issue 1, 2015, pages 3-33
Fabio Sigrist, Hans R. Kuensch, Werner A. Stahel, "spate: An R Package for Spatio-Temporal Modeling with a Stochastic Advection-Diffusion Process.", Journal of Statistical Software, Volume 63, Number 14, 2015, pages 1-23, URL http://www.jstatsoft.org/v63/i14/
Function that returns the color scale for 'image()'. This function is a simplification of the function 'tim.colors()' from the 'fields' package.
cols()cols()
A vector with colors.
Fabio Sigrist
Fields Development Team (2006). fields: Tools for Spatial Data. National Center for Atmospheric Research, Boulder, CO. URL http://www.cgd.ucar.edu/Software/Fields.
Forward Filtering Backward Sampling algorithm for sampling from the
joint full conditional of the hidden state of a linear, Gaussian state
space model. To be more specific, one samples from where
is specified through
and
ffbs(y, lp, G, Sigma, H, Omega, N = dim(y)[2],T = dim(y)[1], NF = dim(G)[1], lglk = FALSE, BwSp = TRUE, filt = FALSE)ffbs(y, lp, G, Sigma, H, Omega, N = dim(y)[2],T = dim(y)[1], NF = dim(G)[1], lglk = FALSE, BwSp = TRUE, filt = FALSE)
y |
Observed data in an T x N matrix with columns and rows corresponding to time and space, respectively. |
lp |
Mean (linear predictor) in an T x N matrix with columns and rows corresponding to time and space, respectively. |
G |
Propagator matrix of the latent process |
Sigma |
Innovation covariance matrix of the latent process |
H |
Observation matrix relating y to |
Omega |
Covariance matrix of the observation error |
N |
Number of points in space. |
T |
Number of points in time. |
NF |
Dimension of the latent process |
lglk |
Logical; if 'TRUE' the value of the log-likelihood is returned as well. |
BwSp |
Logical; if 'TRUE' a sample from the full conditional of |
filt |
Logical; if 'TRUE' the filtered values for |
In the context of the SPDE, are the Fourier coefficients.
A list with entries (depending on whether 'lglk', 'BwSp', 'filt' are 'TRUE' or 'FALSE'):
simAlpha |
A T x N matrix with a sample from the full conditional
of latent process |
ll |
The evaluated log-likelihood, |
mtt |
A T x N matrix with the mean of the full conditional of latent process |
Fabio Sigrist
Forward Filtering Backward Sampling algorithm for sampling from the
joint full conditional of the coefficients and for
evaluation of the log-likelihood.
ffbs.spectral(w=NULL,wFT=NULL,spec=NULL,Gvec=NULL,tau2=NULL,par=NULL,n,T,lglk=FALSE, BwSp=TRUE,NF=n*n,indCos=(1:((n*n-4)/2)*2+3),ns=4,nu=1,dt=1)ffbs.spectral(w=NULL,wFT=NULL,spec=NULL,Gvec=NULL,tau2=NULL,par=NULL,n,T,lglk=FALSE, BwSp=TRUE,NF=n*n,indCos=(1:((n*n-4)/2)*2+3),ns=4,nu=1,dt=1)
w |
Observed data or latent process w (depending on which data model is used) in an T x n*n matrix with columns and rows (points on a grid stacked into a vector) corresponding to time and space, respectively. |
wFT |
Vector of length T*n*n containing the real Fourier transform of 'w'. |
spec |
Spectrum of the innovations |
Gvec |
The propagator matrix G in vector format obtained from 'get.G.vec'. If 'Gvec' is not given, it is constructed based on 'par'. |
tau2 |
Measurement error variance tau2. If 'NULL'; tau2=par[9]. |
par |
Vector of parameters for the SPDE in the following order: rho_0, sigma^2, zeta, rho_1, gamma, alpha, mu_x, mu_y, tau^2. If 'spec' and 'Gvec' are given, 'par' will not be used. |
n |
Number of grid points on each axis. n*n is the total number of spatial points. |
T |
Number of points in time. |
lglk |
Logical; if 'TRUE' the value of the log-likelihood is returned as well. |
BwSp |
Logical; if 'TRUE' a sample from the full conditional of |
NF |
Number of Fourier functions used. |
indCos |
Vector of integers indicating the position cosine terms in the 1:NF real Fourier functions. The first 'ns' cosine wavenumbers in 'wave' are not included in 'indCos'. |
ns |
Number of real Fourier functions that have only a cosine and no sine term. 'ns' is maximal 4. |
nu |
Smoothness parameter of the Matern covariance function for the innovations. By default this equals 1 corresponding to the Whittle covariance function. |
dt |
Temporal lag between two time points. By default, this equals 1. |
A list with entries (depending on whether 'lglk' are 'BwSp' are 'TRUE' or 'FALSE'):
simAlpha |
A T x n*n matrix with a sample from the full conditional
of latent process |
ll |
The evaluated log-likelihood, |
Fabio Sigrist
Function for obtaining the spectral propagator matrix G of the vector autoregressive model for the Fourier coefficients.
get.propagator(wave, indCos, zeta, rho1, gamma, alpha, muX, muY, dt = 1, ns=4)get.propagator(wave, indCos, zeta, rho1, gamma, alpha, muX, muY, dt = 1, ns=4)
wave |
Spatial wavenumbers. |
indCos |
Vector of integers indicating the position of columns in 'wave' of wavenumbers of cosine terms. |
zeta |
Damping parameter |
rho1 |
Range parameter of the diffusion term |
gamma |
Parameter that determines the amount of anisotropy in the diffusion term |
alpha |
Parameter that determines the direction of anisotropy in the diffusion term |
muX |
X component of the drift vector. |
muY |
Y component of the drift vector. |
dt |
Temporal lag between two time points. By default, this equals 1. |
ns |
Number of real Fourier functions that have only a cosine and no sine term. 'ns' is maximal 4. |
Propagator matrix G.
Fabio Sigrist
##For illustration, four grid points on each axis n <- 4 wave <- wave.numbers(n) G <- get.propagator(wave=wave$wave,indCos=wave$indCos,zeta=0.5, rho1=0.1, gamma=2, alpha=pi/4, muX=0.2, muY=-0.15,dt=1,ns=4) round(G,digits=2) ## View(round(G,digits=2)) ##An example n <- 50 spec <- matern.spec(wave=spate.init(n=n,T=1)$wave,n=n,rho0=0.05,sigma2=1,norm=TRUE) alphat <- sqrt(spec)*rnorm(n*n) ##Propagate initial state wave <- wave.numbers(n) G <- get.propagator(wave=wave$wave,indCos=wave$indCos,zeta=0.5, rho1=0.02, gamma=2, alpha=pi/4, muX=0.2, muY=0.2,dt=1,ns=4) alphat1 <- G%*%alphat opar <- par(no.readonly = TRUE) par(mfrow=c(1,2)) image(1:n,1:n,matrix(real.fft(alphat,n=n,inv=FALSE),nrow=n),main="Whittle field",xlab="",ylab="",col=cols()) image(1:n,1:n,matrix(real.fft(alphat1,n=n,inv=FALSE),nrow=n),main="Propagated field",xlab="",ylab="",col=cols()) par(opar) # Reset par() settings##For illustration, four grid points on each axis n <- 4 wave <- wave.numbers(n) G <- get.propagator(wave=wave$wave,indCos=wave$indCos,zeta=0.5, rho1=0.1, gamma=2, alpha=pi/4, muX=0.2, muY=-0.15,dt=1,ns=4) round(G,digits=2) ## View(round(G,digits=2)) ##An example n <- 50 spec <- matern.spec(wave=spate.init(n=n,T=1)$wave,n=n,rho0=0.05,sigma2=1,norm=TRUE) alphat <- sqrt(spec)*rnorm(n*n) ##Propagate initial state wave <- wave.numbers(n) G <- get.propagator(wave=wave$wave,indCos=wave$indCos,zeta=0.5, rho1=0.02, gamma=2, alpha=pi/4, muX=0.2, muY=0.2,dt=1,ns=4) alphat1 <- G%*%alphat opar <- par(no.readonly = TRUE) par(mfrow=c(1,2)) image(1:n,1:n,matrix(real.fft(alphat,n=n,inv=FALSE),nrow=n),main="Whittle field",xlab="",ylab="",col=cols()) image(1:n,1:n,matrix(real.fft(alphat1,n=n,inv=FALSE),nrow=n),main="Propagated field",xlab="",ylab="",col=cols()) par(opar) # Reset par() settings
Function for obtaining the spectral propagator matrix G of the vector autoregressive model for the Fourier coefficients in vector form.
get.propagator.vec(wave, indCos, zeta, rho1, gamma, alpha, muX, muY, dt = 1,ns=4)get.propagator.vec(wave, indCos, zeta, rho1, gamma, alpha, muX, muY, dt = 1,ns=4)
wave |
Spatial wavenumbers. |
indCos |
Vector of integers indicating the position of columns in 'wave' of wavenumbers of cosine terms. |
zeta |
Damping parameter |
rho1 |
Range parameter of the diffusion term |
gamma |
Parameter that determines the amount of anisotropy in the diffusion term |
alpha |
Parameter that determines the direction of anisotropy in the diffusion term |
muX |
X component of the drift vector. |
muY |
Y component of the drift vector. |
dt |
Temporal lag between two time points. By default, this equals 1. |
ns |
Number of real Fourier functions that have only a cosine and no sine term. 'ns' is maximal 4. |
A list with three elements 'G11C', 'G11', and 'G12'. The first element contains a vector of length 'ns' which corresponds to the diagonal propagator of the cosin-only terms. The second element contains the remaining diagonal entries of G, i.e., the diagonal entries of the propagator for the cosine / sine pairs. Note that for each pair, only one value is taken since the diagonal elements for both the cosin and sine terms are equal. The third element is a vector with the off-diagonal terms of the propagator for the cosine / sine pairs.
Fabio Sigrist
##For illustration, four grid points on each axis n <- 4 wave <- wave.numbers(n) G <- get.propagator(wave=wave$wave,indCos=wave$indCos,zeta=0.5, rho1=0.1, gamma=2,alpha=pi/4, muX=0.2, muY=-0.15,dt=1,ns=4) diag(G)[1:4] diag(G[wave$indCos,wave$indCos]) diag(G[wave$indCos,wave$indCos+1]) get.propagator.vec(wave=wave$wave,indCos=wave$indCos,zeta=0.5, rho1=0.1, gamma=2,alpha=pi/4, muX=0.2, muY=-0.15,dt=1,ns=4)##For illustration, four grid points on each axis n <- 4 wave <- wave.numbers(n) G <- get.propagator(wave=wave$wave,indCos=wave$indCos,zeta=0.5, rho1=0.1, gamma=2,alpha=pi/4, muX=0.2, muY=-0.15,dt=1,ns=4) diag(G)[1:4] diag(G[wave$indCos,wave$indCos]) diag(G[wave$indCos,wave$indCos+1]) get.propagator.vec(wave=wave$wave,indCos=wave$indCos,zeta=0.5, rho1=0.1, gamma=2,alpha=pi/4, muX=0.2, muY=-0.15,dt=1,ns=4)
Returns the matrix that applies the two-dimensional real Fourier transform.
get.real.dft.mat(wave, indCos, ns = 4, n)get.real.dft.mat(wave, indCos, ns = 4, n)
wave |
Matrix of size 2 x NF with spatial wavenumbers. NF is the number of Fourier functions. |
indCos |
Vector of integers indicating the position of columns in 'wave' of wavenumbers of cosine terms. |
ns |
Number of real Fourier functions that have only a cosine and no sine term. 'ns' is maximal 4. |
n |
Number of grid points on each axis. n x n is the total number of spatial points. |
A matrix that applies the two-dimensional real Fourier transform.
Fabio Sigrist
##Example nr. 1: sampling from a Matern field n <- 50 spateFT <- spate.init(n=n,T=1) spec <- matern.spec(wave=spateFT$wave,n=n,rho0=0.05,sigma2=1,norm=TRUE) Phi <- get.real.dft.mat(wave=spateFT$wave, indCos=spateFT$indCos, n=n) sim <- Phi %*% (sqrt(spec)*rnorm(n*n)) image(1:n,1:n,matrix(sim,nrow=n),main="Sample from Matern field",xlab="",ylab="") ##Example nr. 2: image reconstruction n <- 50##Number of points on each axis ##Low-dimensional: only 41 Fourier functions spateFT <- spate.init(n=n,T=17,NF=45) Phi.LD <- get.real.dft.mat(wave=spateFT$wave, indCos=spateFT$indCos, ns=spateFT$ns, n=n) ##Mid-dimensional: 545 (of potentially 2500) Fourier functions spateFT <- spate.init(n=n,T=17,NF=101) Phi.MD <- get.real.dft.mat(wave=spateFT$wave, indCos=spateFT$indCos, ns=spateFT$ns, n=n) ##High-dimensional: all 2500 Fourier functions spateFT <- spate.init(n=n,T=17,NF=2500) Phi.HD <- get.real.dft.mat(wave=spateFT$wave, indCos=spateFT$indCos, ns=spateFT$ns, n=n) ##Define image image <- rep(0,n*n) for(i in 1:n){ for(j in 1:n){ image[(i-1)*n+j] <- cos(5*(i-n/2)/n*pi)*sin(5*(j)/n*pi)*(1-abs(i/n-1/2)-abs(j/n-1/2)) } } opar <- par(no.readonly = TRUE) par(mfrow=c(2,2),mar=c(2,3,2,1)) image(1:n, 1:n, matrix(image, nrow = n),col = cols(),xlab="",ylab="",main="Original image") ##Aply inverse Fourier transform, dimension reduction, and Fourier transform spec.LD <- t(Phi.LD) %*% image image.LD <- Phi.LD %*% spec.LD spec.MD <- t(Phi.MD) %*% image image.MD <- Phi.MD %*% spec.MD spec.HD <- t(Phi.HD) %*% image image.HD <- Phi.HD %*% spec.HD image(1:n, 1:n, matrix(image.LD, nrow = n),col = cols(), xlab="",ylab="",main="45 of 2500 Fourier terms") image(1:n, 1:n, matrix(image.MD, nrow = n),col = cols(), xlab="",ylab="",main="101 of 2500 Fourier terms") image(1:n, 1:n, matrix(image.HD, nrow = n),col = cols(), xlab="",ylab="",main="All 2500 Fourier terms") par(opar) # Reset par() settings##Example nr. 1: sampling from a Matern field n <- 50 spateFT <- spate.init(n=n,T=1) spec <- matern.spec(wave=spateFT$wave,n=n,rho0=0.05,sigma2=1,norm=TRUE) Phi <- get.real.dft.mat(wave=spateFT$wave, indCos=spateFT$indCos, n=n) sim <- Phi %*% (sqrt(spec)*rnorm(n*n)) image(1:n,1:n,matrix(sim,nrow=n),main="Sample from Matern field",xlab="",ylab="") ##Example nr. 2: image reconstruction n <- 50##Number of points on each axis ##Low-dimensional: only 41 Fourier functions spateFT <- spate.init(n=n,T=17,NF=45) Phi.LD <- get.real.dft.mat(wave=spateFT$wave, indCos=spateFT$indCos, ns=spateFT$ns, n=n) ##Mid-dimensional: 545 (of potentially 2500) Fourier functions spateFT <- spate.init(n=n,T=17,NF=101) Phi.MD <- get.real.dft.mat(wave=spateFT$wave, indCos=spateFT$indCos, ns=spateFT$ns, n=n) ##High-dimensional: all 2500 Fourier functions spateFT <- spate.init(n=n,T=17,NF=2500) Phi.HD <- get.real.dft.mat(wave=spateFT$wave, indCos=spateFT$indCos, ns=spateFT$ns, n=n) ##Define image image <- rep(0,n*n) for(i in 1:n){ for(j in 1:n){ image[(i-1)*n+j] <- cos(5*(i-n/2)/n*pi)*sin(5*(j)/n*pi)*(1-abs(i/n-1/2)-abs(j/n-1/2)) } } opar <- par(no.readonly = TRUE) par(mfrow=c(2,2),mar=c(2,3,2,1)) image(1:n, 1:n, matrix(image, nrow = n),col = cols(),xlab="",ylab="",main="Original image") ##Aply inverse Fourier transform, dimension reduction, and Fourier transform spec.LD <- t(Phi.LD) %*% image image.LD <- Phi.LD %*% spec.LD spec.MD <- t(Phi.MD) %*% image image.MD <- Phi.MD %*% spec.MD spec.HD <- t(Phi.HD) %*% image image.HD <- Phi.HD %*% spec.HD image(1:n, 1:n, matrix(image.LD, nrow = n),col = cols(), xlab="",ylab="",main="45 of 2500 Fourier terms") image(1:n, 1:n, matrix(image.MD, nrow = n),col = cols(), xlab="",ylab="",main="101 of 2500 Fourier terms") image(1:n, 1:n, matrix(image.HD, nrow = n),col = cols(), xlab="",ylab="",main="All 2500 Fourier terms") par(opar) # Reset par() settings
Auxilary function for the conversion between the complex FFT and the real Fourier transform.
index.complex.to.real.dft(n)index.complex.to.real.dft(n)
n |
Number of points on each axis. n x n is the total number of spatial points. |
A a list of indices used for the conversion between the complex FFT and the real Fourier transform.
Fabio Sigrist
Spectrum of the innovation term epsilon.
innov.spec(wave,n,ns=4,rho0,sigma2,zeta,rho1,alpha,gamma,nu=1,dt=1,norm=TRUE)innov.spec(wave,n,ns=4,rho0,sigma2,zeta,rho1,alpha,gamma,nu=1,dt=1,norm=TRUE)
wave |
Spatial wavenumbers. |
n |
Number of grid points on each axis. n x n is the total number of spatial points. |
ns |
Number of real Fourier functions that have only a cosine and no sine term. 'ns' is maximal 4. |
rho0 |
Range of the Matern covariance funtion for the innovation term epsilon |
sigma2 |
Marginal variance of the Matern covariance funtion for the innovation term epsilon |
zeta |
Damping parameter |
rho1 |
Range parameter of the diffusion term |
alpha |
Parameter that determines the direction of anisotropy in the diffusion term |
gamma |
Parameter that determines the amount of anisotropy in the diffusion term |
nu |
Smoothness parameter of the Matern covariance function for the innovations. By default, this equals 1 corresponding to the Whittle covariance function. |
dt |
Temporal lag between two time points. By default, this equals 1. |
norm |
logical; if 'TRUE' the spectrum is multiplied by n*n so that after applying the real Fourier transform 'real.FFT' one has the correct normalization. |
Vector with the spectrum of the integrated innovation term epsilon hat.
Fabio Sigrist
n <- 100 spec <- innov.spec(wave=spate.init(n=n,T=1)$wave,n=n,rho0=0.05,sigma2=0.5,zeta=0.5, rho1=0.05,alpha=pi/4,gamma=2,norm=TRUE) sim <- real.fft(sqrt(spec)*rnorm(n*n),n=n,inv=FALSE) image(1:n,1:n,matrix(sim,nrow=n),main="Sample from the integrated stochastic innovation",xlab="",ylab="",col=cols())n <- 100 spec <- innov.spec(wave=spate.init(n=n,T=1)$wave,n=n,rho0=0.05,sigma2=0.5,zeta=0.5, rho1=0.05,alpha=pi/4,gamma=2,norm=TRUE) sim <- real.fft(sqrt(spec)*rnorm(n*n),n=n,inv=FALSE) image(1:n,1:n,matrix(sim,nrow=n),main="Sample from the integrated stochastic innovation",xlab="",ylab="",col=cols())
Calculates the linear predictor.
lin.pred(x, beta)lin.pred(x, beta)
x |
Covariates in an array of dimensions p x T X N, where p denotes the number of covariates, T the number of time points, and N the number of spatial points. |
beta |
Coefficients of covariates in a vector of length p. |
Matrix of dimension T x N with linear predictors.
Fabio Sigrist
Evaluates the log-likelihood of the hyperparameters given the data (Gaussian case) or given the latent variable w (in the Tobit case).
loglike(par=NULL,w=NULL,wFT=NULL,x=NULL,spec=NULL,Gvec=NULL,tau2=NULL,n,T, NF=n*n,indCos=(1:((n*n-4)/2)*2+3),ns=4,nu=1,dt=1,logScale=FALSE, logInd=c(1,2,3,4,5,9),negative=FALSE)loglike(par=NULL,w=NULL,wFT=NULL,x=NULL,spec=NULL,Gvec=NULL,tau2=NULL,n,T, NF=n*n,indCos=(1:((n*n-4)/2)*2+3),ns=4,nu=1,dt=1,logScale=FALSE, logInd=c(1,2,3,4,5,9),negative=FALSE)
par |
Vector of parameters for the SPDE in the following order: rho_0, sigma^2, zeta, rho_1, gamma, alpha, mu_x, mu_y, tau^2, regression coefficients beta. rho_0 and sigma^2 are the range and marginal variance of the Whittle covariance funtion for the innovation term epsilon. zeta is the damping parameter. rho_1, gamma, and alpha parametrize the diffusion matrix with rho_1 being a range parameter, gamma and alpha determining the amount and the direction, respectively, of anisotropy. mu_x and mu_y are the two components of the drift vector. tau^2 denotes the variance of nugget effect or measurment error. Subsequently in par are the regression coefficients beta, if there are covariates. |
w |
Matrix of size T x N, where T and N denote the number of points in time and space. In the case of a Gaussian data model, w contains the observed values, with the Tobit model, w denotes the latent normal variable. |
wFT |
A vector with the (discrete) Fourier transform of the observed or latent w, depending on which data model is used. Note that, in contrast to w, this needs to be in stacked vector format. Use 'TSmat.to.vect'. |
x |
Covariates in an array of dimensions p x T X N, where p denotes the number of covariates, T the number of time points, and n the number of spatial points. |
spec |
A vector containing the spectrum of the innovation term epsilon. If 'spec' is not given, it is constructed based on 'par'. |
Gvec |
The propagator matrix G in vector format obtained from 'get.G.vec'. If 'Gvec' is not given, it is constructed based on 'par'. |
tau2 |
Measurement error variance tau2. If 'NULL'; tau2=par[9]. |
n |
Number of grid points on each axis. n x n is the total number of spatial points. |
T |
Number of points in time. |
NF |
Number of Fourier functions. |
indCos |
Vector of integers indicating the position of wavenumbers of cosine-only terms. |
ns |
Number of real Fourier functions that have only a cosine and no sine term. 'ns' is maximal 4. |
nu |
Smoothness parameter of the Matern covariance function for the innovations. By default this equals 1 corresponding to the Whittle covariance function. |
dt |
Temporal lag between two time points. By default, this equals 1. |
logScale |
logical; if 'TRUE' the parameters specified in 'logInd' are on the logarithmic scale. This is used for constraining parameters to be positive. |
logInd |
Vector of integers indicating which parameters are on the log-scale. |
negative |
logical; if 'TRUE' the negative log-likelihood is returned otherwise the positive log-likelihood is returned. |
Value of the log-likelihood evaluated at 'par'.
Fabio Sigrist
n <- 20 T <- 20 ##Specify hyper-parameters par <- c(rho0=0.1,sigma2=0.2,zeta=0.5,rho1=0.1,gamma=2,alpha=pi/4,muX=0.2,muY=-0.2,tau2=0.01) ##Simulate data spateSim <- spate.sim(par=par,n=n,T=T,seed=4) w <- spateSim$w ##Initial values for optim. This takes a couple of seconds. parI <- c(rho0=0.2,sigma2=0.1,zeta=0.25,rho1=0.01,gamma=1,alpha=0.3,muX=0,muY=0,tau2=0.005) logInd=c(1,2,3,4,5,9) ##Transform to log-scale parI[logInd] <- log(parI[logInd]) ##Fourier transform needs to be done only once wFT <- real.fft.TS(w,n=n,T=T) ##ML estimation using optim, takes a couple of seconds ##Load the precomputed object a line below to save time ##spateMLE <- optim(par=parI,loglike,control=list(trace=TRUE,maxit=1000),wFT=wFT,method="L-BFGS-B", ## lower=c(-10,-10,-10,-10,-10,0,-0.5,-0.5,-10), ## upper=c(10,10,10,10,10,pi/2,0.5,0.5,10),negative=TRUE, ## logScale=TRUE,hessian=TRUE,n=n,T=T) data("spateMLE") mle <- spateMLE$par mle[logInd] <- exp(mle[logInd]) sd=sqrt(diag(solve(spateMLE$hessian))) MleConfInt <- data.frame(array(0,c(4,9))) colnames(MleConfInt) <- names(par) rownames(MleConfInt) <- c("True","Estimate","Lower","Upper") MleConfInt[1,] <- par MleConfInt[2,] <- mle MleConfInt[3,] <- spateMLE$par-2*sd MleConfInt[4,] <- spateMLE$par+2*sd MleConfInt[c(3,4),logInd] <- exp(MleConfInt[c(3,4),logInd]) cat("\n") round(MleConfInt,digits=4)n <- 20 T <- 20 ##Specify hyper-parameters par <- c(rho0=0.1,sigma2=0.2,zeta=0.5,rho1=0.1,gamma=2,alpha=pi/4,muX=0.2,muY=-0.2,tau2=0.01) ##Simulate data spateSim <- spate.sim(par=par,n=n,T=T,seed=4) w <- spateSim$w ##Initial values for optim. This takes a couple of seconds. parI <- c(rho0=0.2,sigma2=0.1,zeta=0.25,rho1=0.01,gamma=1,alpha=0.3,muX=0,muY=0,tau2=0.005) logInd=c(1,2,3,4,5,9) ##Transform to log-scale parI[logInd] <- log(parI[logInd]) ##Fourier transform needs to be done only once wFT <- real.fft.TS(w,n=n,T=T) ##ML estimation using optim, takes a couple of seconds ##Load the precomputed object a line below to save time ##spateMLE <- optim(par=parI,loglike,control=list(trace=TRUE,maxit=1000),wFT=wFT,method="L-BFGS-B", ## lower=c(-10,-10,-10,-10,-10,0,-0.5,-0.5,-10), ## upper=c(10,10,10,10,10,pi/2,0.5,0.5,10),negative=TRUE, ## logScale=TRUE,hessian=TRUE,n=n,T=T) data("spateMLE") mle <- spateMLE$par mle[logInd] <- exp(mle[logInd]) sd=sqrt(diag(solve(spateMLE$hessian))) MleConfInt <- data.frame(array(0,c(4,9))) colnames(MleConfInt) <- names(par) rownames(MleConfInt) <- c("True","Estimate","Lower","Upper") MleConfInt[1,] <- par MleConfInt[2,] <- mle MleConfInt[3,] <- spateMLE$par-2*sd MleConfInt[4,] <- spateMLE$par+2*sd MleConfInt[c(3,4),logInd] <- exp(MleConfInt[c(3,4),logInd]) cat("\n") round(MleConfInt,digits=4)
Maps non-gridded data to a grid based on the coordinates supplied. Cells with no data are NA. For cells with more than one data point, the average is taken.
map.obs.to.grid(n,y.non.grid,coord,lengthx=NULL,lengthy=NULL)map.obs.to.grid(n,y.non.grid,coord,lengthx=NULL,lengthy=NULL)
y.non.grid |
Observed data in an T x N matrix with columns and rows corresponding to time and space, respectively. The coordinates of each observation point need to be specified in 'coord'. |
coord |
Matrix of dimension N x 2 with coordinates of the N observation points. Based on to these coordinates, each observation location is then mapped to a grid cell. |
lengthx |
Use together with 'coord' to specify the length of the x-axis. This is usefull if the observations lie in a rectangular area instead of a square. The length needs to be at least as large as the largest x-distance in 'coord. |
lengthy |
Use together with 'coord' to specify the length of the y-axis. This is usefull if the observations lie in a rectangular area instead of a square. The length needs to be at least as large as the largest y-distance in 'coord. |
n |
Number of point per axis of the square into which the points are mapped. In total, the process is modeled on a grid of size n*n. |
The function returns data in an T x n^2 matrix with columns and rows corresponding to time and space, respectively. Cells with no data are NA. For cells with more than one data point, the average is taken.
Fabio Sigrist
## See code of 'spate.mcmc'.## See code of 'spate.mcmc'.
Spectrum of the Matern covariance function. Note that the spectrum is
renormalized, by dividing with the sum over all frequencies so that
they sum to one, so that
is the marginal variance no matter how many
wavenumbers are included.
matern.spec(wave, n, ns=4, rho0, sigma2, nu = 1, norm = TRUE)matern.spec(wave, n, ns=4, rho0, sigma2, nu = 1, norm = TRUE)
wave |
Spatial wavenumbers. |
n |
Number of grid points on each axis. n x n is the total number of spatial points. |
ns |
Integer indicating the number of cosine-only terms. Maximally this is 4. |
rho0 |
Range parameter. |
sigma2 |
Marginal variance parameter. |
nu |
Smoothness parameter of the Matern covariance function. By default this equals 1 corresponding to the Whittle covariance function. |
norm |
logical; if 'TRUE' the spectrum is multiplied by n*n so that after applying the real Fourier transform 'real.FFT' one has the correct normalization. |
The Matern covariance function is of the form
with 'd' being the Euclidean distance between two points and K_nu(.) a modified Bessel function. Its spectrum is given by
where 'w' is a spatial wavenumber.
Vector with the spectrum of the Matern covariance function.
Fabio Sigrist
n <- 100 spec <- matern.spec(wave=spate.init(n=n,T=1)$wave,n=n,rho0=0.05,sigma2=1,norm=TRUE) sim <- real.fft(sqrt(spec)*rnorm(n*n),n=n,inv=FALSE) image(1:n,1:n,matrix(sim,nrow=n),main="Sample from a Gaussian process with Matern covariance function",xlab="",ylab="",col=cols())n <- 100 spec <- matern.spec(wave=spate.init(n=n,T=1)$wave,n=n,rho0=0.05,sigma2=1,norm=TRUE) sim <- real.fft(sqrt(spec)*rnorm(n*n),n=n,inv=FALSE) image(1:n,1:n,matrix(sim,nrow=n),main="Sample from a Gaussian process with Matern covariance function",xlab="",ylab="",col=cols())
Auxilary function for summarizing MCMC output and illustrating the posterior distributions.
mcmc.summary(data, probs = c(0.025, 0.5, 0.975), mean = FALSE)mcmc.summary(data, probs = c(0.025, 0.5, 0.975), mean = FALSE)
data |
Matrix of size p x Nmc where p denotes the number of parameters and Nmc the number of MCMC samples. |
probs |
Vector of quantiles that should be computed for each parameter. |
mean |
logical; if 'TRUE' the mean of the posterior distributions is computed as well. |
Matrix with quantiles and the mean of the posterior distributions.
Fabio Sigrist
data("spateMCMC") mcmc.summary(spateMCMC$Post,mean=TRUE)data("spateMCMC") mcmc.summary(spateMCMC$Post,mean=TRUE)
Default prior for direction of anisotropy in diffusion parameter alpha. A uniform prior on [0,pi/4] is used.
Palpha(alpha, log = FALSE)Palpha(alpha, log = FALSE)
alpha |
A quantile |
log |
Indicates whether the logarithm should be calculated or not |
Value of (log) prior at 'alpha'.
Fabio Sigrist
Default prior for amount of anisotropy in diffusion parameter gamma. A uniform prior on log(gamma) over the interval [1/100,100] is used.
Pgamma(gamma, log = FALSE)Pgamma(gamma, log = FALSE)
gamma |
A quantile |
log |
Indicates whether the logarithm should be calculated or not |
Value of (log) prior at 'gamma'.
Fabio Sigrist
Default prior for transformation parameter of the Tobit model. A locally constant, improper prior on the positive real line is used.
Plambda(lambda, log = FALSE)Plambda(lambda, log = FALSE)
lambda |
A quantile |
log |
Indicates whether the logarithm should be calculated or not |
Value of (log) prior at 'lambda'.
Fabio Sigrist
Plots trace plots, pair plots, the posterior of the hyperparameters and the posterior of the latent spatio-temporal process.
## S3 method for class 'spateMCMC' plot(x,..., trace = TRUE, hist = TRUE, medianHist=TRUE, pairs = FALSE,ask = TRUE, ToFile = FALSE, path = NULL,file = NULL,true=NULL,BurnInAdaptive=NULL, postProcess = FALSE)## S3 method for class 'spateMCMC' plot(x,..., trace = TRUE, hist = TRUE, medianHist=TRUE, pairs = FALSE,ask = TRUE, ToFile = FALSE, path = NULL,file = NULL,true=NULL,BurnInAdaptive=NULL, postProcess = FALSE)
x |
A 'spateMCMC' object obtained from 'spate.mcmc'. |
... |
Arguments to be passed to 'spate.plot' in case 'postProcess=TRUE' is selected. |
trace |
logical; if 'TRUE' trace plots are made |
hist |
logical; if 'TRUE' histograms of the posterior distributions for the hyper-parameters are plotted |
medianHist |
logical; if 'TRUE' medians are added to the histograms. |
pairs |
logical; if 'TRUE' trace plots are made |
ask |
logical; if 'TRUE' (and the R session is interactive) the user is asked for input, before a new figure is drawn. |
ToFile |
logical; if 'TRUE' the plots are save to a file. |
path |
The path. |
file |
The file name. |
true |
The true value of the parameters (for simulation only). |
BurnInAdaptive |
The number of samples used as burn-in before starting the adaptive estimation of Metropolis-Hastings proposal covariance for the hyper-parameters. |
postProcess |
logical; if 'TRUE' the posterior of the spatio-temporal process xi is plotted as well. |
Plots illustrating a fitted model saved in a 'spateMCMC' object.
Fabio Sigrist
data("spateMCMC") plot(spateMCMC,medianHist=TRUE,pairs=TRUE)data("spateMCMC") plot(spateMCMC,medianHist=TRUE,pairs=TRUE)
This is the plotting function for 'spateSim' objects. It calles the function 'spate.plot()'.
## S3 method for class 'spateSim' plot(x,..., plotXi =TRUE,plotW = FALSE)## S3 method for class 'spateSim' plot(x,..., plotXi =TRUE,plotW = FALSE)
x |
'spateSim' object obtained from 'spate.sim'. |
... |
Arguments to be passed to 'spate.plot' |
plotXi |
Logical; if 'TRUE' the latent process 'xi' is plotted. |
plotW |
Logical; if 'TRUE' the observed process 'w' is plotted. |
Plots illustrating the simulated space-time field.
Fabio Sigrist
spateSim <-spate.sim(par=c(rho0=0.1,sigma2=0.2,zeta=0.5,rho1=0.1,gamma=2, alpha=pi/4,muX=0.2,muY=-0.2,tau2=0.01),n=50,T=9) plot(spateSim)spateSim <-spate.sim(par=c(rho0=0.1,sigma2=0.2,zeta=0.5,rho1=0.1,gamma=2, alpha=pi/4,muX=0.2,muY=-0.2,tau2=0.01),n=50,T=9) plot(spateSim)
Default prior for x-component of drift vector mu. A uniform prior on the interval [-0.5,0.5] is used.
Pmux(mux, log = FALSE)Pmux(mux, log = FALSE)
mux |
A quantile |
log |
Indicates whether the logarithm should be calculated or not |
Value of (log) prior at 'mux'.
Fabio Sigrist
Default prior for y-component of drift vector mu. A uniform prior on the interval [-0.5,0.5] is used.
Pmuy(muy, log = FALSE)Pmuy(muy, log = FALSE)
muy |
A quantile |
log |
Indicates whether the logarithm should be calculated or not |
Value of (log) prior at 'muy'.
Fabio Sigrist
Histogram of posterior distributions.
post.dist.hist(data, true=NULL, breaks = 20, mean = FALSE, median = TRUE)post.dist.hist(data, true=NULL, breaks = 20, mean = FALSE, median = TRUE)
data |
Matrix of size p x Nmc where p denotes the number of parameters and Nmc the number of MCMC samples. |
true |
The true value of the parameters (for simulation only). |
breaks |
Parameter for 'hist()' function. |
mean |
logical; if 'TRUE' the mean is added to the histogram. |
median |
logical; if 'TRUE' the median is added to the histogram. |
Histograms illustrating posterior distributions.
Fabio Sigrist
Default prior for range parameter rho0 of stochastic source-sink term epsilon. A uniform prior on [0,100] is used.
Prho0(rho0, log = FALSE)Prho0(rho0, log = FALSE)
rho0 |
A quantile |
log |
Indicates whether the logarithm should be calculated or not |
Value of (log) prior at 'rho0'.
Fabio Sigrist
Default prior for range parameter rho1 of diffusive term. A uniform prior on [0,100] is used.
Prho1(rho1, log = FALSE)Prho1(rho1, log = FALSE)
rho1 |
A quantile. |
log |
Indicates whether the logarithm should be calculated or not |
Value of (log) prior at 'rho1'.
Fabio Sigrist
Print function for spateMCMC objects.
## S3 method for class 'spateMCMC' print(x,...)## S3 method for class 'spateMCMC' print(x,...)
x |
A 'spateMCMC' object obtained from 'spate.mcmc'. |
... |
not used. |
Fabio Sigrist
Print function for 'spateSim' objects.
## S3 method for class 'spateSim' print(x,...)## S3 method for class 'spateSim' print(x,...)
x |
'spateSim' object obtained from 'spate.sim'. |
... |
Arguments to be passed to 'spate.plot' |
Fabio Sigrist
spateSim <-spate.sim(par=c(rho0=0.1,sigma2=0.2,zeta=0.5,rho1=0.1,gamma=2, alpha=pi/4,muX=0.2,muY=-0.2,tau2=0.01),n=50,T=9) spateSimspateSim <-spate.sim(par=c(rho0=0.1,sigma2=0.2,zeta=0.5,rho1=0.1,gamma=2, alpha=pi/4,muX=0.2,muY=-0.2,tau2=0.01),n=50,T=9) spateSim
Function that propagates the vector 'alphat'. This is equivalent to multiplying 'alphat' with the propagator matrix G. It is a lot faster though, due to the block-diagonal structure of G. This is a wrapper function of a C function.
propagate.spectral(alphat,spateFT=NULL,n=NULL,Gvec=NULL,par=NULL)propagate.spectral(alphat,spateFT=NULL,n=NULL,Gvec=NULL,par=NULL)
alphat |
A vector of spectral coefficients. |
spateFT |
A 'spateFT' obtained from 'spate.init'. Either this or 'n' needs to be given. |
n |
Number of points on each axis. n x n is the total number of spatial points. Either this or 'spateFT' needs to be given. |
Gvec |
The propagator matrix G in vector format obtained from 'get.propagator.vec'. If 'Gvec' is not given, it is constructed based on 'par'. |
par |
Parameters for the SPDE in the following order: rho_0, sigma^2, zeta, rho_1, gamma, alpha, mu_x, mu_y, tau^2. If 'Gvec' is not given, 'par' needs to be given. |
A vector of propagated coefficients G*alphat.
Fabio Sigrist
n <- 50 spec <- matern.spec(wave=spate.init(n=n,T=1)$wave,n=n,rho0=0.05,sigma2=1,norm=TRUE) alphat <- sqrt(spec)*rnorm(n*n) ##Propagate initial state wave <- wave.numbers(n) Gvec <- get.propagator.vec(wave=wave$wave,indCos=wave$indCos,zeta=0.1,rho1=0.02,gamma=2, alpha=pi/4,muX=0.2,muY=0.2,dt=1,ns=4) alphat1 <- propagate.spectral(alphat,n=n,Gvec=Gvec) opar <- par(no.readonly = TRUE) par(mfrow=c(1,2)) image(1:n,1:n,matrix(real.fft(alphat,n=n,inv=FALSE),nrow=n),main="Whittle field",xlab="",ylab="",col=cols()) image(1:n,1:n,matrix(real.fft(alphat1,n=n,inv=FALSE),nrow=n),main="Propagated field",xlab="",ylab="",col=cols()) par(opar) # Reset par() settingsn <- 50 spec <- matern.spec(wave=spate.init(n=n,T=1)$wave,n=n,rho0=0.05,sigma2=1,norm=TRUE) alphat <- sqrt(spec)*rnorm(n*n) ##Propagate initial state wave <- wave.numbers(n) Gvec <- get.propagator.vec(wave=wave$wave,indCos=wave$indCos,zeta=0.1,rho1=0.02,gamma=2, alpha=pi/4,muX=0.2,muY=0.2,dt=1,ns=4) alphat1 <- propagate.spectral(alphat,n=n,Gvec=Gvec) opar <- par(no.readonly = TRUE) par(mfrow=c(1,2)) image(1:n,1:n,matrix(real.fft(alphat,n=n,inv=FALSE),nrow=n),main="Whittle field",xlab="",ylab="",col=cols()) image(1:n,1:n,matrix(real.fft(alphat1,n=n,inv=FALSE),nrow=n),main="Propagated field",xlab="",ylab="",col=cols()) par(opar) # Reset par() settings
Default prior for marginal variance parameter sigma2 (=sigma^2) of the stochastic source-sink term epsilon. A uniform, improper prior on sigma (P[sigma] propto 1 or P[sigma2] propto 1/tau) is used.
Psigma2(sigma2, log = FALSE)Psigma2(sigma2, log = FALSE)
sigma2 |
A quantile |
log |
Indicates whether the logarithm should be calculated or not |
Value of (log) prior at 'sigma2'.
Fabio Sigrist
Default prior for measurment error or small scale variation tau2 (nugget effect). A uniform, improper prior on tau (P[tau] propto 1 or P[tau2] propto 1/tau) is used.
Ptau2(tau2, log = FALSE)Ptau2(tau2, log = FALSE)
tau2 |
A quantile |
log |
Indicates whether the logarithm should be calculated or not |
Value of (log) prior at tau2.
Fabio Sigrist
Default prior for damping parameter zeta. A uniform, improper prior on the positive real line (P[zeta] propto 1) is used.
Pzeta(zeta, log = FALSE)Pzeta(zeta, log = FALSE)
zeta |
A quantile |
log |
Indicates whether the logarithm should be calculated or not |
Value of (log) prior at 'zeta'.
Fabio Sigrist
Fast calculation of the real Fourier transform. This is a wrapper function for a C function which uses the complex FFT function from the 'fftw3' library.
real.fft(w,n,inv=TRUE,indFFT=NULL)real.fft(w,n,inv=TRUE,indFFT=NULL)
w |
A spatial field in a stacked vector of length N=n^2. |
n |
Number of grid points on each axis. n x n is the total number of spatial points. |
inv |
Indicates whether the inverse Fourier transform should be calculated or not. |
indFFT |
A list of containing vectors of natural numbers representing indices used to transform between the real and the complex Fourier transform. |
A vector of length n*n containing the real (inverse) Fourier transformation of 'w'.
Fabio Sigrist
n <- 100 spec <- matern.spec(wave=spate.init(n=n,T=1)$wave,n=n,rho0=0.05,sigma2=1,norm=TRUE) sim <- real.fft(sqrt(spec)*rnorm(n*n),n=n,inv=FALSE) image(1:n,1:n,matrix(sim,nrow=n),main="Sample from Matern field",xlab="",ylab="")n <- 100 spec <- matern.spec(wave=spate.init(n=n,T=1)$wave,n=n,rho0=0.05,sigma2=1,norm=TRUE) sim <- real.fft(sqrt(spec)*rnorm(n*n),n=n,inv=FALSE) image(1:n,1:n,matrix(sim,nrow=n),main="Sample from Matern field",xlab="",ylab="")
This function calculates the two-dimensional real Fourier transform of a space-time field. This is a wrapper function for a C function which uses the complex FFT function from the 'fftw3' library. In contrast to using T times the function 'real.FFT', R needs to communicate with C only once and not T times which saves computational time.
real.fft.TS(w,n,T,inv=TRUE,indFFT=NULL)real.fft.TS(w,n,T,inv=TRUE,indFFT=NULL)
w |
Spatio-temporal field in a stacked vector of length T x N. Stacking is done first over space and then time. E.g., the first N=n^2 entries contain the spatial field at time t=1. Note that the spatial field itself is stacked as well, i.e., each spatial field is in a vector of length N=n^2. |
n |
Number of grid points on each axis. n x n is the total number of spatial points. |
T |
Number of time points. |
inv |
Indicates whether the inverse Fourier transform should be calculated or not. |
indFFT |
A list of containing vectors of natural numbers representing indices used to transform between the real and the complex Fourier transform. |
A vector of length T x N containing the (inverse) Fourier transform of 'w'.
Fabio Sigrist
n <- 100 T <- 4 spec <- matern.spec(wave=spate.init(n=n,T=1)$wave,n=n,rho0=0.05,sigma2=1,norm=TRUE) specsim <- matrix(0,nrow=T,ncol=n*n) for(t in 1:T) specsim[t,] <- rnorm(n*n)*sqrt(spec) maternsim <- vect.to.TSmat(real.fft.TS(TSmat.to.vect(specsim),n=n,T=T,inv=FALSE),T=T) opar <- par(no.readonly = TRUE) par(mfrow=c(2,2)) for(t in 1:T) image(1:n,1:n,matrix(maternsim[t,],nrow=n), main="Sample from Matern field",xlab="",ylab="") par(opar) # Reset par() settingsn <- 100 T <- 4 spec <- matern.spec(wave=spate.init(n=n,T=1)$wave,n=n,rho0=0.05,sigma2=1,norm=TRUE) specsim <- matrix(0,nrow=T,ncol=n*n) for(t in 1:T) specsim[t,] <- rnorm(n*n)*sqrt(spec) maternsim <- vect.to.TSmat(real.fft.TS(TSmat.to.vect(specsim),n=n,T=T,inv=FALSE),T=T) opar <- par(no.readonly = TRUE) par(mfrow=c(2,2)) for(t in 1:T) image(1:n,1:n,matrix(maternsim[t,],nrow=n), main="Sample from Matern field",xlab="",ylab="") par(opar) # Reset par() settings
Sample from the full conditional of the Fourier coefficients.
sample.four.coef(w=NULL,wFT=NULL,spec=NULL,Gvec=NULL,tau2=NULL,par=NULL,n,T, NF=n*n,indCos=(1:((n*n-4)/2)*2+3),ns=4,nu=1,dt=1)sample.four.coef(w=NULL,wFT=NULL,spec=NULL,Gvec=NULL,tau2=NULL,par=NULL,n,T, NF=n*n,indCos=(1:((n*n-4)/2)*2+3),ns=4,nu=1,dt=1)
w |
Observed data or latent process w (depending on which data model is used) in an T x n*n matrix with columns and rows (points on a grid stacked into a vector) corresponding to time and space, respectively. |
wFT |
Vector of length T*n*n containing the real Fourier transform of 'w'. |
spec |
Spectrum of the innovations |
Gvec |
The propagator matrix G in vector format obtained from 'get.G.vec'. If 'Gvec' is not given, it is constructed based on 'par'. |
tau2 |
Measurement error variance tau2. If 'NULL'; tau2=par[9]. |
par |
Vector of parameters for the SPDE in the following order: rho_0, sigma^2, zeta, rho_1, gamma, alpha, mu_x, mu_y, tau^2. If 'spec' and 'Gvec' are given, 'par' will not be used. |
n |
Number of grid points on each axis. n*n is the total number of spatial points. |
T |
Number of points in time. |
NF |
Number of Fourier functions used. |
indCos |
Vector of integers indicating the position cosine terms in the 1:NF real Fourier functions. The first 'ns' cosine wavenumbers in 'wave' are not included in 'indCos'. |
ns |
Number of real Fourier functions that have only a cosine and no sine term. 'ns' is maximal 4. |
nu |
Smoothness parameter of the Matern covariance function for the innovations. By default this equals 1 corresponding to the Whittle covariance function. |
dt |
Temporal lag between two time points. By default, this equals 1. |
A T x n*n matrix with a sample from the full conditional
of latent process .
Fabio Sigrist
##Specifications for simulated example n <- 50 T <- 4 par <- c(rho0=0.1,sigma2=0.2,zeta=0.5,rho1=0.1,gamma=2,alpha=pi/4,muX=0.2,muY=-0.2,tau2=0.01) spateSim <- spate.sim(par=par,n=n,T=T,seed=4) w <- spateSim$w ##Sample from full conditional Nmc <- 50 alphaS <- array(0,c(T,n*n,Nmc)) wFT <- real.fft.TS(w,n=n,T=T) for(i in 1:Nmc){ alphaS[,,i] <- sample.four.coef(wFT=wFT,par=par,n=n,T=T,NF=n*n) } ##Mean from full conditional alphaMean <- apply(alphaS,c(1,2),mean) xiMean <- real.fft.TS(alphaMean,n=n,T=T,inv=FALSE) opar <- par(no.readonly = TRUE) par(mfrow=c(2,4),mar=c(1,1,1,1)) for(t in 1:4) image(1:n,1:n,matrix(w[t,],nrow=n),xlab="",ylab="",col=cols(), main=paste("w(",t,")",sep=""),xaxt='n',yaxt='n') for(t in 1:4) image(1:n,1:n,matrix(xiMean[t,],nrow=n),xlab="",ylab="",col=cols(), main=paste("xiPost(",t,")",sep=""),xaxt='n',yaxt='n') par(opar) # Reset par() settings##Specifications for simulated example n <- 50 T <- 4 par <- c(rho0=0.1,sigma2=0.2,zeta=0.5,rho1=0.1,gamma=2,alpha=pi/4,muX=0.2,muY=-0.2,tau2=0.01) spateSim <- spate.sim(par=par,n=n,T=T,seed=4) w <- spateSim$w ##Sample from full conditional Nmc <- 50 alphaS <- array(0,c(T,n*n,Nmc)) wFT <- real.fft.TS(w,n=n,T=T) for(i in 1:Nmc){ alphaS[,,i] <- sample.four.coef(wFT=wFT,par=par,n=n,T=T,NF=n*n) } ##Mean from full conditional alphaMean <- apply(alphaS,c(1,2),mean) xiMean <- real.fft.TS(alphaMean,n=n,T=T,inv=FALSE) opar <- par(no.readonly = TRUE) par(mfrow=c(2,4),mar=c(1,1,1,1)) for(t in 1:4) image(1:n,1:n,matrix(w[t,],nrow=n),xlab="",ylab="",col=cols(), main=paste("w(",t,")",sep=""),xaxt='n',yaxt='n') for(t in 1:4) image(1:n,1:n,matrix(xiMean[t,],nrow=n),xlab="",ylab="",col=cols(), main=paste("xiPost(",t,")",sep=""),xaxt='n',yaxt='n') par(opar) # Reset par() settings
Auxilary function for constructing 'spateFT' objects which are used for the two-dimensional Fourier transform.
spate.init(n,T,NF=n*n)spate.init(n,T,NF=n*n)
n |
Number of points on each axis. n x n is the total number of spatial points. |
T |
Number of temporal points. |
NF |
This integer specifies the number of Fourier functions. If NF<n*n, dimension reduction is obtained. In this case, Fourier functions with wavenumbers closest to the origin (0,0) are first included. If a given 'NF' implies a basis with anisotropy, i.e., higher frequencies in one direction than in another, this is automatically corrected by using the next higher integer NF' such that the basis has the same resolution in all directions. |
A 'spateFT' object. This is a list with
wave |
a matrix containing the wavenumbers |
indCos |
a vector indicating the position of the cosine terms (excluding the frist 'ns') |
ns |
an integer indicating the number of cosine-only terms |
indFFT |
a list of indices used for the conversion between the complex FFT and the real Fourier transform. |
n |
number of points on each axis |
T |
number of points in time |
Fabio Sigrist
MCMC algorithm for fitting the model.
spate.mcmc(y,coord=NULL,lengthx=NULL,lengthy=NULL,Sind=NULL,n=NULL, IncidenceMat=FALSE,x=NULL,SV=c(rho0=0.2,sigma2=0.1, zeta=0.25,rho1=0.2,gamma=1,alpha=0.3,muX=0,muY=0,tau2=0.005), betaSV=rep(0,dim(x)[1]),RWCov=NULL,parh=NULL,tPred=NULL, sPred=NULL,P.rho0=Prho0,P.sigma2=Psigma2,P.zeta=Pzeta,P.rho1=Prho1, P.gamma=Pgamma,P.alpha=Palpha,P.mux=Pmux,P.muy=Pmuy,P.tau2=Ptau2, lambdaSV=1,sdlambda=0.01,P.lambda=Plambda,DataModel="Normal", DimRed=FALSE,NFour=NULL,indEst=1:9,Nmc=10000,BurnIn =1000, path=NULL,file=NULL,SaveToFile=FALSE,PlotToFile=FALSE, FixEffMetrop=TRUE,saveProcess=FALSE,Nsave=200,seed=NULL, Padding=FALSE,adaptive=TRUE,NCovEst=500,BurnInCovEst=500, MultCov=0.5,printRWCov=FALSE,MultStdDevLambda=0.75, Separable=FALSE,Drift=!Separable,Diffusion=!Separable, logInd=c(1,2,3,4,5,9),nu=1,plotTrace=TRUE, plotHist=FALSE,plotPairs=FALSE,trueVal=NULL, plotObsLocations=FALSE,trace=TRUE,monitorProcess=FALSE, tProcess=NULL,sProcess=NULL)spate.mcmc(y,coord=NULL,lengthx=NULL,lengthy=NULL,Sind=NULL,n=NULL, IncidenceMat=FALSE,x=NULL,SV=c(rho0=0.2,sigma2=0.1, zeta=0.25,rho1=0.2,gamma=1,alpha=0.3,muX=0,muY=0,tau2=0.005), betaSV=rep(0,dim(x)[1]),RWCov=NULL,parh=NULL,tPred=NULL, sPred=NULL,P.rho0=Prho0,P.sigma2=Psigma2,P.zeta=Pzeta,P.rho1=Prho1, P.gamma=Pgamma,P.alpha=Palpha,P.mux=Pmux,P.muy=Pmuy,P.tau2=Ptau2, lambdaSV=1,sdlambda=0.01,P.lambda=Plambda,DataModel="Normal", DimRed=FALSE,NFour=NULL,indEst=1:9,Nmc=10000,BurnIn =1000, path=NULL,file=NULL,SaveToFile=FALSE,PlotToFile=FALSE, FixEffMetrop=TRUE,saveProcess=FALSE,Nsave=200,seed=NULL, Padding=FALSE,adaptive=TRUE,NCovEst=500,BurnInCovEst=500, MultCov=0.5,printRWCov=FALSE,MultStdDevLambda=0.75, Separable=FALSE,Drift=!Separable,Diffusion=!Separable, logInd=c(1,2,3,4,5,9),nu=1,plotTrace=TRUE, plotHist=FALSE,plotPairs=FALSE,trueVal=NULL, plotObsLocations=FALSE,trace=TRUE,monitorProcess=FALSE, tProcess=NULL,sProcess=NULL)
y |
Observed data in an T x N matrix with columns and rows corresponding to time and space (observations on a grid stacked into a vector), respectively. By default, at each time point, the observations are assumed to lie on a square grid with each axis scaled so that it has unit length. |
coord |
If specified, this needs to be a matrix of dimension N x 2 with coordinates of the N observation points. Observations in 'y' can either be on a square grid or not. If not, the coordinates of each observation point need to be specified in 'coord'. According to these coordinates, each observation location is then mapped to a grid cell. If 'coord' is not specified, the observations in 'y' are assumed to lie on a square grid with each axis scaled so that it has unit length. |
lengthx |
Use together with 'coord' to specify the length of the x-axis. This is usefull if the observations lie in a rectangular area instead of a square. The length needs to be at least as large as the largest x-distance in 'coord. |
lengthy |
Use together with 'coord' to specify the length of the y-axis. This is usefull if the observations lie in a rectangular area instead of a square. The length needs to be at least as large as the largest y-distance in 'coord. |
Sind |
Vector of indices of grid cells where observations are made, in case, the observation are not made at every grid cell. Alternatively, the coordinates of the observation locations can be specfied in 'coord'. |
n |
Number of point per axis of the square into which the points are mapped. In total, the process is modeled on a grid of size n*n. |
IncidenceMat |
Logical; if 'TRUE' an incidence matrix relating the latent process to observation locations is used. This is only recommended to use when the observations are relatively low-dimensional and when the latent process is modeled in a reduced dimensional space as well. |
x |
Covariates in an array of dimensions p x T X N, where p denotes the number of covariates, T the number of time points, and N the number of spatial points. |
SV |
Starting values for parameters. Parameters for the SPDE in the following order: rho_0, sigma^2, zeta, rho_1, gamma, alpha, mu_x, mu_y, tau^2. rho_0 and sigma^2 are the range and marginal variance of the Matern covariance funtion for the innovation term epsilon. zeta is the damping parameter. rho_1, gamma, and alpha parametrize the diffusion matrix with rho_1 being a range parameter, gamma and alpha determining the amount and the direction, respectively, of anisotropy. mu_x and mu_y are the two components of the drift vector. tau^2 denotes the nugget effect or measurment error. |
betaSV |
Starting values for regression coefficients. |
RWCov |
Covariance matrix of the proposal distribution used in the random walk Metropolis-Hastings step for the hyper-parameters. |
parh |
Only used in prediction mode. If 'parh' is not 'NULL', this indicates that 'spate.mcmc' is used for making predictions at locations (tPred,sPred) instead of applying the traditional MCMC algorithm. In case 'parh' is not 'NULL', it is a Npar x Nsim matrix containing Nsim samples from the posterior of the Npar parameters. This argument is used by the wrapper function 'spate.predict'. |
tPred |
Time points where predictions are made.This needs to be a vector if predictions are made at multiple times. For instance, if T is the number of time points in the data 'y', then tPred=c(T+1, T+2) means that predictions are made at time 'T+1' and 'T+2'. This argument is used by the wrapper function 'spate.predict'. |
sPred |
Vector of indices of grid cells (positions of locations in the stacked spatial vector) where predictions are made. This argument is used by the wrapper function 'spate.predict'. |
P.rho0 |
Function specifying the prior for rho0. |
P.sigma2 |
Function specifying the prior for sigma2. |
P.zeta |
Function specifying the prior for zeta. |
P.rho1 |
Function specifying the prior for rho1. |
P.gamma |
Function specifying the prior for gamma. |
P.alpha |
Function specifying the prior for alpha. |
P.mux |
Function specifying the prior for mux. |
P.muy |
Function specifying the prior for muy. |
P.tau2 |
Function specifying the prior for tau2. |
lambdaSV |
Starting value for transformation parameter lambda in the Tobit model. |
sdlambda |
Standard deviation of the proposal distribution used in the random walk Metropolis-Hastings step for lambda. |
P.lambda |
Function specifying the prior for lambda. |
DataModel |
Specifies the data model. "Normal" or "SkewTobit" are available options. |
DimRed |
Logical; if 'TRUE' dimension reduction is applied. This means that not the full number (n*n) of Fourier functions is used but rather only a reduced dimensional basis of dimension 'NFour'. |
NFour |
If 'DimRed' is 'TRUE', this specifies the number of Fourier functions. |
indEst |
A vector of numbers specifying which for which parameters the posterior should be computed and which should be held fix (at their starting value). If the corresponding to the index of rho_0, sigma^2, zeta, rho_1, gamma, alpha, mu_x, mu_y, tau^2 is present in the vector, the parameter will be estimated otherwise not. Default is indEst=1:9 which means that one samples from the posterior for all parameters. |
Nmc |
Number of MCMC samples. |
BurnIn |
Length of the burn-in period. |
path |
Path, in case plots and / or the spateMCMC object should be save in a file. |
file |
File name, in case plots and / or the spateMCMC object should be save in a file. |
SaveToFile |
Indicates whether the spateMCMC object should be save in a file. |
PlotToFile |
Indicates whether the MCMC output analysis plots should be save in a file. |
FixEffMetrop |
The fixed effects, i.e., the regression coefficients, can either be sampled in a Gibbs step or updated together with the hyperparameters in the Metropolis-Hastings step. The latter is the default and recommended option since correlations between fixed effects and the random process can result in slow mixing. |
saveProcess |
Logical; if 'TRUE' samples from the posterior of the latent spatio-temporal process xi are saved. |
Nsave |
Number of samples from the posterior of the latent spatio-temporal process xi that should be save. |
seed |
Seed for random generator. |
Padding |
Indicates whether padding is applied or not. If the range parameters are large relative to the domain, this is recommended since otherwise spurious periodicity can occur. |
adaptive |
Indicates whether an adaptive Metropolis-Hastings algorithm is used or not. If yes, the proposal covariance matrix 'RWCov' is adaptively estimated during the algorithm and tuning does not need to be done by hand. |
NCovEst |
Minimal number of samples to be used for estimating the proposal matrix. |
BurnInCovEst |
Burn-in period for estimating the proposal matrix. |
MultCov |
Numeric used as multiplier for the adaptively estimated proposal cocariance matrix 'RWCov' of the hyper-parameters. I.e., the estimated covariance matrix is multiplied by 'MultCov'. |
printRWCov |
Logical, if 'TRUE' the estimated proposal cocariance matrix is printed each time. |
MultStdDevLambda |
Numeric used as multiplier for the adaptively estimated proposal standard deviation of the Tobit transformation parameter lambda. I.e., the estimated standard deviation is multiplied by 'MultStdDevLambda'. |
Separable |
Indicates whether a separable model, i.e., no transport / drift and no diffusion, should be estimated. |
Drift |
Indicates whether a drift term should be included. |
Diffusion |
Indicates whether a diffusion term should be included. |
logInd |
Indicates which parameters are sampled on the log-scale. Default is logInd=c(1, 2, 3, 4, 5, 9) corresponding to rho_0, sigma2, zeta, rho_1, gamma, and tau^2. |
nu |
Smoothness parameter of the Matern covariance function for the innovations. By default this equals 1 corresponding to the Whittle covariance function. |
plotTrace |
Indicates whether trace plots are made. |
plotHist |
Indicates whether histograms of the posterior distributions are made. |
plotPairs |
Indicates whether scatter plots of the hyper-parameters and the regression coefficients are made. |
trueVal |
In simulations, true values can be supplied for comparison with the MCMC output. |
plotObsLocations |
Logical; if 'TRUE' the observations locations are ploted together with the grid cells. |
trace |
Logical; if 'TRUE' tracing information on the progress of the MCMC algorithm is produced. |
monitorProcess |
Logical; if 'TRUE' in addition to the trace plots of the hyper-parameters, the mixing properties of the latent process xi=Phi*alpha is monitored. This is done by plotting the current sample of the process. More specifically, the time series at locations 'sProcess' and the spatial fieldd at time points 'tProcess'. |
tProcess |
To be secified if 'monitorProcess=TRUE'. Time points at which spatial fields of the sampled process should be plotted. |
sProcess |
To be secified if 'monitorProcess=TRUE'. Locations at which time series of the sampled process should be plotted. |
The function returns a 'spateMCMC' object with, amongst others, the following entries
Post |
Matrix containing samples from the posterior of the hyper-parameters and the regression coefficient |
xiPost |
Array with samples from the posterior of the spatio-temporal process |
RWCov |
(Estimated) proposal covariance matrix |
Fabio Sigrist
##Specify hyper-parameters par <- c(rho0=0.1,sigma2=0.2,zeta=0.5,rho1=0.1,gamma=2,alpha=pi/4,muX=0.2,muY=-0.2,tau2=0.01) ##Simulate data spateSim <- spate.sim(par=par,n=20,T=20,seed=4) w <- spateSim$w ##Below is an example to illustrate the use of the MCMC algorithm. ##In practice, more samples are needed for a sufficiently large effective sample size. ##The following takes a couple of minutes. ##Load the precomputed object some lines below to save time. ##spateMCMC <- spate.mcmc(y=w,x=NULL,SV=c(rho0=0.2,sigma2=0.1, ## zeta=0.25,rho1=0.2,gamma=1,alpha=0.3,muX=0,muY=0,tau2=0.005), ## RWCov=diag(c(0.005,0.005,0.05,0.005,0.005,0.001,0.0002,0.0002,0.0002)), ## Nmc=10000,BurnIn=2000,seed=4,Padding=FALSE,plotTrace=TRUE,NCovEst=500, ## BurnInCovEst=500,trueVal=par,saveProcess=TRUE) ##spateMCMC ##plot(spateMCMC.fit,true=par,postProcess=TRUE) ##Instead of waiting, you can also use this precomputed object data("spateMCMC") spateMCMC plot(spateMCMC,true=par,medianHist=FALSE)##Specify hyper-parameters par <- c(rho0=0.1,sigma2=0.2,zeta=0.5,rho1=0.1,gamma=2,alpha=pi/4,muX=0.2,muY=-0.2,tau2=0.01) ##Simulate data spateSim <- spate.sim(par=par,n=20,T=20,seed=4) w <- spateSim$w ##Below is an example to illustrate the use of the MCMC algorithm. ##In practice, more samples are needed for a sufficiently large effective sample size. ##The following takes a couple of minutes. ##Load the precomputed object some lines below to save time. ##spateMCMC <- spate.mcmc(y=w,x=NULL,SV=c(rho0=0.2,sigma2=0.1, ## zeta=0.25,rho1=0.2,gamma=1,alpha=0.3,muX=0,muY=0,tau2=0.005), ## RWCov=diag(c(0.005,0.005,0.05,0.005,0.005,0.001,0.0002,0.0002,0.0002)), ## Nmc=10000,BurnIn=2000,seed=4,Padding=FALSE,plotTrace=TRUE,NCovEst=500, ## BurnInCovEst=500,trueVal=par,saveProcess=TRUE) ##spateMCMC ##plot(spateMCMC.fit,true=par,postProcess=TRUE) ##Instead of waiting, you can also use this precomputed object data("spateMCMC") spateMCMC plot(spateMCMC,true=par,medianHist=FALSE)
Generates a figure or an animation of a spatio-temporal field.
spate.plot(xi,nx=NULL,whichT=NULL,format="ImgTogether",ToFile=FALSE,path=NULL, file=NULL,indScale=FALSE,main=NULL,mfrow=NULL, imagesize=c(1000,1000),zlim=NULL,breaks=NULL,...)spate.plot(xi,nx=NULL,whichT=NULL,format="ImgTogether",ToFile=FALSE,path=NULL, file=NULL,indScale=FALSE,main=NULL,mfrow=NULL, imagesize=c(1000,1000),zlim=NULL,breaks=NULL,...)
xi |
A spatio-temporal field stored in an T x N matrix with columns and rows corresponding to time and space, respectively. |
nx |
Integer specifying the number of points on the x-axis. If 'NULL', a quadratic grid is assumed. |
whichT |
Vector of integers specifying the time points that are plotted. If 'NULL', all time points are plotted. |
format |
A string specifying how the spatio-temporal field should be ploted. "ImgTogether" produces one single plot containing all spatial fields at all time points. With "ImgSeparate", the spatial fields at each time point are plotted in separate plots. |
ToFile |
Indicates whether the output should be saved to a file. |
path |
Path indicating where to save the file. |
file |
File name. |
indScale |
Indicates whether the color scale for the spatial plots is the same for all time points (indScale=FALSE) or separate for each time point (indScale=TRUE). |
main |
Titles for the plots. Can be either be NULL or a character vector of length equal to the number of time points or 1. |
mfrow |
See 'par'. Can be either NULL or an integer vector of length two. If it is NULL, the functions determines mfrow automatically. |
imagesize |
The size of the .jpeg image if ToFile=TRUE. |
zlim |
Graphical parameter to be passed to 'image'. Determines the scale on the z-axis of the plots. If 'indScale=FALSE' one can specify the common scale on the z-axis of the plots through this argument. |
breaks |
Graphical parameter to be passed to 'image'. |
... |
Other graphical parameters that are passed to 'image' and 'par'. |
Plots illustrating a space-time field.
Fabio Sigrist
spateSim <- spate.sim(par=c(rho0=0.1,sigma2=0.2,zeta=0.5,rho1=0.1,gamma=2, alpha=pi/4,muX=0.2,muY=-0.2,tau2=0.01),n=50,T=9) spate.plot(spateSim$xi)spateSim <- spate.sim(par=c(rho0=0.1,sigma2=0.2,zeta=0.5,rho1=0.1,gamma=2, alpha=pi/4,muX=0.2,muY=-0.2,tau2=0.01),n=50,T=9) spate.plot(spateSim$xi)
Obtain samples from predictive distribution in space and time given the posterior of the hyperparameters.
spate.predict(y,tPred,sPred=NULL,xPred=NULL,yPred=NULL,spateMCMC,Nsim=200, BurnIn=5,coord=NULL,lengthx=NULL,lengthy=NULL,Sind=NULL, n=NULL,IncidenceMat=FALSE,x=NULL,DataModel="Normal", DimRed=FALSE,NFour=NULL,seed=NULL,nu =1,trace=FALSE)spate.predict(y,tPred,sPred=NULL,xPred=NULL,yPred=NULL,spateMCMC,Nsim=200, BurnIn=5,coord=NULL,lengthx=NULL,lengthy=NULL,Sind=NULL, n=NULL,IncidenceMat=FALSE,x=NULL,DataModel="Normal", DimRed=FALSE,NFour=NULL,seed=NULL,nu =1,trace=FALSE)
y |
Observed data in an T x N matrix with columns and rows corresponding to time and space, respectively. |
x |
Covariates in an array of dimensions p x T X N, where p denotes the number of covariates, T the number of time points, and N the number of spatial points. |
tPred |
Time points where predictions are made.This needs to be a vector if predictions are made at multiple times. For instance, if T is the number of time points in the data 'y', then tPred=c(T+1, T+2) means that predictions are made at time 'T+1' and 'T+2'. If 'xPred' and 'yPred' are empty, then predictions are made at all spatial points for each time point in 'tPred'. Otherwise 'xPred' and 'yPred', or 'sPred, need to have the same length as 'tPred', and predictions are made at the points (tPred,xPred,yPred), ore (tPred, sPred), respectively. |
sPred |
Vector of indices of grid cells (positions of locations in the stacked spatial vector) where predictions are made. This is an alternative to specifying the coordinates 'xPred' and 'yPred'. |
xPred |
Vector of x-coordinates of spatial points where predictions are made. This is an alternative to specifying the grid cell in 'sPred'. |
yPred |
Vector of y-coordinates of spatial points where predictions are made. This is an alternative to specifying the grid cell in 'sPred'. |
spateMCMC |
'spateMCMC' object obtained from 'spate.mcmc' containing the posterior of the hyper-parameters and information on the model used. |
Nsim |
Number of samples used to characterize the predictive distribution. |
BurnIn |
Length of burn-in period. |
coord |
If specified, this needs to be a matrix of dimension N x 2 with coordinates of the N observation points. Observations in 'y' can either be on a square grid or not. If not, the coordinates of each observation point need to be specified in 'coord'. According to these coordinates, each observation location is then mapped to a grid cell. If 'coord' is not specified, the observations in 'y' are assumed to lie on a square grid with each axis scaled so that it has unit length. |
lengthx |
Use together with 'coord' to specify the length of the x-axis. This is usefull if the observations lie in a rectangular area instead of a square. The length needs to be at least as large as the largest x-distance in 'coord. |
lengthy |
Use together with 'coord' to specify the length of the y-axis. This is usefull if the observations lie in a rectangular area instead of a square. The length needs to be at least as large as the largest y-distance in 'coord. |
Sind |
Vector of indices of grid cells where observations are made, in case, the observation are not made at every grid cell. Alternatively, the coordinates of the observation locations can be specfied in 'coord'. |
n |
Number of point per axis of the square into which the points are mapped. In total, the process is modeled on a grid of size n*n. |
IncidenceMat |
Logical; if 'TRUE' and incidence matrix relating the latent process to observation locations. This is only recommended to use when the observations are relatively low-dimensional and when the latent process is modeled in a reduced dimensional spaceas well. |
DataModel |
Specifies the data model. "Normal" or "SkewTobit". |
DimRed |
Logical; if 'TRUE' dimension reduction is applied. This means that not the full number (n*n) of Fourier functions is used but rather only a reduced dimensional basis of dimension 'NFour'. |
NFour |
If 'DimRed' is 'TRUE', this specifies the number of Fourier functions. |
seed |
Seed for random generator. |
nu |
Smoothness parameter of the Matern covariance function for the innovations. By default this equals 1 corresponding to the Whittle covariance function. |
trace |
Logical; if 'TRUE' tracing information on the progress of the MCMC algorithm is produced. |
Depending on whether 'xPred' and 'yPred' are empty or not, either
(i) |
an array of size t x s x Nsim where the first index is for time, the second for space, and the third for the number of samples 'Nsim' |
or
(ii) |
a matrix of size length(tPred) x Nsim |
Fabio Sigrist
par <- c(rho0=0.1,sigma2=0.2,zeta=0.5,rho1=0.1,gamma=2,alpha=pi/4,muX=0.2,muY=-0.2,tau2=0.01) ##Simulate data spateSim <- spate.sim(par=par,n=20,T=20,seed=4) w <- spateSim$w data("spateMCMC") ##Make predictions. Takes a couple of seconds predict <- spate.predict(y=w,tPred=(17:25),spateMCMC=spateMCMC,Nsim =200, BurnIn=10,DataModel="Normal") Pmean <- apply(predict,c(1,2),mean) Psd <- apply(predict,c(1,2),sd) opar <- par(no.readonly = TRUE) par(mfrow=c(2,2)) image(1:20,1:20,matrix(w[19,],nrow=20),main="Observed field at t=19",xlab="x",ylab="y") image(1:20,1:20,matrix(Pmean[3,],nrow=20),main="Fitted field at t=19",xlab="x",ylab="y") image(1:20,1:20,matrix(w[20,],nrow=20),main="Observed field at t=20",xlab="x",ylab="y") image(1:20,1:20,matrix(Pmean[4,],nrow=20),main="Fitted field at t=20",xlab="x",ylab="y") par(mfrow=c(3,3)) zlim=c(min(Pmean),max(Pmean)) for(i in 1:9){ image(1:20,1:20,matrix(Pmean[i,],nrow=20),zlim=zlim, main=paste("Mean t=",i+16,sep=""),xlab="x",ylab="y") } par(mfrow=c(3,3)) zlim=c(min(Psd),max(Psd)) for(i in 1:9){ image(1:20,1:20,matrix(Psd[i,],nrow=20),zlim=zlim, main=paste("Std.dev. t=",i+16,sep=""),xlab="x",ylab="y") par(opar) # Reset par() settings }par <- c(rho0=0.1,sigma2=0.2,zeta=0.5,rho1=0.1,gamma=2,alpha=pi/4,muX=0.2,muY=-0.2,tau2=0.01) ##Simulate data spateSim <- spate.sim(par=par,n=20,T=20,seed=4) w <- spateSim$w data("spateMCMC") ##Make predictions. Takes a couple of seconds predict <- spate.predict(y=w,tPred=(17:25),spateMCMC=spateMCMC,Nsim =200, BurnIn=10,DataModel="Normal") Pmean <- apply(predict,c(1,2),mean) Psd <- apply(predict,c(1,2),sd) opar <- par(no.readonly = TRUE) par(mfrow=c(2,2)) image(1:20,1:20,matrix(w[19,],nrow=20),main="Observed field at t=19",xlab="x",ylab="y") image(1:20,1:20,matrix(Pmean[3,],nrow=20),main="Fitted field at t=19",xlab="x",ylab="y") image(1:20,1:20,matrix(w[20,],nrow=20),main="Observed field at t=20",xlab="x",ylab="y") image(1:20,1:20,matrix(Pmean[4,],nrow=20),main="Fitted field at t=20",xlab="x",ylab="y") par(mfrow=c(3,3)) zlim=c(min(Pmean),max(Pmean)) for(i in 1:9){ image(1:20,1:20,matrix(Pmean[i,],nrow=20),zlim=zlim, main=paste("Mean t=",i+16,sep=""),xlab="x",ylab="y") } par(mfrow=c(3,3)) zlim=c(min(Psd),max(Psd)) for(i in 1:9){ image(1:20,1:20,matrix(Psd[i,],nrow=20),zlim=zlim, main=paste("Std.dev. t=",i+16,sep=""),xlab="x",ylab="y") par(opar) # Reset par() settings }
Generates one sample from the Gaussian process specified through the SPDE.
spate.sim(par,n,T,seed=NULL,StartVal=NULL,nu=1)spate.sim(par,n,T,seed=NULL,StartVal=NULL,nu=1)
par |
Vector of parameters for the SPDE in the following order: rho_0, sigma^2, zeta, rho_1, gamma, alpha, mu_x, mu_y, tau^2. rho_0 and sigma^2 are the range and marginal variance of the Matern covariance funtion for the innovation term epsilon. zeta is the damping parameter. rho_1, gamma, and alpha parametrize the diffusion matrix with rho_1 being a range parameter, gamma and alpha determining the amount and the direction, respectively, of anisotropy. mu_x and mu_y are the two components of the drift vector. tau^2 denotes the variance of nugget effect or measurment error. |
n |
Number of grid points on each axis. n x n is the total number of spatial points. |
T |
Number of points in time. |
seed |
Seed for random number generator. |
StartVal |
A starting value (field) for the SPDE can be defined. This is the spatial field at the initial time that get propagated forward by the SPDE. The starting fields needs to be a stacked vector of lengths n x n (number of spatial points). Use 'as.vector()' to convert a spatial matrix to a vector. |
nu |
Smoothness parameter of the Matern covariance function for the innovations. By default this equals 1 corresponding to the Whittle covariance function. |
A list containing a simulated spatio-temporal field xi with covariance structure as defined through the SPDE, a simulated observation field w obtained by adding a measurement error, and the simulated Fourier coefficients. The last two are returned only on demand.
Fabio Sigrist
StartVal <- rep(0,100^2) StartVal[75*100+75] <- 1000 par <- c(rho0=0.05,sigma2=0.7^2,zeta=-log(0.99),rho1=0.06, gamma=3,alpha=pi/4,muX=-0.1,muY=-0.1,tau2=0.00001) spateSim <- spate.sim(par=par,n=100,T=6,StartVal=StartVal,seed=1) plot(spateSim,mfrow=c(2,3),mar=c(2,2,2,2),indScale=TRUE, cex.axis=1.5,cex.main=2)StartVal <- rep(0,100^2) StartVal[75*100+75] <- 1000 par <- c(rho0=0.05,sigma2=0.7^2,zeta=-log(0.99),rho1=0.06, gamma=3,alpha=pi/4,muX=-0.1,muY=-0.1,tau2=0.00001) spateSim <- spate.sim(par=par,n=100,T=6,StartVal=StartVal,seed=1) plot(spateSim,mfrow=c(2,3),mar=c(2,2,2,2),indScale=TRUE, cex.axis=1.5,cex.main=2)
Precalculated 'spateMCMC' object containing a fitted model (MCMC output) obtained from 'spate.mcmc'.
spateMCMCspateMCMC
Precalculated maximum likelihood estimate using 'optim' and the function 'loglike'.
spateMLEspateMLE
Summary function for 'spateSim' objects.
## S3 method for class 'spateSim' summary(object,...)## S3 method for class 'spateSim' summary(object,...)
object |
'spateSim' object obtained from 'spate.sim()'. |
... |
not used. |
Fabio Sigrist
spateSim <- spate.sim(par=c(rho0=0.1,sigma2=0.2,zeta=0.5,rho1=0.1,gamma=2, alpha=pi/4,muX=0.2,muY=-0.2,tau2=0.01),n=50,T=9) summary(spateSim)spateSim <- spate.sim(par=c(rho0=0.1,sigma2=0.2,zeta=0.5,rho1=0.1,gamma=2, alpha=pi/4,muX=0.2,muY=-0.2,tau2=0.01),n=50,T=9) summary(spateSim)
Full conditional for transformation parameter lambda of the Tobit model. This is used in the Metropolis-Hastings step of the MCMC algorithm.
tobit.lambda.log.full.cond(y, z, tau2, lambda)tobit.lambda.log.full.cond(y, z, tau2, lambda)
y |
Observed data. |
z |
Latent Gaussian variable. |
tau2 |
Value of variance (corresponds to nugget effect). |
lambda |
Value of transformation parameter lambda. |
Value of evaluated full conditional for transformation parameter lambda.
Fabio Sigrist
Trace plots for MCMC output analysis.
trace.plot(data, true = NULL, BurnIn = NULL,BurnInAdaptive=NULL)trace.plot(data, true = NULL, BurnIn = NULL,BurnInAdaptive=NULL)
data |
A p x Nmc data.frame of matrix where p denotes the number of parameters and Nmc the number of Monte Carlo samples. |
true |
The true value of the parameters (for simulation only). |
BurnIn |
The number of samples used as burn-in if the burn-in has not yet been removed from the sample. |
BurnInAdaptive |
The number of samples used as burn-in before starting the adaptive estimation of Metropolis-Hastings proposal covariance for the hyper-parameters. |
Trace plots.
Fabio Sigrist
data <- matrix(rnorm(1200),nrow=6) opar <- par(no.readonly = TRUE) par(mfrow=c(2,3)) trace.plot(data,true=rep(0,6)) par(opar) # Reset par() settingsdata <- matrix(rnorm(1200),nrow=6) opar <- par(no.readonly = TRUE) par(mfrow=c(2,3)) trace.plot(data,true=rep(0,6)) par(opar) # Reset par() settings
Converts a time-space matrix with columns and rows corresponding to time and space into a stacked N*T vector.
TSmat.to.vect(mat)TSmat.to.vect(mat)
mat |
A T x N matrix with columns and rows corresponding to time and space, respectively. |
A vector of stacked values. Stacking is done first over space and then time.
Fabio Sigrist
vect <- 1:12 mat <- vect.to.TSmat(vect,T=3)##Convert vector to matrix TSmat.to.vect(mat)##Convert matrix to vector.vect <- 1:12 mat <- vect.to.TSmat(vect,T=3)##Convert vector to matrix TSmat.to.vect(mat)##Convert matrix to vector.
Converts a stacked N*T vector into a time-space matrix with columns and rows corresponding to time and space, respectively.
vect.to.TSmat(vect, T = 1)vect.to.TSmat(vect, T = 1)
vect |
A vector of stacked values. Stacking is done first over space and then time. |
T |
Number of time points. |
A T x N matrix with columns and rows corresponding to time and space, respectively.
Fabio Sigrist
vect <- 1:12 vect vect.to.TSmat(vect,T=3)vect <- 1:12 vect vect.to.TSmat(vect,T=3)
Calculates the Eucledian norm of a vector
vnorm(v)vnorm(v)
v |
Vector. |
The Eucledian norm of the vector 'v'.
Fabio Sigrist
v <- c(1,2) vnorm(v)v <- c(1,2) vnorm(v)
Returns wave numbers used in real Fourier transform.
wave.numbers(n)wave.numbers(n)
n |
Number of grid points on each axis. n x n is the total number of spatial points. |
Returns a list with
wave |
A 2 x n^2 matrix with wavenumbers used in the real Fourier transform. The first four columns contain the wavenumbers that are only used by cosine terms and not by sine terms. Subsequent columns alternate between wavenumbers of cosine and sine terms. |
indCos |
Vector of integers indicating the position of columns in 'wave' of wavenumbers of cosine terms. The first four cosine wavenumbers in 'wave' are not included in 'indCos'. |
Fabio Sigrist