computeW: Obtain quantities associated with W
In mapn/DAbayes: Bayesian Statistical Model for Climate-Change Detection and Attribution

Description Usage Arguments Value Author(s) Examples

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.

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

`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

a list of 7 elements

Pulong Ma <[email protected]>

#################### 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)