Nothing
#' Compute the Robust g test statistic
#'
#' @description
#' The Robust g test statistic is computed for testing for periodicity.
#'
#' @param y Vector containing the series.
#' @param freqSet The set to search frequencies.
#'
#' @return Vector of length 2 containing the Robust g test statistic
#' and the frequency for the maximum periodogram.
#'
#' @details \code{An} refers to searching frequencies in the range
#' twice as dense as the region of the fourier frequencies, i.e.,
#' \eqn{An = {0,1/n,\ldots,floor((2 n - 1)/2)/n}}.
#' \code{En} refers to searching frequencies in
#' \eqn{En = {j/101 | j=1,\dots,50 and j/101 \ge 1/n}}.
#'
#' @author Yuanhao Lai
#'
#' @references Ahdesmaki, M., Lahdesmaki, H., Pearson, R., Huttunen, H.,
#' and Yli-Harja O.(2005).
#' \emph{BMC Bioinformatics} \bold{6}:117.
#' \url{http://www.biomedcentral.com/1471-2105/6/117}
#'
#' @keywords internal
# These programs are from the R package 'GeneCycle', but they are modified a
# little bit for simplification. Compute the Robust g test statistic
GetFitRobustG <- function(y,freqSet = c("An","En")) {
freqSet <- match.arg(freqSet)
Dpgram <- rankSpectrum(y,freqSet)
Dpgram <- Dpgram[-1,]
maxL <- which.max(Dpgram[,2])
gr <- Dpgram[maxL,2]/sum(Dpgram[,2])
ans <- c(gr,Dpgram[maxL,1]) #statistic and frequency
names(ans) <- c("gr","freq")
ans
}
####################################################
# This function implements the robust, rank-based spectral
# estimator introduced in Pearson et al. 2003
#####################################################
#' The robust rank-based spectral
#'
#' @description This function implements the robust, rank-based spectral.
#'
#' @param x Series.
#' @param freqSet The set to search frequencies.
#'
#' @return A numerical value of the robust rank-based spectral.
#'
#' @keywords internal
rankSpectrum <- function(x,freqSet = c("An","En"))
{
freqSet <- match.arg(freqSet)
##############################################
# Some adjustable parameters
##############################################
# Length of the original sequence
n <- length(x)
# Let us define the maximum lag for the correlation coefficient:
maxM <- n-2
##############################################
# Correlation coefficient
##############################################
# Reserve space
Rsm <- matrix(NA, nrow = maxM+1, ncol = 1)
nonMissing <- rep(TRUE,n)
nmi <- 1:n # indices for nonmissing values
# Mean removal
x <- x- mean(x)
# Run through all the lags
for (lags in 0:maxM)
{
# Modified Spearman's method
indexes <- 1:(n-lags) # Initial indices
ends <- n
# Values in both the original and shifted vectors must be present:
temp <- (nonMissing[1:(ends-lags)] + nonMissing[(lags+1):ends]) >= 2
indexPresent <- which(temp)
indexes <- indexes[indexPresent] # The indices that are present in
# both sequences
Rsm[lags+1] <- ifelse(
length(indexes)<=1 ,
0 ,
cor(x[indexes], x[indexes+lags], method="spearman" ) * length( x[indexes] )/n)
if(is.na(Rsm[lags+1])) Rsm[lags+1]<-0
}
# Zero-padding
# Length of the zero-padded one-sided "rho"(=Rsm)
if(freqSet=="An"){
zp <- 2*length(x)
Rsm[(length(Rsm)+1):zp] <- 0
fftemp <- fft(Rsm)[1:floor(zp/2)]
freq <- (1:floor(zp/2)-1)/zp
}else if(freqSet=="En"){
zp <- 101
Rsm <- c(Rsm,rep(0,zp-length(Rsm)))
fftemp <- fft(Rsm)
fftemp <- c(fftemp[1],fftemp[(ceiling(zp/n)+1):51])
freq <- c(0,(ceiling(zp/n):50)/101)
}
# The following implementation is as in (Ahdesmaki, Lahdesmaki et al., 2005)
Ssm <- abs( 2*Re(fftemp) - Rsm[1] )
# Return the spectral content, frequencies [0,pi)
ans <- cbind(freq,Ssm)
colnames(ans) <- c("frequency", "robustPeriodogram")
return(ans)
}
# n <- 1000
# x <- rnorm(n,sd =5)
# rx <- rankSpectrum(x)
# max(rx)
# sum(rx)
#
# var(x)
# mad(x)
#
# var(1:n)
#
# xl <- seq(10,100,1)
# yl <- sapply(xl,FUN = function(n){
# rx <- rankSpectrum(rnorm(n))
# sum(rx)
# })
#
# plot(xl,yl)
#
# summary(lm(yl~xl))
#
#
# yl2 <- sapply(xl,FUN = function(n){
# var(1:n)
# })
#
# plot(xl,yl2)
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.