options(markdown.HTML.stylesheet = "path/to/a/custom/style.css")
source("../R/utilities.R")
Rcpp::sourceCpp("../src/utilities.cpp")
library(MASS)

Create data under H0 and HA

set.seed(2019)
n = 30
k = 3
nullX = as.matrix(rnorm(30), ncol=1)
nullY = MASS::mvrnorm(n, rep(0,3), diag(3))
set.seed(2019)
n = 30
k = 3
Y = matrix(NA, n, k)
X = rnorm(30)
rho = X * 0.2
rho = rho - min(rho)
X = matrix(X, ncol=1)
for (i in 1:n){
  Sigma = matrix(rho[i], k, k)
  diag(Sigma) = 1
  Y[i,] = MASS::mvrnorm(1, rep(0,k), Sigma)
}

Score test under H0

q = get_score_c(nullX, nullY[,1], nullY[,2])
p = 1-pchisq(q, 1)
print(paste("The score statistic is", round(q,2), 
            "and the p-value is", round(p,2)))

Score test under HA

q = get_score_c(X, Y[,1], Y[,2])
p = 1-pchisq(q, 1)
print(paste("The score statistic is", round(q,2), 
            "and the p-value is", round(p,2)))

Degree test under H0

d = get_degree_c(as.matrix(nullX, ncol=1), nullY[,1], nullY[,2:3])
p = get_p_from_degree_c(nullY[,1], nullY[,2:3], d)
print(paste("The degree statistic is", round(d,2), 
            "and the p-value is", round(p,2)))

Degree test under HA

d = get_degree_c(as.matrix(X, ncol=1), Y[,1], Y[,2:3])
p = get_p_from_degree_c(Y[,1], Y[,2:3], d)
print(paste("The degree statistic is", round(d,2), 
            "and the p-value is", round(p,2)))


tk382/covarianceBP documentation built on July 14, 2020, 4:01 p.m.