The RobPer-package

Description

Calculates periodograms based on (robustly) fitting periodic functions to light curves and other irregulary observed time series and detects high periodogram bars.

Details

Light curves occur in astroparticle physics and are irregularly sampled times series (t[i], y[i]), i=1,…,n, or (t[i], y[i], s[i]), i=1,…,n, consisting of unequally spaced observation times t[1],…,t[n], observed values y[1],…,y[n] and possibly measurement accuracies s[1],…,s[n]. The pattern of the observation times t[i] may be periodic with sampling period ps. The observed values y[i] may possibly contain a periodic fluctuation yf[i] with fluctuation period pf. One is interested in finding pf. The measurement accuracies s[i] give information about how precise the y[i] were measured. They can be interpreted as estimates for the standard deviations of the observed values. For more details see Thieler et al. (2013) or Thieler, Fried and Rathjens (2016).

This package includes three main functions: RobPer calculates the periodogram, possibly taking into account measurement accuracies. With betaCvMfit, outlying periodogram bars (indicating a period) can be detected. This function bases on robustly fitting a distribution using Cramér-von-Mises (CvM) distance minimization (see also Clarke, McKinnon and Riley 2012). The function tsgen can be used to generate artificial light curves. For more details about the implementation see Thieler, Fried and Rathjens (2016).

A preliminary version of this package is used in Thieler et al. (2013). The FastS-function and the FastTau-function presented here are slightly changed versions of R-Code published in Salibian-Barrera and Yohai (2006) and Salibian-Barrera, Willems and Zamar (2008).

The financial support of the DFG (SFB 876 "Providing Information by Resource-Constrained Data Analysis", project C3, and GK 1032 "Statistische Modellbildung") is gratefully acknowledged. We thank the ITMC at TU Dortmund University for providing computer resources on LiDO.

Author(s)

Anita M. Thieler, Jonathan Rathjens and Roland Fried, with contributions from Brenton R. Clarke (see betaCvMfit), Matias Salibian-Barrera, Gert Willems and Victor Yohai (see FastS and FastTau) and Uwe Ligges (see TK95).

Maintainer: Jonathan Rathjens <jonathan.rathjens@tu-dortmund.de>

References

Clarke, B. R., McKinnon, P. L. and Riley, G. (2012): A Fast Robust Method for Fitting Gamma Distributions. Statistical Papers, 53 (4), 1001-1014

Salibian-Barrera, M. and Yohai, V. (2006): A Fast Algorithm for S-Regression Estimates. Journal of Computational and Graphical Statistics, 15 (2), 414-427

Salibian-Barrera, M., Willems, G. and Zamar, R. (2008): The Fast-tau Estimator for Regression. Journal of Computational and Graphical Statistics, 17 (3), 659-682

Thieler, A. M., Backes, M., Fried, R. and Rhode, W. (2013): Periodicity Detection in Irregularly Sampled Light Curves by Robust Regression and Outlier Detection. Statistical Analysis and Data Mining, 6 (1), 73-89

Thieler, A. M., Fried, R. and Rathjens, J. (2016): RobPer: An R Package to Calculate Periodograms for Light Curves Based on Robust Regression. Journal of Statistical Software, 69 (9), 1-36, <doi:10.18637/jss.v069.i09>

Examples

# Generate a disturbed light curve:
set.seed(22)
lightcurve <- tsgen(ttype="sine",ytype="peak" , pf=7, redpart=0.1, s.outlier.fraction=0.1,
    interval=TRUE, npoints=200, ncycles=25, ps=20, SNR=3, alpha=0)

# Plotting the light curve (vertical bars show measurement accuracies)
plot(lightcurve[,1], lightcurve[,2], pch=16, cex=0.5, xlab="t", ylab="y", 
    main="a Light Curve")
rect(lightcurve[,1], lightcurve[,2]+lightcurve[,3], lightcurve[,1], 
    lightcurve[,2]-lightcurve[,3])

# The lightcurve has a period of 7:
plot(lightcurve[,1]%%7, lightcurve[,2], pch=16, cex=0.5, xlab="t", ylab="y",
    main="Phase Diagram of a Light Curve")
rect(lightcurve[,1]%%7, lightcurve[,2]+lightcurve[,3], lightcurve[,1]%%7, 
    lightcurve[,2]-lightcurve[,3])

# Calculate a periodogram of a light curve:
PP <- RobPer(lightcurve, model="splines", regression="huber", weighting=FALSE, 
    var1=FALSE, periods=1:50)

# Searching for extremely high periodogram bars:
betavalues <- betaCvMfit(PP)
crit.val <- qbeta((0.95)^(1/50),shape1=betavalues[1], shape2=betavalues[2])

hist(PP, breaks=20, freq=FALSE, ylim=c(0,100), xlim=c(0,0.08), col=8, main ="")
betafun <- function(x) dbeta(x, shape1=betavalues[1], shape2=betavalues[2])
curve(betafun, add=TRUE, lwd=2)
abline(v=crit.val, lwd=2)

# alternatives for fitting beta distributions:
# method of moments:
par.mom <- betaCvMfit(PP, rob=FALSE, CvM=FALSE)
myf.mom <- function(x) dbeta(x, shape1=par.mom[1], shape2=par.mom[2])
curve(myf.mom, add=TRUE, lwd=2, col="red")
crit.mom <- qbeta((0.95)^(1/50),shape1=par.mom[1], shape2=par.mom[2])
abline(v=crit.mom, lwd=2, col="red")

# robust method of moments
par.rob <- betaCvMfit(PP, rob=TRUE, CvM=FALSE)
myf.rob <- function(x) dbeta(x, shape1=par.rob[1], shape2=par.rob[2])
curve(myf.rob, add=TRUE, lwd=2, col="blue")
crit.rob <- qbeta((0.95)^(1/50),shape1=par.rob[1], shape2=par.rob[2])
abline(v=crit.rob, lwd=2, col="blue")

legend("topright", fill=c("black","red","blue"), 
    legend=c("CvM", "moments", "robust moments"), bg="white")
box()

# Detect fluctuation period:
plot(1:50, PP, xlab="Trial Period", ylab="Periodogram", type="l", 
    main="Periodogram fitting periodic splines using M-regression (Huber function)")
abline(h=crit.val, lwd=2)
text(c(7,14), PP[c(7,14)], c(7,14), adj=1, pos=4)
axis(1, at=7, labels=expression(p[f]==7))

# Comparison with non-robust periodogram
# (see package vignette, section 5.1 for further graphical analysis)
PP2 <- RobPer(lightcurve, model="splines", regression="L2",
    weighting=FALSE, var1=FALSE, periods=1:50)
betavalues2 <- betaCvMfit(PP2)
crit.val2   <- qbeta((0.95)^(1/50),shape1=betavalues2[1], shape2=betavalues2[2])

plot(1:50, PP2, xlab="Trial Period", ylab="Periodogram", type="l", 
    main="Periodogram fitting periodic splines using L2-regression")
abline(h=crit.val2, lwd=2)

Results

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(RobPer)
Loading required package: robustbase
Loading required package: quantreg
Loading required package: SparseM

Attaching package: 'SparseM'

The following object is masked from 'package:base':

    backsolve

Loading required package: splines
Loading required package: BB
Loading required package: rgenoud
##  rgenoud (Version 5.7-12.4, Build Date: 2015-07-19)
##  See http://sekhon.berkeley.edu/rgenoud for additional documentation.
##  Please cite software as:
##   Walter Mebane, Jr. and Jasjeet S. Sekhon. 2011.
##   ``Genetic Optimization Using Derivatives: The rgenoud package for R.''
##   Journal of Statistical Software, 42(11): 1-26. 
##

> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/RobPer/RobPer-package.rd_%03d_medium.png", width=480, height=480)
> ### Name: RobPer-package
> ### Title: The RobPer-package
> ### Aliases: RobPer-package
> 
> ### ** Examples
> 
> # Generate a disturbed light curve:
> set.seed(22)
> lightcurve <- tsgen(ttype="sine",ytype="peak" , pf=7, redpart=0.1, s.outlier.fraction=0.1,
+     interval=TRUE, npoints=200, ncycles=25, ps=20, SNR=3, alpha=0)
  Successful convergence.
> 
> # Plotting the light curve (vertical bars show measurement accuracies)
> plot(lightcurve[,1], lightcurve[,2], pch=16, cex=0.5, xlab="t", ylab="y", 
+     main="a Light Curve")
> rect(lightcurve[,1], lightcurve[,2]+lightcurve[,3], lightcurve[,1], 
+     lightcurve[,2]-lightcurve[,3])
> 
> # The lightcurve has a period of 7:
> plot(lightcurve[,1]%%7, lightcurve[,2], pch=16, cex=0.5, xlab="t", ylab="y",
+     main="Phase Diagram of a Light Curve")
> rect(lightcurve[,1]%%7, lightcurve[,2]+lightcurve[,3], lightcurve[,1]%%7, 
+     lightcurve[,2]-lightcurve[,3])
> 
> # Calculate a periodogram of a light curve:
> PP <- RobPer(lightcurve, model="splines", regression="huber", weighting=FALSE, 
+     var1=FALSE, periods=1:50)
> 
> # Searching for extremely high periodogram bars:
> betavalues <- betaCvMfit(PP)
> crit.val <- qbeta((0.95)^(1/50),shape1=betavalues[1], shape2=betavalues[2])
> 
> hist(PP, breaks=20, freq=FALSE, ylim=c(0,100), xlim=c(0,0.08), col=8, main ="")
> betafun <- function(x) dbeta(x, shape1=betavalues[1], shape2=betavalues[2])
> curve(betafun, add=TRUE, lwd=2)
> abline(v=crit.val, lwd=2)
> 
> # alternatives for fitting beta distributions:
> # method of moments:
> par.mom <- betaCvMfit(PP, rob=FALSE, CvM=FALSE)
> myf.mom <- function(x) dbeta(x, shape1=par.mom[1], shape2=par.mom[2])
> curve(myf.mom, add=TRUE, lwd=2, col="red")
> crit.mom <- qbeta((0.95)^(1/50),shape1=par.mom[1], shape2=par.mom[2])
> abline(v=crit.mom, lwd=2, col="red")
> 
> # robust method of moments
> par.rob <- betaCvMfit(PP, rob=TRUE, CvM=FALSE)
> myf.rob <- function(x) dbeta(x, shape1=par.rob[1], shape2=par.rob[2])
> curve(myf.rob, add=TRUE, lwd=2, col="blue")
> crit.rob <- qbeta((0.95)^(1/50),shape1=par.rob[1], shape2=par.rob[2])
> abline(v=crit.rob, lwd=2, col="blue")
> 
> legend("topright", fill=c("black","red","blue"), 
+     legend=c("CvM", "moments", "robust moments"), bg="white")
> box()
> 
> # Detect fluctuation period:
> plot(1:50, PP, xlab="Trial Period", ylab="Periodogram", type="l", 
+     main="Periodogram fitting periodic splines using M-regression (Huber function)")
> abline(h=crit.val, lwd=2)
> text(c(7,14), PP[c(7,14)], c(7,14), adj=1, pos=4)
> axis(1, at=7, labels=expression(p[f]==7))
> 
> # Comparison with non-robust periodogram
> # (see package vignette, section 5.1 for further graphical analysis)
> PP2 <- RobPer(lightcurve, model="splines", regression="L2",
+     weighting=FALSE, var1=FALSE, periods=1:50)
> betavalues2 <- betaCvMfit(PP2)
> crit.val2   <- qbeta((0.95)^(1/50),shape1=betavalues2[1], shape2=betavalues2[2])
> 
> plot(1:50, PP2, xlab="Trial Period", ylab="Periodogram", type="l", 
+     main="Periodogram fitting periodic splines using L2-regression")
> abline(h=crit.val2, lwd=2)
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>