R functions for scalable subsampling aggregation

##########################################################
SSUBBAG = 
function(x, beta=0.8, c1=1,  statistic, ...)
{ # scalable subsampling aggregation
  # algorithm uses h=c1*b; c1=1 yields non-overlapping blocks
      n <- length(x)
      b <- round(0.5*n^beta)
      h <- round( c1*b)
      q <- floor( (n-b)/h )+1
      thetahat <- c(1:q) * 0
      A=matrix(c(1:(b*q) ),nrow= b, ncol=q ) 
      for(i in 1:q) {
A[,i] <- as.vector(x[((i - 1) * h + 1): ((i - 1) * h + b) ])}
      thetahat=apply(A,2, statistic, ...)
      # apply function mimicks distributed with q machines
      SS=mean(thetahat)
      SS
}
##########################################################
density.atzero = 
function(x, undersmooth=F)
{
aa=density.estimate(x,undersmooth)
bb=lsfit(c(1:512),aa$x)
bb=bb$coefficients
c0=floor(-bb[1]/bb[2])

return(aa$y[c0])
}
##########################################################
density.estimate =
function(x, undersmooth=F)
{
n <- length(x)
if(undersmooth==F){bb = 1.5*sd(x)*n^(-0.2)} #2
if(undersmooth==T){bb = 1.5*sd(x)*n^(-0.25)} #1
      # bandwidth choice
aa= density(x, bw = bb, kernel = "gaussian" )
return(aa)
# can do: plot(aa$x, aa$y, type="lines")
}
##########################################################
spec.density.atzero = 
function(x, undersmooth=F){
aa=spec.parzen(x,undersmooth)
return(aa[1,2])
}
##########################################################
spec.parzen =
function(x, undersmooth=F)
{
# estimate spectral density of time series x using Parzen window
# smoothness is inversely proportional to parameter M

n <- length(x)# w is a subdivision of the interval [0, pi]
              # near 0 is low freq., near pi is highest freq.
if(undersmooth==F){M=1+7*floor(n^0.2)}  
if(undersmooth==T){M=1+7*floor(n^0.25)} 

w <- (3.1416 * c(0:(n  - 1)))/n 
f <- c(1:length(w)) * 0
L <- lambda.parzen(M)
xcov <- acf(x, lag.max = (M - 1), type = "covariance", plot = F)$acf
xcov <- as.vector(xcov)
for(i in 1:length(w)) {
cosines <- cos(w[i] * c(0:(M - 1)))
f[i] <-  0.1591546 * ( - xcov[1] + 2 * sum(L * xcov * cosines)) 
}
fw <- cbind(w, f)# to plot, just do: plot(fw, type="l")
fw
}
##########################################################
lambda.parzen =
function(M)
{y=c(1:M)-1
L=y
for (i in 1:M) {L[i]= parzen.taper (y[i]/M)}
L }
##########################################################
parzen.taper = 
function(x)
{x <- abs(x)
val <- 0
if(x < .5)
{val <- 1 - 6*x^2 + 6*x^3} 
else
{if(x < 1) {val <- 2*(1-x)^3}
}
return(val)
}
##########################################################