Source code for astroML.time_series.ACF

"""
Auto-correlation functions
"""
import numpy as np
from scipy import fftpack
from astropy.timeseries import LombScargle


[docs]def ACF_scargle(t, y, dy, n_omega=2 ** 10, omega_max=100): """Compute the Auto-correlation function via Scargle's method Parameters ---------- t : array_like times of observation. Assumed to be in increasing order. y : array_like values of each observation. Should be same shape as t dy : float or array_like errors in each observation. n_omega : int (optional) number of angular frequencies at which to evaluate the periodogram default is 2^10 omega_max : float (optional) maximum value of omega at which to evaluate the periodogram default is 100 Returns ------- ACF, t : ndarrays The auto-correlation function and associated times """ t = np.asarray(t) y = np.asarray(y) if y.shape != t.shape: raise ValueError("shapes of t and y must match") dy = np.asarray(dy) * np.ones(y.shape) d_omega = omega_max * 1. / (n_omega + 1) omega = d_omega * np.arange(1, n_omega + 1) # recall that P(omega = 0) = (chi^2(0) - chi^2(0)) / chi^2(0) # = 0 # compute P and shifted full-frequency array ls = LombScargle(t, y, dy) frequency = np.asarray(omega) / (2 * np.pi) P = ls.power(frequency) P = np.concatenate([[0], P, P[-2::-1]]) # compute PW, the power of the window function lsw = LombScargle(t, np.ones(len(t)), dy, center_data=False) PW = lsw.power(frequency) PW = np.concatenate([[0], PW, PW[-2::-1]]) # compute the inverse fourier transform of P and PW rho = fftpack.ifft(P).real rhoW = fftpack.ifft(PW).real ACF = fftpack.fftshift(rho / rhoW) / np.sqrt(2) N = len(ACF) dt = 2 * np.pi / N / (omega[1] - omega[0]) t = dt * (np.arange(N) - N // 2) return ACF, t
[docs]def ACF_EK(t, y, dy, bins=20): """Auto-correlation function via the Edelson-Krolik method Parameters ---------- t : array_like times of observation. Assumed to be in increasing order. y : array_like values of each observation. Should be same shape as t dy : float or array_like errors in each observation. bins : int or array_like (optional) if integer, the number of bins to use in the analysis. if array, the (nbins + 1) bin edges. Default is bins=20. Returns ------- ACF : ndarray The auto-correlation function and associated times err : ndarray the error in the ACF bins : ndarray bin edges used in computation """ t = np.asarray(t) y = np.asarray(y) if y.shape != t.shape: raise ValueError("shapes of t and y must match") if t.ndim != 1: raise ValueError("t should be a 1-dimensional array") dy = np.asarray(dy) * np.ones(y.shape) # compute mean and standard deviation of y w = 1. / dy / dy w /= w.sum() mu = np.dot(w, y) sigma = np.std(y, ddof=1) dy2 = dy[:, None] dt = t - t[:, None] UDCF = ((y - mu) * (y - mu)[:, None] / np.sqrt((sigma ** 2 - dy ** 2) * (sigma ** 2 - dy2 ** 2))) # determine binning bins = np.asarray(bins) if bins.size == 1: dt_min = dt.min() dt_max = dt.max() bins = np.linspace(dt_min, dt_max + 1E-10, bins + 1) ACF = np.zeros(len(bins) - 1) M = np.zeros(len(bins) - 1) for i in range(len(bins) - 1): flag = (dt >= bins[i]) & (dt < bins[i + 1]) M[i] = flag.sum() ACF[i] = np.sum(UDCF[flag]) ACF /= M return ACF, np.sqrt(2. / M), bins