R/plotsHelper.R

Defines functions mscToCDFquantiles cdfToMSqCoh paxpt7 FtoMSC C2toF trnrm llftr7 lftr3p lplotDefault lplotSpec plotFtest

##     The multitaper R package
##     Multitaper and spectral analysis package for R
##     Copyright (C) 2013 Karim Rahim 
##
##     Written by Karim Rahim.
##
##     This file is part of the multitaper package for R.
##     http://cran.r-project.org/web/packages/multitaper/index.html
## 
##     The multitaper package is free software: you can redistribute it and 
##     or modify it under the terms of the GNU General Public License as 
##     published by the Free Software Foundation, either version 2 of the 
##     License, or any later version.
##
##     The multitaper package is distributed in the hope that it will be 
##     useful, but WITHOUT ANY WARRANTY; without even the implied warranty 
##     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##     GNU General Public License for more details.
##
##     You should have received a copy of the GNU General Public License
##     along with multitaper.  If not, see <http://www.gnu.org/licenses/>.
##
##
##     If you wish to report bugs please contact the author:
## 
##     Karim Rahim
##     [email protected]


## siglines degrees of freedom correction Oct 4, 2012 karim

##################################################################
##
##  .plotFtest
##
##  Takes a mtm object, and plots the harmonic F-test statistic 
##  (obj$Ftest).
##
##################################################################
.plotFtest <- function(x, 
                       ftbase=1.01, 
                       siglines=NULL, 
                       xlab="Frequency", 
                       ...) {

    if(is.null(x$mtm$Ftest) || !("Ftest" %in% class(x))) {
      stop(paste("Ftest not computed for given mtm object!"))
    }

    ## correct notes and warnings with devel R version to upload to cran.
    ##arglist <- list(...)
    ##log <- arglist$log ##match.call(expand.dots = )$log
    ##ylab <- arglist$ylab ##match.call(expand.dots = )$ylab
    log <- match.call(expand.dots = )$log
    ylab <- match.call(expand.dots = )$ylab
    
    if(is.null(ylab)) ylab <- "Harmonic F-test Statistic"
    
    ylog = "n"
    if(is.null(log) || log == "yes") {
        ylog = "y"
    }

    ftestVals = x$mtm$Ftest
    ftestVals[ftestVals < ftbase] <- ftbase
    ftmax <- max(ftestVals)

    .lplotDefault(x$freq, ftestVals, log=ylog, ylab=ylab, xlab=xlab, 
                 ylim=c(ftbase,ftmax), type="l", ...)
    
    ## add siglines if defined
    if(!is.null(siglines)) {
        for(j in 1:length(siglines)) {
            if(is.numeric(siglines[j]) && 0.80 <= siglines[j] && 1.000000 >= siglines[j]) {
                ## degree of freedom correction P&W page 499 changed to 2, 2*k-2 date Sept 30 2012
                sig0 <- qf(siglines[j],2,2*x$mtm$k-2) 
                abline(h=sig0, col="red", lty=2, lwd=1) 
                mtext(paste(siglines[j]*100,"%",sep=""), side=4, line=0, at=ceiling(sig0), col="red")
            }
        } ## end for
    } ## end logical
}

## this is a hack method to strip depreciated hidden parameters.
## This suggestion is from Gavin Simpson's blog and he traces it
## to a suggestion from Brian Ripley
## see: http://ucfagls.wordpress.com/2011/07/23/
## local hidden plotting routine to strip depreciated parameter
## dT from arguments lost
.lplotSpec <- function(x, ..., dT) {
    ## should call plot.spec prior to 3.1 dev.
    plot(x, ...)
}

## this is currently used in F--test plots.
## modified Feb 2015 to compile for cran
.lplotDefault <- function(x, y, ..., dT) {
    ## should call plot default 
    plot(x, y, ...)
}

## utilities functions added for plot.mtm.coh
## fortran versions exist and could be cleaned up and implemented...
## nhi is nfreqs
## c
## c x       Input Data
## c xv      Variance estimate of X
## c xp      Plotting Array, Returned
## c         xp(.,1) mean - 1 standard deviation
## c         xp(.,2) smoothed x
## c         xp(.,3) mean + 1 standard deviation
## c nehl    Number each half for smoothing limits
## -- for smoothing variance
## c nehc    Number Each Half for Center
## -- for smoothing data
## c slo     Low end symmetry - "even","odd", "extend"
## c shi   High end symmetry
.lftr3p <- function(x, xv, nhi, nehl, nehc, slo, shi) {
    ip <- 2

    xp <- matrix(NA, nhi, 3)

    xp[,3] <- .llftr7(xv, nhi, "hi", "even", "extend", nehl, ip)
    xp[,2] <- .llftr7(x, nhi, "-", slo, shi, nehc, ip)

    xp[,1] <- xp[,2] - sqrt(xp[,3])
    xp[,3] <- xp[,2] + sqrt(xp[,3])

    return(xp)
}


##modified djt jk
##Programmers note. currently uses jkcoh7 which containes segment averaging
##Segment averaging can be looked into and implemented
##Segment averaging can likely be vectorized or should it be
##implemented in C?
## Dave sets ip=2 in the file ts/lftr3p.f
.llftr7 <- function(x,nhi,lohi,slo,shi,neh,ip) {

    nlo <- 1
    y <- array(NA, nhi)
    z <- array(NA, nhi+2*neh)

    zNlo <- nlo + neh
    zNhi <- nhi + neh
    zHi <- zNhi + neh

    ##  Generate Weights
    fw <- as.double(neh + 1)
    wt <- ((1-((-neh:neh)/fw))*(1+((-neh:neh)/fw)))**ip
    cwt <- sum(wt)
    wt <-  wt/cwt
    
    ## Move data to working (z) array, extend ends,default
    innerSeq <- zNlo:zNhi
    z[innerSeq] <- x
    lowerSeq <- 1:neh
    z[lowerSeq] <- x[1]
    upperSeq <- (zNhi+1):zHi
    z[upperSeq] <- x[nhi]
    
    ## Low End
    
    if(tolower(slo) == "even") {
        z[rev(lowerSeq)] <-  z[(zNlo+1):(zNlo+neh)]
    }

    
    if(tolower(slo) == "odd") {
        ## not tested...
        z[rev(lowerSeq)] <-  2.0*z[zNlo] - z[(zNlo+1):(zNlo+neh)]
    }

    ## High End
    ## Programmer's note: The numbers map to numbers on do loops
    ## in Thomson's fortran code.
    if(tolower(shi) == "even") {
        for( j in  1:neh) { ## 850
            z[zNhi+j] <- z[zNhi-j]
        } ## 850 
    }

    if(tolower(shi) == "odd") {
        for( j  in 1:neh) { ## 860
            z[zNhi + j] =  2.*z[zNhi] - z[zNhi-j]
        } ## 860
    }
    
    ## High Limit, supress local Minima
    if(tolower(lohi) ==  "hi") {
        for( n in  (nlo+1):(nhi-1) ) { ## 1400
            if( (x[n] < x[n-1]) &&  (x[n]  < x[n+1]) ) {
                zNoff <- n +neh
                z[zNoff] <- (x[n-1]+x[n+1])/2.0
            }
        } ## 1400
    }
    ## Low Limit, supress local Maxima
    if(tolower(lohi) == "lo") {
        for( n in  (nlo+1):(nhi-1) ) {
            if( (x[n] > x[n-1] ) &&  (x[n] > x[n+1]) ) {
                zNoff <- n +neh
                z[n] = (x[n-1]+x[n+1])/2.0
            }
        } ## 1500
    }

    zOffSetSeq <- 1:(2*neh+1)
    for (n in nlo:nhi) { ## 2000 
        y[n] <- sum(wt*z[zOffSetSeq])
        zOffSetSeq <- zOffSetSeq +1
    } ## 2000
    return(y)
    ##checks out on first test
    
}


### functions to convert the coherence to the transformed coherence
## normalizing constant 2k-2
.trnrm <- function(k) sqrt(2*k-2)

## These formulae are based on:
## "Jackknifed error estimates for spectra, coherences,
## and transfer functions"
## by Thomson, DJ and Chave, AD
## Advances in Spectrum Estimation
##

## coherence to quantiles of the CDF
.C2toF <-  function(xx, trnrm_) {
    return( trnrm_*log((1.0+sqrt(xx))/(1.0-sqrt(xx)))/2.0 )
}

##quantiles to MSC
.FtoMSC <- function(ff, trnrm_) tanh(ff/trnrm_)**2


## odinlibs nplot function....
## c
## c       Cumulative probability points in 1-2-5 sequence
## c ndata Number data points; output approximately from 1/ndata
## c       to 1 - 1/ndata
## c nmax  Maximum number of Outputs = Dimension of out,cout,Qnorm
## c       nmax approx > 6* log10(ndata)
## c
## c out   Cumulative distribution
## c cout  Character*8 version of CDF
## c Qnorm Quantiles of Standard Normal at CDF
## c nout  Number of output points
## c
## using jkcoh defaults...
.paxpt7 <- function(ndata=2000, nmax=40) {
    ndec <-  round(log10(max(11,ndata)));
    nout = 6*ndec -1;
    if(nout > nmax) {
        return;
    }
    n <- 0;
    out <- array(NA, nout)
    for(m in seq(-ndec, -1, 1)) {
        for(k in c(1,2,5)) {
            n <- n +1;
            v <- as.double(k*10**m);
            out[n] <- v;
            out[nout +1 -n] <- as.double(1.0 - v);
        }
    }
    return(list(out=out, Qnorm=qnorm(out),nout=nout));
}

.cdfToMSqCoh <- function(cdf, k) {
    fnavm <- as.double(k-1);
    return(1.0  - (1.0 - cdf)**(1.0/fnavm));
}

## used in coherence plot--added May, 2013.
.mscToCDFquantiles <- function(msc, k) {
    1 - (1-msc)^(k-1)
}
## end utilities added mainly for plot.mtm.coh

Try the multitaper package in your browser

Any scripts or data that you put into this service are public.

multitaper documentation built on Nov. 17, 2017, 4:02 a.m.