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