# Function to compute the test statistic u
#' @import fda
compute.u <- function(group1 = NULL, group2 = NULL, dif = NULL, basis, dependent = F){
if(is.null(dif)){
if(any(c(is.null(group1), is.null(group2)))) stop("something's wrong")
else
dif <- group1 - group2
}
stopifnot(ncol(dif$coef) > 0)
if(is.fd(dif))
N <- ncol(dif$coefs)
else
N <- NCOL(dif)
m <- function(t) abs(fda::eval.fd(t, fda::mean.fd(dif)))
s <- function(t) fda::eval.fd(t, fda::sd.fd(dif))
mint <- integrate(m, dif$basis$rangeval[1], dif$basis$rangeval[2])$value
sint <- integrate(s, dif$basis$rangeval[1], dif$basis$rangeval[2])$value
stat <- mint/(sint/sqrt(N))
if(!dependent){
bas2 <- function(t){
sum(fda::eval.basis(t, basis)^2)
}
bas2 <- Vectorize(bas2)
tointegrate <- function(t) sqrt(bas2(t))
intdenom <- integrate(tointegrate, dif$basis$rangeval[1], dif$basis$rangeval[2])$value
sigma <- sint/intdenom
c <- 1
}
else{
# estimate Sigma from data
sigma <- cov(t(dif$coefs))
## WE DO NOT NEED THIS CHUNK OF CODE ANY MORE - DISCARD IF NO ERRORS OCCUR IN FUTURE!
# compute scaling constant
# intexpect <- function(sigma, bas){
# nbas <- bas$nbasis
# sum1 <- function(t, bas, sigma){
# ret <- sum(fda::eval.basis(bas, t)^2*diag(sigma))
# ret
# }
#
# sum2 <- function(t, bas, sigma){
# sumup <- 0
# for(k in 1:nbas){
# for(m in 1:nbas){
# if(k < m) sumup <- sumup + fda::eval.basis(bas, t)[, k] * fda::eval.basis(bas, t)[, m]*sigma[k, m]
# }
# }
# ret <-2*sumup
# ret
# }
#
# sqrt_to_int <- function(t, bas, sigma){
# ret <- sqrt(sum1(t, bas, sigma) + sum2(t, bas, sigma))
# return(ret)
# }
#
# neval <- 100
# sqrteval <- rep(0, neval)
# times <- seq(bas$rangeval[1], bas$rangeval[2], length.out = neval)
# for(i in 1:length(times)){
# sqrteval[i] <- sqrt_to_int(times[i], bas, sigma)
# }
#
# # plot(times, sqrteval, ylim = c(0, 0.5))
#
# ret <- 0
# for(i in 2:neval){
# ret <- ret + 0.5*(times[2] - times[1])*(sqrteval[i] + sqrteval[i-1])
# }
# ret
# }
#
# int <- intexpect(sigma, bas = basis)
# c <- (sint / int)^2
## END OF NOT NEEDED CODE
c <- 1
# scale Sigma
sigma <- sigma*c
}
return(list(stat = stat, sigma = sigma, c = c))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.