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,\ldots,n} or (t_i, y_i, s_i)_{i=1,\ldots,n} consisting of unequally spaced observation times t_1, \ldots, t_n, observed values y_1, \ldots, y_n and possibly measurement accuracies s_1, \ldots, s_n. The pattern of the observation times t_i may be periodic with sampling period p_s. The observed values y_i may possibly contain a periodic fluctuation y_{f;i} with fluctuation period p_f. One is interested in finding p_f. 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)

S-Regression using the Fast-S-Algorithm

Description

Performs S-Regression using the Fast-S-Algorithm.

Usage

FastS(x, y, Scontrol=list(int = FALSE, N = 100, kk = 2, tt = 5, b= .5,
 cc = 1.547, seed=NULL), beta_gamma)

Arguments

Details

The Fast-S-Algorithm to efficiently perform S-Regression was published by Salibian-Barrera and Yohai (2006). It bases on starting with a set of N parameter candidates, locally optimizing them, but only with kk iterations, optimizing the tt best candidates to convergence and then choosing the best parameter candidate. The rho-function used is the biweight function with tuning parameter cc, the value b is set to the expected value of the rho-function applied to the residuals. The default cc=1.547 and b=.5 is chosen following Rousseeuw and Yohai (1984) to obtain an approximative breakdown point of 0.5. When setting int to TRUE, this adds an intercept column to the design matrix. For more details see Salibian-Barrera and Yohai (2006) or Thieler, Fried and Rathjens (2016).

The R-function FastS used in RobPer is a slightly changed version of the R-code published in Salibian-Barrera and Yohai (2006). It was changed in order to work more efficiently, especially when fitting step functions, and to specify one parameter candidate in advance. For details see Thieler, Fried and Rathjens (2016).

Value

Author(s)

Matias Salibian-Barrera and Victor Yohai, modified by Anita M. Thieler

References

Rousseeuw, P. J. and Yohai, V. J. (1984): Robust Regression by Means of S-estimators. In Franke, J., Härdle, W. und Martin, D. (eds.): Robust and Nonlinear Time Series Analysis. Berlin New York: Springer, Lecture Notes in Statistics No. 26, 256-272

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

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>

See Also

Applied in RobPer. See FastTau for example.

Tau-Regression using the Fast-tau-Algorithm

Description

Performs tau-Regression using the Fast-tau-Algorithm.

Usage

FastTau(x, y, taucontrol = list(N = 500, kk = 2, tt = 5, rr = 2, approximate = 0),
 beta_gamma)

Arguments

Details

The Fast-tau-Algorithm to efficiently perform tau-Regression was published by Salibian-Barrera, Willems and Zamar (2008). It bases on starting with a set of N parameter candidates, locally optimizing them using kk iterations, then optimizing the tt best candidates to convergence and finally choosing the best parameter candidate. Since calculation of the objective value is computationally expensive, it is possible to approximate it with rr iteration steps when choosing approximate=TRUE. For more details see Salibian-Barrera, Willems and Zamar (2008).

The R-function FastTau used in RobPer is a slightly changed version of the R-code published in Salibian-Barrera, Willems and Zamar (2008). It was changed in order to work more efficiently, especially when fitting step functions, and to specify one parameter candidate in advance. For details see Thieler, Fried and Rathjens (2016).

Value

Author(s)

Matias Salibian-Barrera and Gert Willems, modified by Anita M. Thieler

References

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., 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>

See Also

Applied in RobPer.

Examples

set.seed(22)
# Generate a disturbed light curve
lightcurve <- tsgen(ttype="unif",ytype="sine" , pf=7, redpart=0.1, interval=TRUE,
    npoints=100, ncycles=10, ps=7, SNR=4, alpha=0)
tt <- lightcurve[,1]
y  <- lightcurve[,2]
s  <- rep(1,100)  # unweighted regression

plot(tt, y, type="l", main="Fitting a sine to a disturbed lightcurve")

# Fit the true model (a sine of period 7)... designmatrix:
X <- Xgen(tt, n=100, s, pp=7, design="sine")
# Robust tau-fit:
beta_FastTau <- FastTau(X, y)$beta
# Robust S-fit:
beta_FastS <- FastS(X, y)$beta
# Least squares fit:
beta_lm <- lm(y~0+X)$coeff

# Plot:
sin7_fun <- function(t, beta) beta[1]+ beta[2]*sin(t*2*pi/7)+ beta[3]*cos(t*2*pi/7)
sin_FastTau <- function(t) sin7_fun(t, beta_FastTau)
sin_FastS <- function(t) sin7_fun(t, beta_FastS)
sin_lm <- function(t) sin7_fun(t, beta_lm)
curve(sin_FastTau, col="green", add=TRUE)
curve(sin_FastS, col="blue", add=TRUE, lty=2)
curve(sin_lm, col="red", add=TRUE)

legend("topleft", fill=c("black", "red", "green", "blue"),
    legend=c("Light Curve (disturbed)", "Least Squares Fit", "FastTau Fit", "FastS Fit"))

Data: Light curve from Mrk 421

Description

Gamma ray light curve from Markarian 421.

Usage

Mrk421

Format

A data frame of three variables, with a time series of length 655 appropriate to RobPer.

Details

The data in Mrk421 and Mrk501 have been collected from various original sources, combined, and published by Tluczykont et al. (2010) of the Deutsches Elektronen-Synchrotron.

Their sources are data from the experiments:

Whipple (Kerrick et al. 1995; Schubnell et al. 1996; Buckley et al. 1996; Maraschi et al. 1999)

HEGRA (Aharonian et al. 1999a, 1999b; Krawczynski et al. 2001; Aharonian et al. 2001, 2002, 2003, 2004; Kestel 2003)

CAT (Piron 2000; Piron et al. 2001)

HESS (Aharonian et al. 2005, 2006)

MAGIC (Albert et al. 2008; Donnarumma et al. 2009)

VERITAS (Rebillot et al. 2006; Donnarumma et al. 2009)

Note

See Vignette Section 5.3 for example.

Source

Data kindly provided by the Deutsches Elektronen-Synchrotron, Gamma Astronomy group (see Details).

References

Tluczykont, M., Bernardini, E., Satalecka, K., Clavero, R., Shayduk, M. and Kalekin, O. (2010): Long-term lightcurves from combined united very high energy gamma-ray data. Astronomy & Astrophysics, 524, A48

Their data sources:

Aharonian, F., Akhperjanian, A., Barrio, J., et al. (1999a): The temporal characteristics of the TeV gamma-radiation from Mkn 501 in 1997: I. Data from the stereoscopic imaging atmospheric Cherenkov telescope system of HEGRA. Astronomy & Astrophysics, 342(1), 69

Aharonian, F., Akhperjanian, A., Barrio, J., et al. (1999b): The temporal characteristics of the TeV gamma-emission from Mkn 501 in 1997: II. Results from HEGRA CT1 and CT2. Astronomy & Astrophysics, 349(1), 29

Aharonian, F., Akhperjanian, A., Barrio, J., et al. (2001): The TeV Energy Spectrum of Markarian 501 Measured with the Stereoscopic Telescope System of HEGRA during 1998 and 1999. The Astrophysical Journal, 546(2), 898

Aharonian, F., Akhperjanian, A., Beilicke, M., et al. (2002): Variations of the TeV energy spectrum at different flux levels of Mkn 421 observed with the HEGRA system of Cherenkov telescopes. Astronomy & Astrophysics, 393(1), 89

Aharonian, F., Akhperjanian, A., Beilicke, M., et al. (2003): TeV gamma-ray light curve and energy spectrum of Mkn 421 during its 2001 flare as measured with HEGRA CT1. Astronomy & Astrophysics, 410(3), 813

Aharonian, F., Akhperjanian, A., Beilicke, M., et al. (2004): The Crab Nebula and Pulsar between 500 GeV and 80 TeV: Observations with the HEGRA Stereoscopic Air Cerenkov Telescopes. The Astrophysical Journal, 614(2), 897

Aharonian, F., Akhperjanian, A., Aye, K., et al. (2005): Observations of Mkn 421 in 2004 with HESS at large zenith angles. Astronomy & Astrophysics, 437(1), 95

Aharonian, F., Akhperjanian, A., Bazer-Bachi, A. R., et al. (2006): Observations of the Crab nebula with HESS. Astronomy & Astrophysics, 457(3), 899

Albert, J., Aliu, E., Anderhub, H., et al. (2008): VHE gamma-Ray Observation of the Crab Nebula and its Pulsar with the MAGIC Telescope. The Astrophysical Journal, 674(2), 1037

Buckley, J. H., Akerlof, C. W., Biller, S., et al. (1996): Gamma-Ray Variability of the BL Lacertae Object Markarian 421. The Astrophysical Journal, 472, L9

Donnarumma, I., Vittorini, V., Vercellone, S., et al. (2009): The June 2008 Flare of Markarian 421 from Optical to TeV Energies. The Astrophysical Journal, 691, L13

Kerrick, A. D., Akerlof, C. W., Biller, S. D., et al. (1995): Outburst of TeV photons from Markarian 421. The Astrophysical Journal, 438, L59

Kestel, M. (2003): TeV gamma-Flux and Spectrum of Markarian 421 in 1999/2000 with Hegra CT1 using refined Analysis Methods. Ph.D. Thesis, Technische Universit\"at M\"unchen

Krawczynski, H., Sambruna, R., Kohnle, A., et al. (2001): Simultaneous X-Ray and TeV Gamma-Ray Observation of the TeV Blazar Markarian 421 during 2000 February and May. The Astrophysical Journal, 559(1), 187

Maraschi, L., Fossati, G., Tavecchio, F., et al. (1999): Simultaneous X-Ray and TeV Observations of a Rapid Flare from Markarian 421. The Astrophysical Journal, 526, L81

Piron, F. (2000): \'Etude des propri\'et\'es spectrales et de la variabilit\'e de l'\'emission Gamma sup\'erieure \'e 250 GeV des blazars observ\'es entre Octobre 1996 et F\'evrier 2000 par le t\'elescope \'e effet Tcherenkov athmosph\'erique C.A.T. Ph.D. Thesis, Universit\'e de Paris-Sud

Piron, F., Djannati-Atai, A., Punch, M., et al. (2001): Temporal and spectral gamma-ray properties of Mkn 421 above 250 GeV from CAT observations between 1996 and 2000. Astronomy & Astrophysics, 374(3), 895

Rebillot, P. F., Badran, H. M., Blaylock, G., et al. (2006): Multiwavelength Observations of the Blazar Markarian 421 in 2002 December and 2003 January. The Astrophysical Journal, 641(2), 740

Schubnell, M. S., Akerlof, C. W., Biller, S., et al. (1996): Very High Energy Gamma-Ray Emission from the Blazar Markarian 421. The Astrophysical Journal, 460, 644

Data: Light curve from Mrk 501

Description

Gamma ray light curve from Markarian 501.

Usage

Mrk501

Format

A data frame of three variables, with a time series of length 210 appropriate to RobPer.

Details

The data in Mrk421 and Mrk501 have been collected from various original sources, combined, and published by Tluczykont et al. (2010) of the Deutsches Elektronen-Synchrotron.

Their sources are data from the experiments:

Whipple (Kerrick et al. 1995; Schubnell et al. 1996; Buckley et al. 1996; Maraschi et al. 1999)

HEGRA (Aharonian et al. 1999a, 1999b; Krawczynski et al. 2001; Aharonian et al. 2001, 2002, 2003, 2004; Kestel 2003)

CAT (Piron 2000; Piron et al. 2001)

HESS (Aharonian et al. 2005, 2006)

MAGIC (Albert et al. 2008; Donnarumma et al. 2009)

VERITAS (Rebillot et al. 2006; Donnarumma et al. 2009)

Note

See Vignette Section 5.3 for example.

Source

Data kindly provided by the Deutsches Elektronen-Synchrotron, Gamma Astronomy group (see Details).

References

Tluczykont, M., Bernardini, E., Satalecka, K., Clavero, R., Shayduk, M. and Kalekin, O. (2010): Long-term lightcurves from combined united very high energy gamma-ray data. Astronomy & Astrophysics, 524, A48

Their data sources:

Aharonian, F., Akhperjanian, A., Barrio, J., et al. (1999a): The temporal characteristics of the TeV gamma-radiation from Mkn 501 in 1997: I. Data from the stereoscopic imaging atmospheric Cherenkov telescope system of HEGRA. Astronomy & Astrophysics, 342(1), 69

Aharonian, F., Akhperjanian, A., Barrio, J., et al. (1999b): The temporal characteristics of the TeV gamma-emission from Mkn 501 in 1997: II. Results from HEGRA CT1 and CT2. Astronomy & Astrophysics, 349(1), 29

Aharonian, F., Akhperjanian, A., Barrio, J., et al. (2001): The TeV Energy Spectrum of Markarian 501 Measured with the Stereoscopic Telescope System of HEGRA during 1998 and 1999. The Astrophysical Journal, 546(2), 898

Aharonian, F., Akhperjanian, A., Beilicke, M., et al. (2002): Variations of the TeV energy spectrum at different flux levels of Mkn 421 observed with the HEGRA system of Cherenkov telescopes. Astronomy & Astrophysics, 393(1), 89

Aharonian, F., Akhperjanian, A., Beilicke, M., et al. (2003): TeV gamma-ray light curve and energy spectrum of Mkn 421 during its 2001 flare as measured with HEGRA CT1. Astronomy & Astrophysics, 410(3), 813

Aharonian, F., Akhperjanian, A., Beilicke, M., et al. (2004): The Crab Nebula and Pulsar between 500 GeV and 80 TeV: Observations with the HEGRA Stereoscopic Air Cerenkov Telescopes. The Astrophysical Journal, 614(2), 897

Aharonian, F., Akhperjanian, A., Aye, K., et al. (2005): Observations of Mkn 421 in 2004 with HESS at large zenith angles. Astronomy & Astrophysics, 437(1), 95

Aharonian, F., Akhperjanian, A., Bazer-Bachi, A. R., et al. (2006): Observations of the Crab nebula with HESS. Astronomy & Astrophysics, 457(3), 899

Albert, J., Aliu, E., Anderhub, H., et al. (2008): VHE gamma-Ray Observation of the Crab Nebula and its Pulsar with the MAGIC Telescope. The Astrophysical Journal, 674(2), 1037

Buckley, J. H., Akerlof, C. W., Biller, S., et al. (1996): Gamma-Ray Variability of the BL Lacertae Object Markarian 421. The Astrophysical Journal, 472, L9

Donnarumma, I., Vittorini, V., Vercellone, S., et al. (2009): The June 2008 Flare of Markarian 421 from Optical to TeV Energies. The Astrophysical Journal, 691, L13

Kerrick, A. D., Akerlof, C. W., Biller, S. D., et al. (1995): Outburst of TeV photons from Markarian 421. The Astrophysical Journal, 438, L59

Kestel, M. (2003): TeV gamma-Flux and Spectrum of Markarian 421 in 1999/2000 with Hegra CT1 using refined Analysis Methods. Ph.D. Thesis, Technische Universit\"at M\"unchen

Krawczynski, H., Sambruna, R., Kohnle, A., et al. (2001): Simultaneous X-Ray and TeV Gamma-Ray Observation of the TeV Blazar Markarian 421 during 2000 February and May. The Astrophysical Journal, 559(1), 187

Maraschi, L., Fossati, G., Tavecchio, F., et al. (1999): Simultaneous X-Ray and TeV Observations of a Rapid Flare from Markarian 421. The Astrophysical Journal, 526, L81

Piron, F. (2000): \'Etude des propri\'et\'es spectrales et de la variabilit\'e de l'\'emission Gamma sup\'erieure \'e 250 GeV des blazars observ\'es entre Octobre 1996 et F\'evrier 2000 par le t\'elescope \'e effet Tcherenkov athmosph\'erique C.A.T. Ph.D. Thesis, Universit\'e de Paris-Sud

Piron, F., Djannati-Atai, A., Punch, M., et al. (2001): Temporal and spectral gamma-ray properties of Mkn 421 above 250 GeV from CAT observations between 1996 and 2000. Astronomy & Astrophysics, 374(3), 895

Rebillot, P. F., Badran, H. M., Blaylock, G., et al. (2006): Multiwavelength Observations of the Blazar Markarian 421 in 2002 December and 2003 January. The Astrophysical Journal, 641(2), 740

Schubnell, M. S., Akerlof, C. W., Biller, S., et al. (1996): Very High Energy Gamma-Ray Emission from the Blazar Markarian 421. The Astrophysical Journal, 460, 644

Periodogram based on (robustly) fitting a periodic function to a light curve

Description

Calculates a periodogram by fitting a periodic function to a light curve, using a possibly robust regression technique and possibly taking into account measurement accuracies. See RobPer-package for more information about light curves. For a lot of more details see Thieler, Fried and Rathjens (2016) and Thieler et al. (2013).

Usage

RobPer(ts, weighting, periods, regression, model, steps = 10, tol = 1e-03,
 var1 = weighting, genoudcontrol = list(pop.size = 50, max.generations = 50,
 wait.generations = 5), LTSopt =TRUE, 
 taucontrol = list(N = 100, kk = 2, tt = 5, rr = 2, approximate = FALSE),
 Scontrol=list(N = ifelse(weighting,200,50), kk = 2, tt = 5, b=.5, cc = 1.547,
 seed = NULL) )

Arguments

Details

For each trial period, a periodic function (defined by model) is fitted to the light curve using regression technique regression. The periodogram bar is the coefficient of determination. In case of model="2step", two different step functions with opposed jumping times are fitted separately and the periodogram bar is the mean of both coefficients of determination. For a lot of more details see Thieler, Fried and Rathjens (2016) and Thieler et al. (2013).

Value

numeric vector: Periodogram bars related to the trial periods.

Note

Performing weighting = FALSE, regression="L2", model="sine" on a equidistantly sampled time series is equivalent to calculating the standard periodogram of Fourier analysis, see Example.

Performing regression="L2", model="sine" is equivalent to calculating a Generalized Lomb-Scargle periodogram (see Zechmeister and Kürster 2009).

Performing regression="L2", model="step" is equivalent to calculating an Epoch Folding (Leahy et al. 1983) or Anaysis of Variance (Schwarzenberg-Czerny 1989) periodogram.

Performing regression="L2", model="2step" is equivalent to calculating a Phase Dispersion Minimization periodogram (Stellingwerf 1978).

A former version of this function is used in Thieler et al. (2013). For more equivalences see there.

Author(s)

Anita M. Thieler, Jonathan Rathjens and Roland Fried

References

Leahy, D. A., Darbro, W., Elsner, R. F., Weisskopf, M. C., Kahn, S., Sutherland, P. G. and Grindlay, J. E. (1983): On Searches for Pulsed Emission with Application to Four Globular Cluster X-ray Sources-NGC 1851, 6441, 6624, and 6712. The Astrophysical Journal, 266 (1), 160-170

Mebane Jr., W. R. and Sekhon, J. S. (2011): Genetic Optimization Using Derivatives: The rgenoud Package for R. Journal of Statistical Software, 42 (11), 1-26

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

Stellingwerf, R. F. (1978): Period Determination Using Phase Dispersion Minimization. The Astrophysical Journal, 224, 953-960

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>

Zechmeister, M. and Kürster, M. (2009): The Generalised Lomb-Scargle Periodogram. A New Formalism for the Floating-Mean and Keplerian Periodograms. Astronomy and Astrophysics, 496 (2), 577-584

See Also

Applies FastS and FastTau, Xgen, examples in RobPer-package and TK95_uneq.

Examples

# For more examples see RobPer-package and TK95_uneq!

# Example to show the equivalence between the periodogram from Fourier analysis
# and the Lomb-Scargle periodogram in case of equidistant sampling and equal weighting:
set.seed(7)
n <- 120
# equidistant time series:
zr <- tsgen(ttype="equi", ytype="const", pf=1, redpart= 0, s.outlier.fraction=0.2, 
    interval=FALSE, npoints=n, ncycles=n, ps=1, SNR=1, alpha=1.5)
# periodogram of Fourier analysis
PP_konv <- spec.pgram(zr[,2], taper = 0,pad = 0, fast = FALSE, demean = TRUE,
    detrend = TRUE, plot = TRUE)
# Lomb-Scargle periodogram - Note: Due to the regression ansatz,
# RobPer is not able to compute period 2 in this case.
PP_new <- RobPer(ts=zr, weighting=FALSE, periods=1/PP_konv$freq,
    regression="L2", model="sine")
plot(PP_konv$freq, PP_konv$spec, ylab="periodogram", xlab="frequency",
    main="Comparison of RobPer(...regression='LS', model='sine') and spec.pgram")
points(PP_konv$freq, PP_new*var(zr[,2])*n/2, type="l")
legend("top",lty=c(1,0), pch=c(-5,1), legend=c("RobPer*var(y)*n/2", "spec.pgram"))
# Due to different ways of computation, the scaled periodograms are not exactly
# identical, but show very similar behavior.

Power law noise generator

Description

Generates an equidistant time series of power law noise according to Timmer and König (1995).

Usage

TK95(N = 1000, alpha = 1.5)

Arguments

Value

numeric vector: The generated time series.

Note

This function is used in Thieler et al. (2013). See also Thieler, Fried and Rathjens (2016).

Author(s)

Anita M. Thieler with contributions of Uwe Ligges

References

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>

Timmer, J. and König, M. (1995) On Generating Power Law Noise. Astronomy and Astrophysics, 300, 707-710

See Also

Applied in tsgen by TK95_uneq.

Examples

set.seed(31)
# Generate power law noise with exponent alpha=1.5:
y <- TK95(N=2000, alpha=1.5)
tt <- seq(along=y)

# Show time series:
plot(tt,y, type="l", main="Power Law Noise", xlab="t", ylab="y")

# Plot Fourier periodogram with log-axes:
temp <- spectrum(y, plot=FALSE)
plot(log(temp$freq), log(temp$spec), main="log-log-Fourier periodogram",
    xlab="log(frequency)", ylab="log(periodogram)")

# A line with slope -alpha for comparison
abline(a=8, b=-1.5, col="red")
text(-2, 12, expression(alpha==1.5), col="red")

Power law noise generator for unequally sampled observation times

Description

Generates power law noise using TK95 according to Timmer and König (1995), with modifications proposed in Uttley, McHardy and Papadakis (2002) for given irregular observation times.

Usage

TK95_uneq(tt, alpha = 1.5)

Arguments

Value

numeric vector: Noise values related to the observation times.

Note

This function is applied in Thieler et al. (2013). See also Thieler, Fried and Rathjens (2016).

Author(s)

Anita M. Thieler

References

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>

Timmer, J. and König, M. (1995) On Generating Power Law Noise. Astronomy and Astrophysics, 300, 707-710

Uttley, P., McHardy, I. M. and Papadakis, I. E. (2002) Measuring the Broad-Band Power Spectra of Active Galactic Nuclei with RXTE. Monthly Notices of the Royal Astronomical Society, 332 (1), 231-250

See Also

Applies TK95, applied in tsgen.

Examples

# Compare with example in TK95 to see that the power law is much more clear in
# equally sampled data!
set.seed(31)
# Generate power law noise with exponent alpha=1.5:
tt <- sampler(ttype="unif", ps=1, ncycles=2000, npoints=2000)
y <- TK95_uneq(tt, alpha=1.5)

# Show time series:
plot(tt,y, type="l", main="Irregular Power Law Noise", xlab="t", ylab="y")

# Plot Lomb-Scargle periodogram with log-axes:
temp <- RobPer(cbind(tt,y,1), weighting=FALSE, model="sine", regression="L2",
    periods=2000/seq(2, 1000, 2))
plot(log(seq(2, 1000, 2)/2000), log(temp), main="log-log-Fourier periodogram",
    xlab="log(frequency)", ylab="log(periodogram)")
title(main= "Power Law not so obvious", cex.main=0.8, line=0.5)

# A line with slope -alpha for comparison
abline(a=-10, b=-1.5, col="red")
text(-5, -1.5, expression(alpha==1.5), col="red")

Designmatrix generator

Description

This function is used to create the designmatrices needed in RobPer to fit periodic functions. See RobPer or Thieler, Fried and Rathjens (2016) for Details.

Usage

Xgen(tt, n, s, pp, design, steps = 10)

Arguments

Value

numeric matrix: Designmatrix.

Note

A former version of this function is used in Thieler et al. (2013).

Author(s)

Anita M. Thieler and Jonathan Rathjens

References

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>

See Also

Applied in RobPer, see FastTau for an example.

Robust fit of a Beta distribution using CvM distance minimization

Description

Robustly fits a Beta distribution to data using Cramér-von-Mises (CvM) distance minimization.

Usage

betaCvMfit(data, CvM = TRUE, rob = TRUE)

Arguments

Details

betaCvMfit fits a Beta distribution to data by minimizing the Cramér-von-Mises distance. Moment estimates of the parameters of the Beta distribution, clipped to positive values, are used as starting values for the optimization process. They are calculated using

\hat a=-\frac{\bar x \cdot (-\bar x + \bar x^2 + \hat s^2)}{\hat s^2},

\hat b= \frac{\hat a - \hat a \bar x}{\bar x}.

These clipped moment estimates can be returned instead of CvM-fitted parameters setting CvM = FALSE.

The Cramér-von-Mises distance is defined as (see Clarke, McKinnon and Riley 2012)

\frac 1n \sum_{i=1}^n \left(F(u_{(i)}) - \frac{i-0.5}{n}\right)^2+ \frac{1}{12n^2},

where u_{(1)},\ldots,u_{(n)} is the ordered sample and F the distribution function of Beta(a,b).

Value

numeric vector: Estimates for the Parameters a,b of a Beta(a,b) distribution with mean a/(a+b).

Note

Adapted from R-Code from Brenton R. Clarke to fit a Gamma distribution (see Clarke, McKinnon and Riley 2012) using Cramér-von-Mises distance minimization. Used in Thieler et al. (2013). See also Thieler, Fried and Rathjens (2016).

Author(s)

Anita M. Thieler, with contributions from Brenton R. Clarke.

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

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>

See Also

See RobPer-package for an example applying betaCvMfit to detect valid periods in a periodogram.

Examples

# data:
set.seed(12)
PP <- c(rbeta(45, shape1=4, shape2=15), runif(5, min=0.8, max=1))
hist(PP, freq=FALSE, breaks=30, ylim=c(0,7), xlab="Periodogram bar")

# true parameters:
myf.true <- function(x) dbeta(x, shape1=4, shape2=15)
curve(myf.true, add=TRUE, lwd=2)

# 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")

# 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")

# CvM distance minimization
par.CvM <- betaCvMfit(PP, rob=TRUE, CvM=TRUE)
myf.CvM <- function(x) dbeta(x, shape1=par.CvM[1], shape2=par.CvM[2])
curve(myf.CvM, add=TRUE, lwd=2, col="green")

# Searching for outliers...
abline(v=qbeta((0.95)^(1/50), shape1=par.CvM[1], shape2=par.CvM[2]), col="green")

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

Disturbing light curve data

Description

Disturbes a light curve replacing measurement accuracies by outliers and/or observed values by atypical values. See RobPer-package for more information about light curves.

Usage

disturber(tt, y, s, ps, s.outlier.fraction = 0, interval)

Arguments

Details

This function disturbes the light curve (t_i,y_i,s_i)_{i=1,\ldots,n} given. It randomly chooses a proportion of s.outlier.fraction measurement accuracies s_i and replaces them by 0.5\min(s_1,\ldots,s_n). In case of interval=TRUE a time interval [t_{start},t_{start}+3p_s] within the intervall [t_1,t_n] is randomly chosen and all observed values belonging to this time interval are replaced by a peak function:

y_i^{changed} = 6 \ \tilde y_{0.9}\ \frac{d_{\mathcal N(t_{start}+1.5p_s, p_s^2)}(t_i) }{ d_{\mathcal N(0,p_s^2)}(0)} \quad \forall \ i \ : \ t_i\in[t_{start}, t_{start}+3p_s],

where d_{\mathcal N(a,b^2)}(x) denotes the density of a normal distribution with mean a and variance b^2 at x.

In case of s.outlier.fraction=0 and interval=FALSE, y and s are returned unchanged.

Value

Note

A former version of this function is used in Thieler et al. (2013). See also Thieler, Fried and Rathjens (2016).

Author(s)

Anita M. Thieler

References

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>

See Also

Applied in tsgen (see there for example).

Noise and measurement accuracy generator for light curves

Description

Generates measurement accuracies, a white noise component depending on them and a second (possibly power law, i.e. red) noise component which does not depend on the measurement accuracies. For more details see tsgen or Thieler, Fried and Rathjens (2016). See RobPer-package for more information about light curves.

Usage

lc_noise(tt, sig, SNR, redpart, alpha = 1.5)

Arguments

Value

Note

A former version of this function is used in Thieler et al. (2013).

Author(s)

Anita M. Thieler and Jonathan Rathjens

References

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>

See Also

Applied in tsgen (see there for an example), applies TK95_uneq.

Generator for irregularly sampled observation times

Description

Generates irregularly sampled observation times with a periodic sampling pattern

Usage

sampler(ttype, npoints, ncycles, ps = 1)

Arguments

Details

sampler generates observation times t_1,\ldots,t_n with a periodic sampling of period p_s. Four distributions are possible: In case of ttype="equi", the t_i are equidistantly sampled with t_i=i\frac{p_sn_s}{n}. For ttype="unif", the observation times are independently drawn form a uniform distribution on [0,n_sp_s]. Both these sampling schemes are aperiodic, the sampling period p_s only influences the length t_n-t_1 of the series of observation times.

For ttype="sine" and ttype="trian", observation cycles z^\star_i are drawn from a uniform distribution on \{1,\ldots,n_s\} and observation phases \varphi^\star_i are drawn from a density

d_{sine}(x)= \sin(2\pi x)+1

(for ttype="sine") or

d_{trian}(x)= 3x, \quad 0\leq x\leq\frac{2}{3},

d_{trian}(x)= 6-6x, \quad \frac{2}{3}<x\leq 1

(for ttype="trian"). The unsorted observation times t^\star_i are then generated using

t^\star_i= \varphi^\star_i+(z^\star_i-1)p_s.

Separately sampling observation cycle and phase was proposed by Hall and Yin (2003). For more details see Thieler, Fried and Rathjens (2016) or Thieler et al. (2013).

Value

numeric vector: Ordered observation times.

Note

To sample from d_{sine}, the function BBsolve, package BB, is used.

A former version of this function is used in Thieler et al. (2013).

Author(s)

Anita M. Thieler and Jonathan Rathjens

References

Hall, P. and Yin, J. (2003): Nonparametric Methods for Deconvolving Multiperiodic Functions. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 65 (4), 869-886

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>

See Also

Applied in tsgen (see there for an example).

Generator for periodic signal in a light curve

Description

Calculates periodically varying values for given observation times.

Usage

signalgen(tt, ytype, pf = 1)

Arguments

Details

The values y_{f;1},\ldots,y_{f;n} with fluctuation period p_f and related to observation times t_1,\ldots,t_n are generated using

y_{f;i}=f\left(\frac{t_i}{p_f}\right), i=1,\ldots,n.

Depending on ytype (see above), f is defined as:

f_{const}(t) = 0,

f_{sine}(t)= \sin\left(\frac{2\pi t}{p_f}\right),

f_{trian}(t)= 3\varphi_{1}(t), \quad 0\leq \varphi_{1}(t)\leq\frac{2}{3},

f_{trian}(t)= 6-6\varphi_{1}(t),\quad \frac{2}{3}<\varphi_{1}(t)\leq 1,

f_{peak}(t)= 9\exp\left(-3p_f^2\left(\varphi_{1}(t)-\frac 23\right)^2\right),\quad 0\leq \varphi_{1}(t)\leq\frac{2}{3},

f_{peak}(t)= 9\exp\left(-12p_f^2\left(\varphi_{1}(t)-\frac 23\right)^2\right),\quad \frac{2}{3}<\varphi_{1}(t)\leq 1,

with \varphi_1(t) = t mod1 = (t-\lfloor t/p_f \rfloor p_f)/p_f = (t%%1)/pf. f_{const} means that there is no (periodic) fluctuation, f_{sine} defines a sine function, f_{trian} defines a triangular shaped periodic function and f_{peak} a periodically repeating peak.

Value

numeric vector: Values y_{f;1},\ldots,y_{f;n}.

Note

This function is used in Thieler et al. (2013). See also Thieler, Fried and Rathjens (2016).

Author(s)

Anita M. Thieler and Jonathan Rathjens

References

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>

See Also

Applied in tsgen (see there for an example).

Data: Light curve from GROJ0422+32

Description

Light curve for gamma ray emission of the source GROJ0422+32.

Usage

star_groj0422.32

Format

A matrix of three columns, with a time series of length 729 appropriate to RobPer.

Details

Data obtained by the BATSE Earth Occultation Monitoring project of the NSSTC, available from https://gammaray.nsstc.nasa.gov/batse/occultation/. The experiments are described in Harmon et al. (2002) and Harmon et al. (2004).

Note

See Vignette Section 5.2 for example.

Source

Data kindly provided by the National Aeronautics and Space Administration (NASA); National Space, Science, and Technology Center (NSSTC); Gamma-Ray Astrophysics Team (see Details).

References

Harmon, B. A., Fishman, G. J., Wilson, C. A., Paciesas, W. S., Zhang, S. N., Finger, M. H., Koshut, T. M., McCollough, M. L., Robinson, C. R. and Rubin, B. C. (2002): The Burst and Transient Source Experiment Earth Occultation Technique. The Astrophysical Journal Supplement Series, 138 (1), 149-183

Harmon, B. A., Wilson, C. A., Fishman, G. J., Connaughton, V., Henze, W., Paciesas, W. S., Finger, M. H., McCollough, M. L., Sahi, M., Peterson, B., Shrader, C. R., Grindlay, J. E. and Barret, D. (2004): The Burst and Transient Source Experiment (BATSE) Earth Occultation Catalog of Low-Energy Gamma-Ray Sources. The Astrophysical Journal Supplement Series, 154 (2), 585-622

Artificial light curve generator

Description

This function generates light curves (special time series) with unequally sampled observation times, different periodicities both in sampling and observed values, with white and power law (red) noise in the observed values and possibly disturbed observations. See RobPer-package for more information about light curves and also Thieler, Fried and Rathjens (2016) for more details in general.

Usage

tsgen(ttype, ytype, pf, redpart, s.outlier.fraction = 0, interval, npoints,
 ncycles, ps, SNR, alpha = 1.5)

Arguments

Details

tsgen generates an artificial light curve consisting of observation times t_1,\ldots,t_n, observation values y_1,\ldots,y_n and measurement accuracies s_1,\ldots,s_n. It calls several subfunctions (see there for details):

sampler is used to sample observation times t_1,\ldots,t_n in the interval [0,n_s*p_s] with a possibly periodic sampling of period p_s.

signalgen generates periodically varying values y_{f;1},\ldots,y_{f;n} at time points t_1,\ldots,t_n with fluctuation period p_f.

lc_noise samples measurement accuracies s_1,\ldots,s_n from a Gamma(3,10)-distribution and a white noise component y_{w;1},\ldots,y_{w;n} with from \mathcal N(0,s_i^2) distributions. A second noise component y_{r;1},\ldots,y_{r;n} does not depend on the s_i. It is generated as red noise, i.e. following a power law with power law index \alpha. For white noise choose \alpha=0, for flicker noise (pink noise) \alpha=1, for brown noise \alpha=2. The power law noise is generated using TK95_uneq and TK95. The noise components are scaled so that the variance of the y_{r;i} has approximately the proportion redpart in the overall noise variance and that SNR is the ratio var(y_f)/var(y_w+y_r). The observed values are set to y_i= y_{f;i}+y_{w;i}+y_{r;i} \forall i.

disturber disturbes the light curve replacing measurement accuracies s_i by outliers (if s.outlier.fraction>0) and observed values y_i by atypical values (if interval=TRUE). In case of s.outlier.fraction=0 and interval=FALSE, the function returns all values unchanged.

Value

Note

Note that the white noise components' variances are exactly s_i^2, so the s_i are no estimates, but true values. In this sense, the measurement accuracies of a generated light curve are more informative than for real light curves, where the measurement accuracies are estimates, see Thieler et al. (2013), where also a former version of this function is applied.

To lower the informativity of the measurement accuracies, set redpart to a strictly positive value, possibly with alpha=0 if no other noise components than white ones are required.

Author(s)

Anita M. Thieler and Jonathan Rathjens

References

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>

See Also

Applies sampler, signalgen, lc_noise, disturber, TK95, TK95_uneq.

Examples

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

# Or do it step by step:
# First sampling observation times:
set.seed(22)
tt <- sampler(ttype="sine", npoints=200, ncycles=100, ps=5)

# show obviously irregular observation times as grey vertical bars on a red time line:
plot(tt, tt*0, type="n", axes=FALSE, xlab="Observation Times", ylab="")
    abline(v=tt, col="grey")
axis(1, pos=0, col="red", col.axis="red")

# Sampling period is 5, look at observation times modulo 10:
hist(tt%%5, xlab="Observation time modulo 5", 
    main="Sine Distribution for the phase (tt modulo 5)", freq=FALSE)
dsin <- function(tt) 0.2*(sin(2*pi*tt/5)+1)
curve(dsin, add=TRUE)

# Then generate periodic fluctuation
yf <- signalgen(tt, ytype="peak", pf=7)

plot(tt, yf, xlab="Observation Times", ylab="Periodic Fluctuation")
plot(tt%%7, yf, main="Phase Diagram (time modulo 7)", 
    xlab="Observation time modulo 7",  ylab="Periodic Fluctuation")

# Add noise and scale signal to the right SNR
temp <- lc_noise(tt,sig=yf, SNR=3, redpart=0.1, alpha=0)
y <- temp$y
s <- temp$s

# Plotting the light curve (vertical bars show measurement accuracies)
plot(tt, y, pch=16, cex=0.5, xlab="t", ylab="y", main="a Light Curve")
rect(tt, y+s, tt, y-s)

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

# replace measurement accuracies by tiny outliers or include a peak
temp <- disturber(tt,y,s,ps=5, s.outlier.fraction=0, interval=FALSE)

# Phase diagram (observation times modulo 10)
plot(tt%%7, temp$y, pch=16, cex=0.5, xlab="t", ylab="y", 
    main="Phase Diagram of a Light Curve")
rect(tt%%7, temp$y+temp$s, tt%%7, temp$y-temp$s)

# The result is the same:
all(cbind(tt,temp$y,temp$s)==lightcurve)