Lomb-Scargle Example

Here is an example of performing a lomb-scargle periodogram on a LINEAR light curve.



In [1]:

    
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

# use seaborn's default plotting styles for matplotlib
import seaborn; seaborn.set()

The Data

The astroML project includes a data fetcher for light curves from the LINEAR survey. We'll use one of the light curves as an example.



In [2]:

    
from astroML.datasets import fetch_LINEAR_sample
data = fetch_LINEAR_sample()









    



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



In [3]:

    
print(data.ids)









    



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



In [4]:

    
lcid = 10022663
t, y, dy = data[lcid].T



In [5]:

    
fig, ax = plt.subplots()
ax.errorbar(t, y, dy, fmt='o');
ax.invert_yaxis()
ax.set_xlabel('MJD')
ax.set_ylabel('magnitude')
ax.set_title('LINEAR object {0}'.format(lcid));

Computing the Periodogram

Now we import the functions to compute the periodogram:



In [6]:

    
# imports are from up one directory
import os, sys; sys.path.append(os.path.abspath('..'))

from multiband_LS.lomb_scargle import LombScargle, LombScargleAstroML

First we'll do the periodogram with a single-term model:



In [7]:

    
model1 = LombScargle(fit_offset=True, Nterms=1).fit(t, y, dy)
print("best period: {0:.4f} days".format(model1.best_period))









    



Finding optimal frequency:
 - Using omega_step = 0.00064
 - Computing periods at 40956 steps from 0.20 to 1.20
Zooming-in on 5 candidate peaks:
 - Computing periods at 1010 steps
best period: 0.6160 days



In [8]:

    
periods = np.linspace(model1.best_period - 0.005,
                      model1.best_period + 0.005, 5000)

P = model1.periodogram(periods)
plt.plot(periods, P)
plt.xlabel('period (days)'); plt.ylabel('Lomb-Scargle Power');

Now with a 6-term truncated Fourier model:



In [9]:

    
model6 = LombScargle(fit_offset=True, Nterms=6).fit(t, y, dy)
print("best period: {0:.4f} days".format(model6.best_period))









    



Finding optimal frequency:
 - Using omega_step = 0.00064
 - Computing periods at 40956 steps from 0.20 to 1.20
Zooming-in on 5 candidate peaks:
 - Computing periods at 1010 steps
best period: 0.6160 days



In [10]:

    
P = model6.periodogram(periods)

plt.plot(periods, P)
plt.xlabel('period (days)'); plt.ylabel('Lomb-Scargle Power');

Let's visualize the best-fit models beside each other:



In [11]:

    
t_fit = np.linspace(0, model6.best_period, 1000)

plt.errorbar(t % model6.best_period, y, dy, fmt='o')
for model in [model1, model6]:
    y_fit = model.predict(t_fit)
    plt.plot(t_fit, y_fit, label="{0} terms".format(model.Nterms))
plt.gca().invert_yaxis()
plt.legend(loc='best');

AstroML Version

For convenience, we've wrapped the astroML version of the lomb-scargle algorithm. Note that this uses a faster algorithm, but is not as general (i.e. it cannot fit multiple terms). It also does not have a predict() functionality



In [12]:

    
model_LS = LombScargleAstroML().fit(t, y, dy)
print("best period: {0:.4f} days".format(model_LS.best_period))









    



Finding optimal frequency:
 - Using omega_step = 0.00064
 - Computing periods at 40956 steps from 0.20 to 1.20
Zooming-in on 5 candidate peaks:
 - Computing periods at 1010 steps
best period: 0.6160 days



In [13]:

    
P = model_LS.periodogram(periods)

plt.plot(periods, P)
plt.xlabel('period (days)'); plt.ylabel('Lomb-Scargle Power');

The results of this will be identical (within machine precision) to the 1-term model above.