In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import seaborn; seaborn.set()
In [2]:
from multiband_LS.lomb_scargle import LombScargle
In [3]:
from astroML.datasets import fetch_LINEAR_sample
data = fetch_LINEAR_sample()
lcid = 10022663
t, y, dy = data[lcid].T
In [4]:
periods = np.linspace(0.6, 0.65, 10000)
omegas = 2 * np.pi / periods
In [5]:
Nterms = 4
model = LombScargle(Nterms=Nterms).fit(t, y, dy)
P = model.power(omegas)
In [6]:
reg = 0.1 * np.ones(2 * Nterms + 1)
reg[:5] = 0 # no regularization on low-order terms
model_reg = LombScargle(Nterms=Nterms, regularization=reg).fit(t, y, dy)
P_reg = model_reg.power(omegas)
In [7]:
plt.plot(periods, P, label='LS')
plt.plot(periods, P_reg, label='regularized LS');
In [9]:
P_best = periods[np.argmax(P)]
omega_best = omegas[np.argmax(P)]
t_fit = np.linspace(0, P_best, 300)
plt.errorbar(t % P_best, y, dy, fmt='o')
plt.plot(t_fit, model.predict(t_fit, omega_best), label='LS')
plt.plot(t_fit, model_reg.predict(t_fit, omega_best), label='regularized LS')
plt.gca().invert_yaxis()
plt.legend();