SysRem

SysRem is an algorithm to remove systematic effects, defined as those which occur similarly (linearly) across a set of observations such as sets of light curves or spectra. The algorithm was described by Tamuz et al. 2005 (MNRAS 356, 1466) in the context of correcting systematic effects in samples of light curves. Repeated application allows to remove more complex effects but the danger of removing actual signal tends to increase with the number of iteration.

Note

The data (or residuals) are stored such that X[::,0] gives the first observation, i.e., the observation is stored in the first column. The data matrix is, therefore, stored here as the transpose of that adopted for the PCA etc..

Note

Compared to the presentation by Tamuz et al., the roles of a and c are exchanged.

Example: Application to mock data

from PyAstronomy import pyasl
import numpy as np
import matplotlib.pylab as plt

# No. of data sets (e.g., light curves)
nds = 50
# No. of data points per data set
nobs = 200

obs, sigs = [], []

# Generate mock data
x = np.linspace(0,1,nobs)-0.5
for i in range(nds):
    # Add noise
    y = np.random.normal(0,0.01+0.0001*i,nobs)
    # Add 3rd degree polynomial
    p = (i, 0.2*i, 0.3*i)
    y += np.polyval(p,x)
    # Add moving Gaussian signal
    y -= 0.02 * np.exp(-(x-(i/nds-0.5))**2/(2*0.02**2))
    # Add growing but stationary Gaussian signal
    y -= (i**2*0.05) * np.exp(-x**2/(2*0.06**2))
    obs.append(y)
    sigs.append(np.ones_like(y)*0.01+0.0001*i)

    plt.plot(x, y, '.-')
plt.title("Mock data")
plt.show()

sr = pyasl.SysRem(obs, sigs)
# First iteration
r, a, c = sr.iterate()
plt.subplot(2,1,1)
plt.title("Residuals after first (top) and second iteration (bottom)")
plt.imshow(r, origin='lower', aspect="auto")
# Second iteration
r, a, c = sr.iterate()
plt.subplot(2,1,2)
plt.imshow(r, origin='lower', aspect="auto")
plt.show()

Example: SYSREM, PCA, and eigenvectors

API

Interface class

class PyAstronomy.pyasl.SysRem(obs, sigs, ms_obs=True, ms_feat=False, a0=None, ms_warn=True)

Implementation of the SysRem algorithm.

SysRem was described by Tamuz et al. 2005 (MNRAS 356, 1466) in the context of correcting systematic effects in samples of light curves, but has been applied in other areas such as planetary atmospheric spectroscopy.

Note

The data (or residuals) are stored such that X[::,0] gives the first observation, i.e., the observation is stored in the first column. The data matrix is, therefore, stored here as the transpose of that adopted for the PCA etc..

Note

Compared to the presentation by Tamuz et al., the roles of a and c are exchanged.

Parameters
obslist of 1d arrays or 2d array

List of observations (e.g., light curves) or 2d array such that obs[::,0] is the first observation.

sigslist of 1d arrays or 2d array

Uncertainties (same format as obs)

a01d array, optional

Starting values for ‘a’ parameter (same length as the observations). If not given, values linearly increasing from 0 to 1 are assumed.

ms_obsboolean, optional

Subtract mean from observations in initial data matrix (columns of data matrix). Default is True.

ms_featboolean, optional

Subtract mean from features in initial data matrix (e.g., spectral bins). Default is False.

ms_warnboolean, optional

If True (default), a warning will be printed if either ms_obs or ms_feat is True so that the data will be manipulated by subtracting mean(s). Set False to suppress warning.

Attributes
rijslist of 2d arrays

List of residual arrays. Each call of the iterate method adds an updated residual array to the list.

aclist of tuples of arrays

The a and c parameters used to obtain the updated residual array from the previous one.

Methods

get_cumulative_model()

Get the cumulative model (sum of all available models)

get_model([last, a, c])

Get model corresponding to a,c vector

iterate([atol, rtol, imax, a0])

A single SysRem iteration: Remove linear systematic effect
get_cumulative_model()

Get the cumulative model (sum of all available models)

get_model(last=True, a=None, c=None)

Get model corresponding to a,c vector

Parameters
lastboolean, optional

If True (default), the model for the last iteration will be returned (a,c must not be given then).

a, carrays, optional

If given, these will be used to calculate the model (last must be set False)

iterate(atol=0.001, rtol=0, imax=1001, a0=None)

A single SysRem iteration: Remove linear systematic effect

Parameters
atol, rtolfloat, optional

Absolute and relative tolerance for a-c iteration required to stop. The pertinent quantity tested is the difference between consecutive models divided by the uncertainties. By default, relative tolerance is ignored.

imaxint

Maximum of iterations. Throws exception if convergence is not reached earlier.

a0array, optional

Starting value for SYSREM model PCA-like component. If not specified, the value defined during instance initialization will be used.

Returns
rij2d array

Updated residual matrix with best-fit model subtracted.

aarray

Best-fit ‘a’ values

carray

Best-fit ‘c’ values

Helper functions

PyAstronomy.pyasl.sysrem_iter_c(rij, sigij, a, sigij2=None)

Estimate ‘c’

Parameters
rij2d array

Residuals (rij[::,i] = ith of nds observation with nobs data points)

sigij2d array

Uncertainties

a1d array

Current estimate of ‘a’ (length is nobs)

sigij22d array, optional

Square of uncertainties. If given, sigij will be ignores.

Returns
c1d array

Estimate of ‘c’

PyAstronomy.pyasl.sysrem_iter_a(rij, sigij, c, sigij2=None)

Estimate ‘a’

Parameters
rij2d array

Residuals (rij[::,i] = ith of nds observation with nobs data points)

sigij2d array

Uncertainties

c1d array

Current estimate of ‘a’ (length is nds)

sigij22d array, optional

Square of uncertainties. If given, sigij will be ignores.

Returns
a1d array

Estimate of ‘a’

PyAstronomy.pyasl.sysrem_data_prepare(obs, sigs, ms_obs=True, ms_feat=False)

Construct data sets for SysRem

Parameters
obslist of arrays

Observations

sigslist of arrays, optional

Uncertainties of the observations

Returns
rij2d array

The residuals (mean-subtracted observations) so that rij[::,i] gives the i-th observation

sm2d array

The uncertainties in the same format as rij

a01d array

A dummy guess for the ‘a’ vector (linear between 0 and 1)

ms_obsboolean, optional

Subtract mean from observations (columns of data matrix). Default is True.

ms_featboolean, optional

Subtract mean from features (e.g., spectral bins). Default is False.

PyAstronomy.pyasl.sysrem_update_rij(rij, a, c)

Subtract current model from residuals

In matrix notation Tamuz et al. write (Eq. 6) R’ = R - c a^T. The operation here implemented corresponds to (R’ = R - c^T)^T -> R’^T = R^T - a c^T (where ^T is the transpose, R is a matrix and a and c are vectors).

Returns
New rij2d array

Updated residuals