Robust Lomb-Scargle


In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

import seaborn as sns
sns.set()

In [2]:
# Create some data
np.random.seed(42)

N = 100

t = 10 * np.random.rand(N)
dy = 0.1 + 0.1 * np.random.rand(N)
y = np.random.normal(np.sin(t) + 0.5 * np.sin(2 * t), dy)

In [3]:
from bombscargle.standard import lomb_scargle, lomb_scargle_huber
from astroML.time_series import lomb_scargle as lomb_scargle_astroml

In [4]:
# Plot the data
plt.errorbar(t, y, dy, fmt='.k', ecolor='gray');



In [5]:
period = np.linspace(0.5, 20, 1000)
omega = 2 * np.pi / period

# Standard Lomb-Scargle result
PLS = lomb_scargle(t, y, dy, omega, Nterms=2,
                   compute_offset=False)

In [6]:
plt.plot(period, PLS)
plt.ylim(0, 1);



In [7]:
# Add some outliers
i = np.random.random(len(y))
mask = (i < 0.2)
y[mask] += 2 * np.random.randn(mask.sum())

In [8]:
# Plot the data
plt.errorbar(t, y, dy, fmt='.k', ecolor='gray');



In [9]:
period = np.linspace(0.5, 20, 1000)
omega = 2 * np.pi / period

# Standard Lomb-Scargle result
PLS = lomb_scargle(t, y, dy, omega, Nterms=1,
                   compute_offset=False)

# Huber Lomb-Scargle result
PLS_huber = lomb_scargle_huber(t, y, dy, omega, Nterms=1, c=5,
                               compute_offset=False)

PLS_astroML = lomb_scargle_astroml(t, y, dy, omega, generalized=False)

In [10]:
plt.plot(period, PLS, label='standard')
plt.plot(period, PLS_huber, label='huber')
#plt.plot(period, PLS_astroML, label='astroML')
plt.legend(frameon=False);


Real Data

Let's use some LINEAR data from astroML:


In [11]:
from astroML.datasets import fetch_LINEAR_sample
data = fetch_LINEAR_sample()


[10003298 10004892 10013411 ...,  9984569  9987252   999528]

In [ ]:
t, y, dy = data.get_light_curve(10022663).T
plt.errorbar(t, y, dy, fmt='o');



In [ ]:
period = np.linspace(0.615, 0.618, 1000)
omega = 2 * np.pi / period

P1 = lomb_scargle(t, y, dy, omega, compute_offset=False)
P2 = lomb_scargle_huber(t, y, dy, omega, compute_offset=False)
P3 = lomb_scargle_astroml(t, y, dy, omega, generalized=False)

In [ ]:
plt.plot(period * 24, P1, label='lomb scargle')
plt.plot(period * 24, P2, label='Huber')
#plt.plot(period * 24, P3, label='astroML version')
plt.legend();

Let's look at a curve we know has outliers...


In [ ]:
t, y, dy = data[1004849].T
plt.errorbar(t, y, dy, fmt='o');

In [ ]:
period = np.linspace(0.2, 0.7, 10000)
omega = 2 * np.pi / period
P = lomb_scargle_astroml(t, y, dy, omega)
plt.plot(period, P);

In [ ]:
period = np.linspace(0.457, 0.46, 1000)
omega = 2 * np.pi / period
P = lomb_scargle_astroml(t, y, dy, omega)
plt.plot(period, P);

In [ ]:
period_best = period[np.argmax(P)]
phase = t % period_best
plt.errorbar(phase, y, dy, fmt='.', ecolor='#AABBFF')
plt.gca().invert_yaxis();

In [ ]:
P1 = lomb_scargle(t, y, dy, omega, compute_offset=False)
P2 = lomb_scargle_huber(t, y, dy, omega, compute_offset=False)

plt.plot(period * 24, P1, label='lomb scargle')
plt.plot(period * 24, P2, label='Huber')
plt.legend();

In [ ]: