computeW: Obtain quantities associated with W

In pulongma/DABayes: Bayesian Statistical Model for Climate-Change Detection and Attribution

Description

This function computes quantities associated with the covariance matrix W, and these quantities are used to evaluate the likelihood function for the proposed statistical model.

Usage

computeW(y, x, B, neighborhood_matrix, gamma_prior)

Arguments

y	an n by N matrix of ensembles of measured climate variable (e.g., temperature) increases
x	a list of m elements, each of which is model-output climate variable (e.g., temperature) increases under a specific forcing scenario
B	an n by r basis function matrix
neighborhood_matrix	an n by n matrix
gamma_prior	a vector of possible discrete values for the prior distribution of gamma

Value

a list of 7 elements

Author(s)

Pulong Ma <[e-mail]>

Examples

#################### simulate data ########################
set.seed(1234)
n <- 30 # number of spatial grid cells on the globe
N <- 10 # number of ensemble members
m <- 3 # number of forcing scenarios
Lj <- c(5, 3, 7) # number of runs for each scenario
L0 <- 8 # number of control runs without any external forcing scenario
trend <- 30
DAdata <- simDAdata(n, N, m, Lj, L0, trend)
# ensembles of the measured variable
y <- DAdata[[1]]
# model outputs for the measured variable under different forcing scenarios
x <- DAdata[[2]]
# model outputs for the measured variable without any external forcing scenario
x0 <- DAdata[[3]]
#################### end of simulation ####################



###########################################################
### preprocessing 
# center the data
y <- y - mean(y)
for(j in 1:m){
  x[[j]] <- x[[j]]-mean(x[[j]])
}
# construct basis function matrix B with principal components
r <- 3
empiricalcovmat <- cov(t(x0))

# compute first r eigenvalues and eigenvectors
temp <- RSpectra::eigs_sym(empiricalcovmat, r, which="LM")
B <- temp$vectors
K_hat <- temp$values


lambda_mean <- log(K_hat)
lambda_var <- var(lambda_mean)/3^2


### pre-computation for W

# possible values for the discrete prior of gamma
gamma_prior <- c(0.5,0.99,1,1.01,2,5)
ng <- length(gamma_prior)

# neighborhood matrix (replace by real one)
neighborhood_matrix <- Matrix::sparseMatrix(1:n, 1:n, x=rep(1,n))
## Not run: 
outW <- computeW(y,x,B,neighborhood_matrix,gamma_prior)

## End(Not run)

pulongma/DABayes documentation built on June 24, 2019, 12:38 a.m.