| Title: | Simulate Auto Regressive Data from Precision Matrices |
|---|---|
| Description: | Using sparse precision matrices and Cholesky factorization simulates data that is auto-regressive. |
| Authors: | Neal Marquez [aut, cre, cph] |
| Maintainer: | Neal Marquez <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 0.1.0 |
| Built: | 2025-02-07 03:26:18 UTC |
| Source: | https://github.com/nmmarquez/ar.matrix |
Functions for creating precision matricies and observations of an AR1 process
Q.AR1(M, sigma, rho, sparse=TRUE, vcov=FALSE) r.AR1(n, M, sigma, rho)Q.AR1(M, sigma, rho, sparse=TRUE, vcov=FALSE) r.AR1(n, M, sigma, rho)
M |
int > 0, number of elements in the AR1 process. |
sigma |
float > 0, pairwise observation standard deviation. |
rho |
float >= 0 & < 1, how correlated pairwise observations are. The function will still run with values outside of the range [0,1) however the stability of the simulation results are not gaurunteed. |
sparse |
bool Should the matrix be of class 'dsCMatrix' |
vcov |
bool If the vcov matrix should be returned instead of the precision matrix. |
n |
int > 0, number of observations to simulate from the GMRF. |
Q.AR1 returns either a precision or variance-covariance function with a AR1 structure.
r.AR1 retrurns a matrix with n rows which are the n observations of a Gaussian Markov random field AR1 process.
require("ggplot2") # simulate AR1 GMRF obs <- r.AR1(100, M=30, sigma=1, rho=.98) # resulting matrix is n x M dim(obs) # subtract off the first time point to more easily observe correlation obs_adj <- obs - obs[,1] # move objects to a data frame ar1_df <- data.frame(obs=c(t(obs_adj)), realization=rep(1:100, each=30), time=rep(1:30, 100)) # plot each realization ggplot(data=ar1_df, aes(time, obs, group=realization, color=realization)) + geom_line()require("ggplot2") # simulate AR1 GMRF obs <- r.AR1(100, M=30, sigma=1, rho=.98) # resulting matrix is n x M dim(obs) # subtract off the first time point to more easily observe correlation obs_adj <- obs - obs[,1] # move objects to a data frame ar1_df <- data.frame(obs=c(t(obs_adj)), realization=rep(1:100, each=30), time=rep(1:30, 100)) # plot each realization ggplot(data=ar1_df, aes(time, obs, group=realization, color=realization)) + geom_line()
Functions for creating precision matricies and observations of a independent identically distributed GMRF process.
Q.iid(M, sigma, sparse=TRUE, vcov=FALSE) r.iid(n, M, sigma)Q.iid(M, sigma, sparse=TRUE, vcov=FALSE) r.iid(n, M, sigma)
M |
int > 0, number of elements in the process. |
sigma |
float > 0, standard deviat |
sparse |
bool Should the matrix be of class 'dsCMatrix' |
vcov |
bool If the vcov matrix should be returned instead of the precision matrix. |
n |
int > 0, number of observations to simulate from the GMRF. |
Q.iid returns either a precision or variance-covariance function with iid structure.
r.iid retrurns a matrix with n rows which are the n observations of a Gaussian Markov random field iid process.
require("leaflet") require("sp") # simulate iid data and attach to spatial polygons data frame US.df@data$data <- c(r.iid(1, M=nrow(US.graph), sigma=1)) # color palette of data pal <- colorNumeric(palette="YlGnBu", domain=US.df@data$data) # see map map1<-leaflet() %>% addProviderTiles("CartoDB.Positron") %>% addPolygons(data=US.df, fillColor=~pal(data), color="#b2aeae", fillOpacity=0.7, weight=0.3, smoothFactor=0.2) %>% addLegend("bottomright", pal=pal, values=US.df$data, title="", opacity=1) map1require("leaflet") require("sp") # simulate iid data and attach to spatial polygons data frame US.df@data$data <- c(r.iid(1, M=nrow(US.graph), sigma=1)) # color palette of data pal <- colorNumeric(palette="YlGnBu", domain=US.df@data$data) # see map map1<-leaflet() %>% addProviderTiles("CartoDB.Positron") %>% addPolygons(data=US.df, fillColor=~pal(data), color="#b2aeae", fillOpacity=0.7, weight=0.3, smoothFactor=0.2) %>% addLegend("bottomright", pal=pal, values=US.df$data, title="", opacity=1) map1
Functions for creating precision matricies and observations of a Leroux CAR(lCAR) process as defined in MacNab 2011. The matrix defines the precision of estimates when observations share connections which are conditionally auto-regressive(CAR).
Q.lCAR(graph, sigma, rho, sparse=TRUE, vcov=FALSE) r.lCAR(n, graph, sigma, rho)Q.lCAR(graph, sigma, rho, sparse=TRUE, vcov=FALSE) r.lCAR(n, graph, sigma, rho)
graph |
matrix, square matrix indicating where two observations are connected (and therefore conditionally auto-regressive). |
sigma |
float > 0, process standard derviation see MacNab 2011. |
rho |
float >= 0 & < 1, how correlated neighbors are. The function will still run with values outside of the range [0,1) however the stability of the simulation results are not gaurunteed. see MacNab 2011. |
sparse |
bool Should the matrix be of class 'dsCMatrix' |
vcov |
bool If the vcov matrix should be returned instead of the precision matrix. |
n |
int > 0, number of observations to simulate from the GMRF. |
Q.lCAR returns either a precision or variance-covariance function with a lCAR structure.
r.lCAR retrurns a matrix with n rows which are the n observations of a Gaussian Markov random field lCAR process.
Y.C. MacNab On Gaussian Markov random fields and Bayesian disease mapping. Statistical Methods in Medical Research. 2011.
require("leaflet") require("sp") # simulate lCAR data and attach to spatial polygons data frame US.df@data$data <- c(r.lCAR(1, graph=US.graph, sigma=1, rho=.99)) # color palette of data pal <- colorNumeric(palette="YlGnBu", domain=US.df@data$data) # see map map1<-leaflet() %>% addProviderTiles("CartoDB.Positron") %>% addPolygons(data=US.df, fillColor=~pal(data), color="#b2aeae", fillOpacity=0.7, weight=0.3, smoothFactor=0.2) %>% addLegend("bottomright", pal=pal, values=US.df$data, title="", opacity=1) map1require("leaflet") require("sp") # simulate lCAR data and attach to spatial polygons data frame US.df@data$data <- c(r.lCAR(1, graph=US.graph, sigma=1, rho=.99)) # color palette of data pal <- colorNumeric(palette="YlGnBu", domain=US.df@data$data) # see map map1<-leaflet() %>% addProviderTiles("CartoDB.Positron") %>% addPolygons(data=US.df, fillColor=~pal(data), color="#b2aeae", fillOpacity=0.7, weight=0.3, smoothFactor=0.2) %>% addLegend("bottomright", pal=pal, values=US.df$data, title="", opacity=1) map1
EXPIREMENTAL. Functions for creating precision matricies and observations of a modified BYM process as defined in MacNab 2011. The matrix defines the precision of estimates when observations share connections which are conditionally auto-regressive(CAR). Because the precision matrix is not symetric the process is not a true GMRF.
Q.mBYM(graph, sigma, rho, vcov=FALSE) r.mBYM(n, graph, sigma, rho)Q.mBYM(graph, sigma, rho, vcov=FALSE) r.mBYM(n, graph, sigma, rho)
graph |
matrix, square matrix indicating where two observations are connected (and therefore conditionally auto-regressive). |
sigma |
float > 0, process standard derviation see MacNab 2011. |
rho |
float >= 0 & < 1, how correlated neighbors are. The function will still run with values outside of the range [0,1) however the stability of the simulation results are not gaurunteed. see MacNab 2011. |
vcov |
bool If the vcov matrix should be returned instead of the precision matrix. |
n |
int > 0, number of observations to simulate from the GMRF. |
Q.mBYM returns either a precision or variance-covariance function with a modified BYM structure.
r.mBYM retrurns a matrix with n rows which are the n observations of a pseudo Gaussian Markov random field of a modified BYM process.
Y.C. MacNab On Gaussian Markov random fields and Bayesian disease mapping. Statistical Methods in Medical Research. 2011.
## Not run: require("leaflet") require("sp") # simulate mBYM data and attach to spatial polygons data frame US.df@data$data <- c(r.mBYM(1, graph=US.graph, sigma=1, rho=.99)) # color palette of data pal <- colorNumeric(palette="YlGnBu", domain=US.df@data$data) # see map map1<-leaflet() %>% addProviderTiles("CartoDB.Positron") %>% addPolygons(data=US.df, fillColor=~pal(data), color="#b2aeae", fillOpacity=0.7, weight=0.3, smoothFactor=0.2) %>% addLegend("bottomright", pal=pal, values=US.df$data, title="", opacity=1) map1 ## End(Not run)## Not run: require("leaflet") require("sp") # simulate mBYM data and attach to spatial polygons data frame US.df@data$data <- c(r.mBYM(1, graph=US.graph, sigma=1, rho=.99)) # color palette of data pal <- colorNumeric(palette="YlGnBu", domain=US.df@data$data) # see map map1<-leaflet() %>% addProviderTiles("CartoDB.Positron") %>% addPolygons(data=US.df, fillColor=~pal(data), color="#b2aeae", fillOpacity=0.7, weight=0.3, smoothFactor=0.2) %>% addLegend("bottomright", pal=pal, values=US.df$data, title="", opacity=1) map1 ## End(Not run)
Functions for creating precision matricies and observations of a proper CAR(pCAR) process as defined in MacNab 2011. The matrix defines the precision of estimates when observations share connections which are conditionally auto-regressive(CAR).
Q.pCAR(graph, sigma, rho, sparse=TRUE, vcov=FALSE) r.pCAR(n, graph, sigma, rho)Q.pCAR(graph, sigma, rho, sparse=TRUE, vcov=FALSE) r.pCAR(n, graph, sigma, rho)
graph |
matrix, square matrix indicating where two observations are connected (and therefore conditionally auto-regressive). |
sigma |
float > 0, process standard derviation see MacNab 2011. |
rho |
float >= 0 & < 1, how correlated neighbors are. The function will still run with values outside of the range [0,1) however the stability of the simulation results are not gaurunteed. see MacNab 2011. |
sparse |
bool Should the matrix be of class 'dsCMatrix' |
vcov |
bool If the vcov matrix should be returned instead of the precision matrix. |
n |
int > 0, number of observations to simulate from the GMRF. |
Q.pCAR returns either a precision or variance-covariance function with a pCAR structure.
r.pCAR retrurns a matrix with n rows which are the n observations of a Gaussian Markov random field pCAR process.
Y.C. MacNab On Gaussian Markov random fields and Bayesian disease mapping. Statistical Methods in Medical Research. 2011.
require("leaflet") require("sp") # simulate pCAR data and attach to spatial polygons data frame US.df@data$data <- c(r.pCAR(1, graph=US.graph, sigma=1, rho=.99)) # color palette of data pal <- colorNumeric(palette="YlGnBu", domain=US.df@data$data) # see map map1<-leaflet() %>% addProviderTiles("CartoDB.Positron") %>% addPolygons(data=US.df, fillColor=~pal(data), color="#b2aeae", fillOpacity=0.7, weight=0.3, smoothFactor=0.2) %>% addLegend("bottomright", pal=pal, values=US.df$data, title="", opacity=1) map1require("leaflet") require("sp") # simulate pCAR data and attach to spatial polygons data frame US.df@data$data <- c(r.pCAR(1, graph=US.graph, sigma=1, rho=.99)) # color palette of data pal <- colorNumeric(palette="YlGnBu", domain=US.df@data$data) # see map map1<-leaflet() %>% addProviderTiles("CartoDB.Positron") %>% addPolygons(data=US.df, fillColor=~pal(data), color="#b2aeae", fillOpacity=0.7, weight=0.3, smoothFactor=0.2) %>% addLegend("bottomright", pal=pal, values=US.df$data, title="", opacity=1) map1
Functions for creating precision matricies and observations of an RW1 process
Q.RW1(M, sigma, sparse=TRUE) r.RW1(n, M, sigma)Q.RW1(M, sigma, sparse=TRUE) r.RW1(n, M, sigma)
M |
int > 0, number of elements in the RW1 process. |
sigma |
float > 0, pairwise observation standard deviation. |
sparse |
bool Should the matrix be of class 'dsCMatrix' |
n |
int > 0, number of observations to simulate from the GMRF. |
Q.RW1 returns a precision matrix with a RW1 structure.
r.RW1 retrurns a matrix with n rows which are the n observations of an Intrinsic Gaussian Markov random field RW1 process.
require("ggplot2") # simulate RW1 GMRF obs <- r.RW1(100, M=30, sigma=1) # resulting matrix is n x M dim(obs) # subtract off the first time point to more easily observe correlation obs_adj <- obs - obs[,1] # move objects to a data frame rw1_df <- data.frame(obs=c(t(obs_adj)), realization=rep(1:100, each=30), time=rep(1:30, 100)) # plot each realization ggplot(data=rw1_df, aes(time, obs, group=realization, color=realization)) + geom_line()require("ggplot2") # simulate RW1 GMRF obs <- r.RW1(100, M=30, sigma=1) # resulting matrix is n x M dim(obs) # subtract off the first time point to more easily observe correlation obs_adj <- obs - obs[,1] # move objects to a data frame rw1_df <- data.frame(obs=c(t(obs_adj)), realization=rep(1:100, each=30), time=rep(1:30, 100)) # plot each realization ggplot(data=rw1_df, aes(time, obs, group=realization, color=realization)) + geom_line()
Functions for creating precision matricies and observations of an RW2 process
Q.RW2(M, sigma, sparse=TRUE) r.RW2(n, M, sigma)Q.RW2(M, sigma, sparse=TRUE) r.RW2(n, M, sigma)
M |
int > 0, number of elements in the RW2 process. |
sigma |
float > 0, pairwise observation standard deviation. |
sparse |
bool Should the matrix be of class 'dsCMatrix' |
n |
int > 0, number of observations to simulate from the GMRF. |
Q.RW2 returns a precision matrix with a RW2 structure.
r.RW2 retrurns a matrix with n rows which are the n observations of an Intrinsic Gaussian Markov random field RW2 process.
require("ggplot2") # simulate RW2 GMRF obs <- r.RW2(100, M=30, sigma=1) # resulting matrix is n x M dim(obs) # move objects to a data frame RW2_df <- data.frame(obs=c(t(obs)), realization=rep(1:100, each=30), time=rep(1:30, 100)) # plot each realization ggplot(data=RW2_df, aes(time, obs, group=realization, color=realization)) + geom_line()require("ggplot2") # simulate RW2 GMRF obs <- r.RW2(100, M=30, sigma=1) # resulting matrix is n x M dim(obs) # move objects to a data frame RW2_df <- data.frame(obs=c(t(obs)), realization=rep(1:100, each=30), time=rep(1:30, 100)) # plot each realization ggplot(data=RW2_df, aes(time, obs, group=realization, color=realization)) + geom_line()
Takes in a square precision matrix, which ideally should be sparse and using Choleski factorization simulates data from a mean 0 process where the inverse of the precision matrix represents the variance-covariance of the points in the process. The resulting simulants represent samples of a Gaussian Markov random field (GMRF).
sim.AR(n, Q)sim.AR(n, Q)
n |
int > 0, number of observations to simulate from the GMRF. |
Q |
matrix, a square precision matrix. |
Matrix object, matrix where each row is a single obsrevation from a GMRF with covariance structure Q^-1.
require("ggplot2") # simulate 2D ar1 process # pairwise correlation rho <- .95 # pairwise variance sigma <- .5 # 2 dimensions of simulations years <- 20 ages <- 10 # kronnecker product to get joint covariance Q2D <- kronecker(Q.AR1(M=years, sigma, rho), Q.AR1(M=ages, sigma, rho)) # simulate the data and place it in a data frame Q2D.df <- data.frame(obs=c(sim.AR(1, Q2D)), age=rep(1:ages, years), year=rep(1:years, each=ages)) # graph results ggplot(data=Q2D.df, aes(year, obs, group=age, color=age)) + geom_line()require("ggplot2") # simulate 2D ar1 process # pairwise correlation rho <- .95 # pairwise variance sigma <- .5 # 2 dimensions of simulations years <- 20 ages <- 10 # kronnecker product to get joint covariance Q2D <- kronecker(Q.AR1(M=years, sigma, rho), Q.AR1(M=ages, sigma, rho)) # simulate the data and place it in a data frame Q2D.df <- data.frame(obs=c(sim.AR(1, Q2D)), age=rep(1:ages, years), year=rep(1:years, each=ages)) # graph results ggplot(data=Q2D.df, aes(year, obs, group=age, color=age)) + geom_line()
Takes in a square precision matrix of a gaussian Markov random field, which ideally should be sparse, and using algorithm 3.1 from Gaussian Markov Random Fields by Rue & Held simulates data from a mean 0 process where the inverse of the precision matrix represents the variance-covariance of the points in the process. The resulting simulants represent samples of a Gaussian Markov random field (GMRF).
sim.GMRF(n, Q, tol = 1e-12)sim.GMRF(n, Q, tol = 1e-12)
n |
int > 0, number of observations to simulate from the GMRF. |
Q |
matrix, a square precision matrix of a GMRF. |
tol |
numeric, The tolerance for the removal of zero eigenvalues. |
Matrix object, matrix where each row is a single obsrevation from a GMRF with precision structure Q.
Rue, H. Held, L. Gaussian Markov Random Fields: Theory and Applications. Chapman and Hall. 2005.
require("ggplot2") # simulate 2D rw1 process # pairwise variance sigma <- .5 # 2 dimensions of simulations years <- 20 ages <- 10 # kronnecker product to get joint precision Q2D <- kronecker(Q.RW1(M=years, sigma), Q.RW1(M=ages, sigma)) # simulate the data and place it in a data frame Q2D.df <- data.frame(obs=c(sim.GMRF(1, Q2D)), age=rep(1:ages, years), year=rep(1:years, each=ages)) # graph results ggplot(data=Q2D.df, aes(year, obs, group=age, color=age)) + geom_line()require("ggplot2") # simulate 2D rw1 process # pairwise variance sigma <- .5 # 2 dimensions of simulations years <- 20 ages <- 10 # kronnecker product to get joint precision Q2D <- kronecker(Q.RW1(M=years, sigma), Q.RW1(M=ages, sigma)) # simulate the data and place it in a data frame Q2D.df <- data.frame(obs=c(sim.GMRF(1, Q2D)), age=rep(1:ages, years), year=rep(1:years, each=ages)) # graph results ggplot(data=Q2D.df, aes(year, obs, group=age, color=age)) + geom_line()
Spatial Polygons data frame with 475 counties from the US states Louisiana, Texas, Mississippi, & Arkansas. FIPS codes for the state and county are provided in the data frame.
A 475x475 matrix where the index corresponds to a row in the US.df Spatial Polygons data frame and the index of the matrix at row i column j is 1 when US.df[i,] and US.df[j,] share a border and 0 when they do not.