Nothing
#' @name SCPT
#' @title Caries out a change-point test for sequences of
#' functional data observed at different spatial locations. Three tests are carried out
#' and three pvalues are returned.
#'
#' @param DATA3 Data for testing. This should be an 3D array with the
#' first coordinate indicating the temporal location on a curve, the
#' second the spatial location, and the third the repitition number/time series iteration.
#' @param loc Two column matrix of spatial locations.
#' @param varprop Proportion of temporal variance explained in first two
#' tests (third test uses all explained variance)
#' @param ncores Number of cores to be used for monte carlo sampling.
#' @param reps_per_core Number of monte carlo samples per core.
#'
#' @details Carries out three different change point tests as outlined in Gromenko, Kokoszka,
#' and Reimherr (2015). P-values are generated using Monte-Carlo. The total MC sample size
#' is given by ncores*reps_per_core. Be warned that the code will not try to use more cores
#' than are available, so users should choose the ncores value appropriately. By default, the
#' code will use all available cores. This will cause a computer to be slugish in other tasks
#' and so the user should be careful when using the default.
#'
#' Data3[i,j,k] indicates an the ith observation at the jth location and on the kth curve. It is
#' important that this order is accurate. The method assumes that the spatio-temporal
#' covariance is separable. Methods such as Constantinou, Kokoszka, and Reimherr (2015) can be
#' used to check this assumption. This allows the spatial and temporal covariances to be
#' estimated separately. The spatial covariance is assumed to be stationary and isotropic. A
#' tapered bsplines estimator is then used to estiamte it nonparametrically. No assumptions are
#' made on the form of the temporal covariance function. It is estimated by pooling appropriately
#' across space.
#'
#' The testing procedures are all based on integrated CUSUM test statistics. The primary difference
#' between the three is in how they handle the temporal covariance. The first statistic uses
#' FPCA to reduce the dimension of the data and normalizes by the corresponding eigenvalues. The
#' second does the same, but does not normalize by the eigenvalue. The third procedure does not use
#' FPCA and does not normalize by the eigenvalues. It can be thought of as the same as procedure 2,
#' but with all dimensions included. The first procedure should be used with caution as it can
#' be unstable for smaller sample sizes/larger numbers of FPCs.
#'
#' @return Returns vector of three values for three change-point tests. The first normalizes each FPC, while the second two do not. The only difference between the second and third is that the third uses all explained variance.
#'
#' @example
#' mean(1:4)
#'
#' @import inline
#' @import Rcpp
#' @import RcppArmadillo
#' @import splines
#' @import parallel
#' @useDynLib scpt
#' @export
spatial_change_test<-function(DATA3,loc,varprop=.85,ncores='all',reps_per_core=100){
#require(parallel); require(RcppArmadillo); require(splines); require(inline)
#source("Simpson.R"); source("Fast_ST.R"); source("Inline.R"); source("Inline_L2.R"); source("Auxiliary_Func.R")
L_per_curve<-dim(DATA3)[1]
N<-dim(DATA3)[2] # number of spatial locations
NY<-dim(DATA3)[3] # number of curves per location
L<-L_per_curve*NY # total temporal observations
DATA.DIFF3<-DATA3[,,2:NY]-DATA3[,,1:(NY-1)]
#c.hat.pr <- get_svar(as.matrix(loc*pi/180), DATA.DIFF3, dim(DATA.DIFF3))
c.hat.pr <- .Call('get_svar',PACKAGE='scpt',as.matrix(loc*pi/180), DATA.DIFF3, dim(DATA.DIFF3))
dd<-c.hat.pr[[1]]
H.hat<-c.hat.pr[[3]]
dd.v<-as.vector(dd)
H.hat.v<-as.vector(H.hat)
#removing zeroes
H.hat.vn<-H.hat.v[dd.v!=0]
dd.vn<-dd.v[dd.v!=0]
#noparametric BSpline estimation
fm1 <- lm(H.hat.vn ~ bs(dd.vn, df = 9))
#truncation
#a<-0.25
#b<-0.27
a<-quantile(dd.vn,.75)
b<-quantile(dd.vn,.95)
slope<-predict(fm1, data.frame(dd.vn = a))/(b-a)
H.final<-matrix(1,N,N)
sp.cov<-predict(fm1, data.frame(dd.vn = dd.vn))
for(i in 1:length(dd.vn)){
H.final[dd.vn[i]==dd]<-sp.cov[i]
if(dd.vn[i] > a & dd.vn[i] < b){
H.final[dd.vn[i]==dd] <- slope*(b-a)-(dd.vn[i]-a)*slope;
}
if(dd.vn[i] > b){
H.final[dd.vn[i]==dd] <- 0;
}
}
H.vec<-as.vector(H.final)
H.vec.n<-H.vec[dd!=0]
dd.n<-dd[dd!=0]
cov.dt<-cbind(dd.n, H.vec.n)
cov.dt<-cov.dt[order(cov.dt[,1]),]
#caclulating weights
c.hm<-matrix(0,N,N)
diag(c.hm)<-c.hat.pr[[2]]
rho.hat.f<-c.hm%*%H.final%*%c.hm
one<-matrix(1, ncol=1, nrow=N)
w<-solve(rho.hat.f^2)%*%one
w<-w/sum(w[,1])
#make_cov_fig()
#temporal covariance
#tvar<-get_tvar(w[,1], c.hat.pr[[2]], DATA.DIFF3, dim(DATA.DIFF3))#may take some time
tvar<-.Call('get_tvar',PACKAGE='scpt',w[,1], c.hat.pr[[2]], DATA.DIFF3, dim(DATA.DIFF3))
FPC<-eigen(tvar[[1]], symmetric=T, only.values = FALSE, EISPACK = FALSE)
val<-FPC$values[FPC$values>0]
val<-val/sum(val)
cvar<-val[1]
for(i in 2:length(val))
cvar<-c(cvar, cvar[i-1]+val[i])
#NFPC<-length(val[cvar<.78])
NFPC<-length(val)
v.hat<-c()
for(i in 1:NFPC)
v.hat<-cbind(v.hat, FPC$vectors[,i]/sqrt(simpson_r(FPC$vectors[,i]^2)))
rho_e<-eigen(rho.hat.f)
rho.hat.f.clip<-(rho_e$vectors[,rho_e$values>0])%*%diag(rho_e$values[rho_e$values>0])%*%t(rho_e$vectors[,rho_e$values>0])
suppressWarnings(Lb<-matrix(t(chol(rho.hat.f.clip,pivot=TRUE)),nrow=N))
cores_all<-detectCores()
if(ncores=='all'||ncores>=cores_all){
XXX<-vector("list",cores_all)
for(i in 1:cores_all)
XXX[[i]]<-reps_per_core
n.cores<-cores_all
}else{XXX<-vector("list",ncores)
for(i in 1:ncores)
XXX[[i]]<-reps_per_core
n.cores<-ncores
}
#XXX<-vector("list",n.cores)
#for(i in 1:n.cores)
# XXX[[i]]<-reps_per_core
# #XXX[[i]]<-250
#tmp<-gn_st(20,Lb,val,w,L_per_curve)
#return(tmp)
#------------------------------------
#generating MC distribution
#------------------------------------
#!!! will use all available cores !!!
#system.time(
#result<-mcmapply(XXX, FUN=gn_st, mc.cores=n.cores,MoreArgs = list(Lb=Lb,val=val,w=w))
result<-mclapply(XXX, FUN=gn_st, mc.cores=n.cores,Lb=Lb,val=val,w=w,L_per_curve = L_per_curve)
#)
#gathering data
SL1<-c()
SL2<-c()
for(i in 1:n.cores){
SL1<-rbind(SL1, result[[i]]$L1)
SL2<-rbind(SL2, result[[i]]$L2)
}
lst<-length(SL2[,1])
nnn<-NFPC
#aaa<-get_stat(DATA3, dim(DATA3), as.matrix(v.hat[,1:nnn]), w[,1], val[1:nnn], nnn)
aaa<-.Call('get_stat',PACKAGE='scpt',DATA3, dim(DATA3), as.matrix(v.hat[,1:nnn]), w[,1], val[1:nnn], nnn)
r1<-aaa[[1]][-1,]
r2<-aaa[[2]][-1,]
L1<-rep(0,nnn)
L2<-rep(0,nnn)
pv1<-rep(0,nnn)
pv2<-rep(0,nnn)
for(i in 1:nnn){
L1[i]<-simpson_r(r1[i,])
L2[i]<-simpson_r(r2[i,])
pv1[i]<-sum((rep(1, lst))[SL1[,i] >= L1[i]])/lst;
pv2[i]<-sum((rep(1, lst))[SL2[,i] >= L2[i]])/lst;
}
spt_for_ret<-which(cumsum(val)>varprop*sum(val))[1]
for_return<-data.frame(pv1[spt_for_ret],pv2[spt_for_ret],pv2[nnn])
colnames(for_return)<-c("Pvalue_Test1", "Pvalue_Test2", "Pvalue_Test2_all")
return(for_return)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.