karthik · December 7, 2012 21:34
diff --git a/mvspec2b_source.R b/mvspec2b_source.R
 #'mvspec2b
 #'
 #' Description not yet complete since function isn't.
 #'@param x matrix of (multi-variate) time-series in columns.  For 'pgram' - which assumes even sampling - x must be a time-series (ts) object with the sampling frequency specified.  Although 'lombscargle' will also work on a ts object, a ts object assumes even sampling and therefore can't be used.  If 'lombscargle' is used on evenly sampled data, virtually the same result as the 'pgram' method will be returned.
 #'@param  type Method used to estimate spectral density.  Options are 'pgram' or 'lombscargle' (default) or 'fft' (from cts package - not checked for correctness).
 #'@param  time_s (optional, required for Lomb-Scargle) Vector of time-points.  time_s points may be unevenly spaced, but for multi-species time-series are assumed to the same for all species.
 #'@param spans  <what param does>
 #'@param kernel <what param does>
 #'@param taper = 0.1 <what param does>
 #'@param pad = 0 <what param does>
 #'@param fast  <what param does>
 #'@param demean <what param does>
 #'@param detrend <what param does>
 #'@param scale T/F(default) - normalize each time series to its standard deviation (should be run with either demean=TRUE or detrend=TRUE).
 #'@param normalize Default is \code{FALSE}
 #'@param logt T/F - Transform data by adding 1 and log-transforming.  This is performed before additional transformations
 #'@param plot <wMusic For Programming
 #'@param na.action = na.pass <what param does>
 #'@param Smethod (optional) Method to stabilize the spectral matrix before inversion in calculating partial coherence.  Specify either 'none','upweighted' (the default..until shrinkage is improved), or'shrinkage' to use either diagonal up-weighting or shrinkage estimation.
 #'@param zeta (optional) diagonal up-weighting parameter.  (Medkour et al. use values of 10^-6 to 10^-2.)
 #'@param LSup <what param does>
 #'@param ... <what param does>
 #'@return In addition to that of stats:spec.pgram, also calulates each frequency's sprectral matrix and partial coherence.
 #'@export
 #'@examples \dontrun{
 #' mvspec2b(foo) # where foo is an object of class AggMean.
 #'}
 mvspec2b <- function(x, type = "lombscargle", time_s = NULL, spans = 3,
    kernel = NULL, taper = 0.1, pad = 0, fast = TRUE, demean = FALSE, detrend = TRUE,
    scale = FALSE, normalize = FALSE, logt = FALSE, plot = FALSE, na.action = na.pass, Smethod = "upweighted",
    zeta = 10^-2, LSup = NULL, ...) {
    message(sprintf("normalize %s, logt %s, demean %s \n", normalize, logt, demean))

    series <- deparse(substitute(x))
    x <- na.action(as.ts(x))
    xfreq <- frequency(x)
    if (type == "pgram" & sum(is.na(x)) == TRUE) {
        x <- na.omit(x)
        if (!is.null(na.action(x))) {
            warning("Some time points were removed due to missing variables.  To avoid losing points, try using \"lomb-scargle\" method instead.")
        }
    }
    ts <- convert_to_lon(x)
    time_s <- attributes(ts$ldata)$time
    if (!is.null(time_s)) {
        tn <- time_s
        xfreq <- 1/min(diff(tn))
    }
    if (type == "lombscargle" && is.null(time_s)) {
        stop("Must provide time-points for Lomb-Scargle method")
    }
    N <- N0 <- nrow(x)
    nser <- ncol(x)
    if (!is.null(spans))
        kernel <- {
            if (is.tskernel(spans))
                spans else kernel("modified.daniell", spans%/%2)
        }
    if (!is.null(kernel) && !is.tskernel(kernel))
        stop("must specify spans or a valid kernel")
    if (logt) {
        message("Going into log")
        x <- log(x + 1)
    }
    if (detrend) {
        if (type != "lombscargle") {
            t <- 1:N - (N + 1)/2
            sumt2 <- N * (N^2 - 1)/12
        }
        if (type == "lombscargle") {
            t <- tn - mean(tn)
            sumt2 <- sum(t^2)
        }
        for (i in 1:ncol(x)) x[, i] <- x[, i] - mean(x[, i], na.rm = TRUE) - sum(x[,
            i] * t, na.rm = TRUE) * t/sumt2
    } else if (demean) {
        x <- sweep(x, 2, colMeans(x, na.rm = TRUE), check.margin = FALSE)
    }
    if (scale) {
        x <- scale(x, center = FALSE, scale = apply(x, 2, sd, na.rm = TRUE))
    }

    if (normalize) {
        x <- apply(x, 2, function(z) { (z - min(z, na.rm = TRUE)) /
        (max(z, na.rm = TRUE) - min(z, na.rm = TRUE)) } )
        }

    x <- spec.taper(x, taper)
    u2 <- (1 - (5/8) * taper * 2)
    u4 <- (1 - (93/128) * taper * 2)
    if (pad > 0) {
        x <- rbind(x, matrix(0, nrow = N * pad, ncol = ncol(x)))
        N <- nrow(x)
    }
    NewN <- if (fast)
        nextn(N) else N
    x <- rbind(x, matrix(0, nrow = (NewN - N), ncol = ncol(x)))
    N <- nrow(x)
    Nspec <- floor(N/2)
    freq <- seq(from = xfreq/N, by = xfreq/N, length = Nspec)
    pgram <- array(NA, dim = c(N, ncol(x), ncol(x)))
    if (type == "pgram") {
        xfft <- mvfft(x)
        for (i in 1:ncol(x)) {
            for (j in 1:ncol(x)) {
                pgram[, i, j] <- xfft[, i] * Conj(xfft[, j])/(N0 * xfreq)
                pgram[1, i, j] <- 0.5 * (pgram[2, i, j] + pgram[N, i, j])
            }
        }
    }
    if (type == "lombscargle") {
                                	# Lomb-Scargle Discrete Fourier Transform (proposed by Scargle 1982, 1989).  This implementation follows Schulz and Stattegger (1997) eqns 1a-e, expcept that it does not implement the Welch-Overlapped-Segment-Averaging method (since the time series we plan to analyse are relatively short).
        freq.temp <- seq(from = xfreq/N, by = xfreq/N, length = N)
        Xk <- array(NA, dim = c(N, ncol(x)))
        for (i in 1:ncol(x)) {
            for (k in 1:length(freq.temp)) {
                omega <- 2 * pi * freq.temp[k]
                tao <- 1/(2 * omega) * atan(sum(sin(2 * omega * tn))/sum(cos(2 *
                  omega * tn)))
                tf <- tn[1] - tao  # assuming all variables are sampled at same time points
                F0 <- (1/sqrt(2)) * exp((0 + (0 + (0 + (0+1i)))) * omega * (tf))
                A <- (sum((cos(omega * (tn - tao)))^2))^(-0.5)
                B <- (sum((sin(omega * (tn - tao)))^2))^(-0.5)
                Xk[k, i] <- F0 * sum(A * x[1:length(tn), i] * cos(omega * (tn - tao)) +
                  (0 + (0 + (0 + (0+1i)))) * B * x[1:length(tn), i] * sin(omega *
                    (tn - tao)), na.rm = TRUE)
            }
        }
        for (i in 1:ncol(x)) {
            for (j in 1:ncol(x)) {
                pgram[, i, j] <- Xk[, i] * Conj(Xk[, j])/(N0 * xfreq)
                pgram[1, i, j] <- 0.5 * (pgram[2, i, j] + pgram[N, i, j])
            }
        }
    }
                # from cts package:
    if (type == "fft") {
        for (i in 1:ncol(x)) {
            for (j in 1:ncol(x)) {
                for (k in 1:length(freq.temp)) omega <- 2 * pi * freq.temp[k]
                pgram[k, i, j] <- ((sum(x[1:length(tn), i] * cos(omega * tn)))^2 +
                  (sum(x[1:length(tn), i] * sin(omega * tn)))^2)/N
            }
        }
    }
    rpgram <- pgram  # save raw periodogram for shrinkage estimation
    if (!is.null(kernel)) {
        for (i in 1:ncol(x)) {
            for (j in 1:ncol(x)) {
                pgram[, i, j] <- kernapply(pgram[, i, j], kernel, circular = TRUE)
            }
        }
        df <- df.kernel(kernel)
        bandwidth <- bandwidth.kernel(kernel)
    } else {
        df <- 2
        bandwidth <- sqrt(1/12)
    }
    df <- df/(u4/u2^2)
    df <- df * (N0/N)
    bandwidth <- bandwidth * xfreq/N
                    # Not sure why following should differ, but they return correct frequencies for each.
    if (type == "pgram")
        {
            pgram <- pgram[2:(Nspec + 1), , , drop = FALSE]
        }  # from spec.pgram
    if (type != "pgram")
        {
            pgram <- pgram[1:(Nspec), , , drop = FALSE]
        }  # from cts:spec.ls
    spec <- matrix(NA, nrow = Nspec, ncol = nser)
    for (i in 1:nser) {
        spec[, i] <- Re(pgram[1:Nspec, i, i])
    }
    if (nser == 1) {
        coh <- phase <- NULL
    } else {
        coh <- phase <- matrix(NA, nrow = Nspec, ncol = nser * (nser - 1)/2)
        for (i in 1:(nser - 1)) {
            for (j in (i + 1):nser) {
                coh[, i + (j - 1) * (j - 2)/2] <- Mod(pgram[, i, j])^2/(spec[, i] *
                  spec[, j])
                phase[, i + (j - 1) * (j - 2)/2] <- Arg(pgram[, i, j])
            }
        }
    }
    for (i in 1:nser) spec[, i] <- spec[, i]/u2
    spec <- drop(spec)
                #=============
    rfxxs <- array(NA, dim = c(nser, nser, Nspec))
    if (Smethod == "shrinkage") {
        rdf <- 2/(u4/u2^2) * (N0/N)
        rbandwidth <- sqrt(1/12) * xfreq/N
        rpgram <- rpgram[2:(Nspec + 1), , , drop = FALSE]
                                # 		rspec <- matrix(NA, nrow = Nspec, ncol = nser)
                                # 		for (i in 1:nser){ rspec[, i] <- Re(rpgram[1:Nspec, i, i]) }
        for (k in 1:Nspec) {
            rfxxs[, , k] <- Re(rpgram[k, , ])
        }
    }
                #========================
                # Stoffer's addition of spectral matrix and MN's addition of partial coherence
    pcoh <- fxx <- array(NA, dim = c(nser, nser, Nspec))
    shrink.wt <- Kappa <- rep(NA, length(Nspec))
    for (k in 1:Nspec) {
        fxx[, , k] <- Re(pgram[k, , ])  # Not sure this should be real part only, but fails for shrinkage method otherwise
        if (nser == 1) {
            pcoh <- shrink.wt <- Kappa <- NULL
        } else {
            pcoh.t <- pcoh(fxx[, , k], k, Smethod = Smethod, zeta = zeta, span = spans,
                rfxxs = rfxxs)
            pcoh[, , k] <- pcoh.t$pcoh
            shrink.wt[k] <- pcoh.t$shrink.wt
            Kappa[k] <- pcoh.t$Kappa
            if (Smethod != "shrinkage") {
                shrink.wt <- NULL
            }
        }
    }
                #========================
    spg.out <- list(freq = freq, spec = spec, coh = coh, phase = phase, pcoh = pcoh,
        kernel = kernel, df = df, bandwidth = bandwidth, n.used = N, fxx = fxx, rfxxs = rfxxs,
        orig.n = N0, series = series, snames = colnames(x), type = type, method = ifelse(!is.null(kernel),
            "Smoothed Periodogram", "Raw Periodogram"), taper = taper, pad = pad,
        detrend = detrend, demean = demean, Smethod = Smethod, zeta = zeta, shrink.wt = shrink.wt,
        Kappa = Kappa)
    class(spg.out) <- "spec"
    return(spg.out)
 }
	#'mvspec2b
	#'
	#' Description not yet complete since function isn't.
	#'@param x matrix of (multi-variate) time-series in columns. For 'pgram' - which assumes even sampling - x must be a time-series (ts) object with the sampling frequency specified. Although 'lombscargle' will also work on a ts object, a ts object assumes even sampling and therefore can't be used. If 'lombscargle' is used on evenly sampled data, virtually the same result as the 'pgram' method will be returned.
	#'@param type Method used to estimate spectral density. Options are 'pgram' or 'lombscargle' (default) or 'fft' (from cts package - not checked for correctness).
	#'@param time_s (optional, required for Lomb-Scargle) Vector of time-points. time_s points may be unevenly spaced, but for multi-species time-series are assumed to the same for all species.
	#'@param spans <what param does>
	#'@param kernel <what param does>
	#'@param taper = 0.1 <what param does>
	#'@param pad = 0 <what param does>
	#'@param fast <what param does>
	#'@param demean <what param does>
	#'@param detrend <what param does>
	#'@param scale T/F(default) - normalize each time series to its standard deviation (should be run with either demean=TRUE or detrend=TRUE).
	#'@param normalize Default is \code{FALSE}
	#'@param logt T/F - Transform data by adding 1 and log-transforming. This is performed before additional transformations
	#'@param plot <wMusic For Programming
	#'@param na.action = na.pass <what param does>
	#'@param Smethod (optional) Method to stabilize the spectral matrix before inversion in calculating partial coherence. Specify either 'none','upweighted' (the default..until shrinkage is improved), or'shrinkage' to use either diagonal up-weighting or shrinkage estimation.
	#'@param zeta (optional) diagonal up-weighting parameter. (Medkour et al. use values of 10^-6 to 10^-2.)
	#'@param LSup <what param does>
	#'@param ... <what param does>
	#'@return In addition to that of stats:spec.pgram, also calulates each frequency's sprectral matrix and partial coherence.
	#'@export
	#'@examples \dontrun{
	#' mvspec2b(foo) # where foo is an object of class AggMean.
	#'}
	mvspec2b <- function(x, type = "lombscargle", time_s = NULL, spans = 3,
	kernel = NULL, taper = 0.1, pad = 0, fast = TRUE, demean = FALSE, detrend = TRUE,
	scale = FALSE, normalize = FALSE, logt = FALSE, plot = FALSE, na.action = na.pass, Smethod = "upweighted",
	zeta = 10^-2, LSup = NULL, ...) {
	message(sprintf("normalize %s, logt %s, demean %s \n", normalize, logt, demean))

	series <- deparse(substitute(x))
	x <- na.action(as.ts(x))
	xfreq <- frequency(x)
	if (type == "pgram" & sum(is.na(x)) == TRUE) {
	x <- na.omit(x)
	if (!is.null(na.action(x))) {
	warning("Some time points were removed due to missing variables. To avoid losing points, try using \"lomb-scargle\" method instead.")
	}
	}
	ts <- convert_to_lon(x)
	time_s <- attributes(ts$ldata)$time
	if (!is.null(time_s)) {
	tn <- time_s
	xfreq <- 1/min(diff(tn))
	}
	if (type == "lombscargle" && is.null(time_s)) {
	stop("Must provide time-points for Lomb-Scargle method")
	}
	N <- N0 <- nrow(x)
	nser <- ncol(x)
	if (!is.null(spans))
	kernel <- {
	if (is.tskernel(spans))
	spans else kernel("modified.daniell", spans%/%2)
	}
	if (!is.null(kernel) && !is.tskernel(kernel))
	stop("must specify spans or a valid kernel")
	if (logt) {
	message("Going into log")
	x <- log(x + 1)
	}
	if (detrend) {
	if (type != "lombscargle") {
	t <- 1:N - (N + 1)/2
	sumt2 <- N * (N^2 - 1)/12
	}
	if (type == "lombscargle") {
	t <- tn - mean(tn)
	sumt2 <- sum(t^2)
	}
	for (i in 1:ncol(x)) x[, i] <- x[, i] - mean(x[, i], na.rm = TRUE) - sum(x[,
	i] * t, na.rm = TRUE) * t/sumt2
	} else if (demean) {
	x <- sweep(x, 2, colMeans(x, na.rm = TRUE), check.margin = FALSE)
	}
	if (scale) {
	x <- scale(x, center = FALSE, scale = apply(x, 2, sd, na.rm = TRUE))
	}

	if (normalize) {
	x <- apply(x, 2, function(z) { (z - min(z, na.rm = TRUE)) /
	(max(z, na.rm = TRUE) - min(z, na.rm = TRUE)) } )
	}

	x <- spec.taper(x, taper)
	u2 <- (1 - (5/8) * taper * 2)
	u4 <- (1 - (93/128) * taper * 2)
	if (pad > 0) {
	x <- rbind(x, matrix(0, nrow = N * pad, ncol = ncol(x)))
	N <- nrow(x)
	}
	NewN <- if (fast)
	nextn(N) else N
	x <- rbind(x, matrix(0, nrow = (NewN - N), ncol = ncol(x)))
	N <- nrow(x)
	Nspec <- floor(N/2)
	freq <- seq(from = xfreq/N, by = xfreq/N, length = Nspec)
	pgram <- array(NA, dim = c(N, ncol(x), ncol(x)))
	if (type == "pgram") {
	xfft <- mvfft(x)
	for (i in 1:ncol(x)) {
	for (j in 1:ncol(x)) {
	pgram[, i, j] <- xfft[, i] * Conj(xfft[, j])/(N0 * xfreq)
	pgram[1, i, j] <- 0.5 * (pgram[2, i, j] + pgram[N, i, j])
	}
	}
	}
	if (type == "lombscargle") {
	# Lomb-Scargle Discrete Fourier Transform (proposed by Scargle 1982, 1989). This implementation follows Schulz and Stattegger (1997) eqns 1a-e, expcept that it does not implement the Welch-Overlapped-Segment-Averaging method (since the time series we plan to analyse are relatively short).
	freq.temp <- seq(from = xfreq/N, by = xfreq/N, length = N)
	Xk <- array(NA, dim = c(N, ncol(x)))
	for (i in 1:ncol(x)) {
	for (k in 1:length(freq.temp)) {
	omega <- 2 * pi * freq.temp[k]
	tao <- 1/(2 * omega) * atan(sum(sin(2 * omega * tn))/sum(cos(2 *
	omega * tn)))
	tf <- tn[1] - tao # assuming all variables are sampled at same time points
	F0 <- (1/sqrt(2)) * exp((0 + (0 + (0 + (0+1i)))) * omega * (tf))
	A <- (sum((cos(omega * (tn - tao)))^2))^(-0.5)
	B <- (sum((sin(omega * (tn - tao)))^2))^(-0.5)
	Xk[k, i] <- F0 * sum(A * x[1:length(tn), i] * cos(omega * (tn - tao)) +
	(0 + (0 + (0 + (0+1i)))) * B * x[1:length(tn), i] * sin(omega *
	(tn - tao)), na.rm = TRUE)
	}
	}
	for (i in 1:ncol(x)) {
	for (j in 1:ncol(x)) {
	pgram[, i, j] <- Xk[, i] * Conj(Xk[, j])/(N0 * xfreq)
	pgram[1, i, j] <- 0.5 * (pgram[2, i, j] + pgram[N, i, j])
	}
	}
	}
	# from cts package:
	if (type == "fft") {
	for (i in 1:ncol(x)) {
	for (j in 1:ncol(x)) {
	for (k in 1:length(freq.temp)) omega <- 2 * pi * freq.temp[k]
	pgram[k, i, j] <- ((sum(x[1:length(tn), i] * cos(omega * tn)))^2 +
	(sum(x[1:length(tn), i] * sin(omega * tn)))^2)/N
	}
	}
	}
	rpgram <- pgram # save raw periodogram for shrinkage estimation
	if (!is.null(kernel)) {
	for (i in 1:ncol(x)) {
	for (j in 1:ncol(x)) {
	pgram[, i, j] <- kernapply(pgram[, i, j], kernel, circular = TRUE)
	}
	}
	df <- df.kernel(kernel)
	bandwidth <- bandwidth.kernel(kernel)
	} else {
	df <- 2
	bandwidth <- sqrt(1/12)
	}
	df <- df/(u4/u2^2)
	df <- df * (N0/N)
	bandwidth <- bandwidth * xfreq/N
	# Not sure why following should differ, but they return correct frequencies for each.
	if (type == "pgram")
	{
	pgram <- pgram[2:(Nspec + 1), , , drop = FALSE]
	} # from spec.pgram
	if (type != "pgram")
	{
	pgram <- pgram[1:(Nspec), , , drop = FALSE]
	} # from cts:spec.ls
	spec <- matrix(NA, nrow = Nspec, ncol = nser)
	for (i in 1:nser) {
	spec[, i] <- Re(pgram[1:Nspec, i, i])
	}
	if (nser == 1) {
	coh <- phase <- NULL
	} else {
	coh <- phase <- matrix(NA, nrow = Nspec, ncol = nser * (nser - 1)/2)
	for (i in 1:(nser - 1)) {
	for (j in (i + 1):nser) {
	coh[, i + (j - 1) * (j - 2)/2] <- Mod(pgram[, i, j])^2/(spec[, i] *
	spec[, j])
	phase[, i + (j - 1) * (j - 2)/2] <- Arg(pgram[, i, j])
	}
	}
	}
	for (i in 1:nser) spec[, i] <- spec[, i]/u2
	spec <- drop(spec)
	#=============
	rfxxs <- array(NA, dim = c(nser, nser, Nspec))
	if (Smethod == "shrinkage") {
	rdf <- 2/(u4/u2^2) * (N0/N)
	rbandwidth <- sqrt(1/12) * xfreq/N
	rpgram <- rpgram[2:(Nspec + 1), , , drop = FALSE]
	# rspec <- matrix(NA, nrow = Nspec, ncol = nser)
	# for (i in 1:nser){ rspec[, i] <- Re(rpgram[1:Nspec, i, i]) }
	for (k in 1:Nspec) {
	rfxxs[, , k] <- Re(rpgram[k, , ])
	}
	}
	#========================
	# Stoffer's addition of spectral matrix and MN's addition of partial coherence
	pcoh <- fxx <- array(NA, dim = c(nser, nser, Nspec))
	shrink.wt <- Kappa <- rep(NA, length(Nspec))
	for (k in 1:Nspec) {
	fxx[, , k] <- Re(pgram[k, , ]) # Not sure this should be real part only, but fails for shrinkage method otherwise
	if (nser == 1) {
	pcoh <- shrink.wt <- Kappa <- NULL
	} else {
	pcoh.t <- pcoh(fxx[, , k], k, Smethod = Smethod, zeta = zeta, span = spans,
	rfxxs = rfxxs)
	pcoh[, , k] <- pcoh.t$pcoh
	shrink.wt[k] <- pcoh.t$shrink.wt
	Kappa[k] <- pcoh.t$Kappa
	if (Smethod != "shrinkage") {
	shrink.wt <- NULL
	}
	}
	}
	#========================
	spg.out <- list(freq = freq, spec = spec, coh = coh, phase = phase, pcoh = pcoh,
	kernel = kernel, df = df, bandwidth = bandwidth, n.used = N, fxx = fxx, rfxxs = rfxxs,
	orig.n = N0, series = series, snames = colnames(x), type = type, method = ifelse(!is.null(kernel),
	"Smoothed Periodogram", "Raw Periodogram"), taper = taper, pad = pad,
	detrend = detrend, demean = demean, Smethod = Smethod, zeta = zeta, shrink.wt = shrink.wt,
	Kappa = Kappa)
	class(spg.out) <- "spec"
	return(spg.out)
	}