Example of a Multi-band Lomb-Scargle Fit

Here's a quick example of doing a multi-band lomb-scargle fit to one of the stripe 82 RR Lyrae stars from Sesar 2010.

First some standard Python imports and settings:



In [1]:

    
from __future__ import print_function, division

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

# Use seaborn plot style defaults
import seaborn; seaborn.set()

Fetching the data

We'll start by fetching observed light curves from the Sesar 2010 Paper:



In [2]:

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

from multiband_LS.data import fetch_rrlyrae
rrlyrae = fetch_rrlyrae()
lcids = np.fromiter(rrlyrae.ids, dtype=int)

Let's plot the first set of light curves



In [3]:

    
t, y, dy, filts = rrlyrae.get_lightcurve(lcids[0], return_1d=True)

for i, filt in enumerate('ugriz'):
    mask = (filts == filt)
    plt.errorbar(t[mask], y[mask], dy[mask], fmt='o', label=filt)

plt.gca().invert_yaxis()
plt.legend(ncol=3, loc='upper left', fontsize=12)
plt.xlabel('observation date'); plt.ylabel('magnitude');

Doing a multi-band fit

The LombScargleMultiband class has the tools needed to compute multiband lomb-scargle powers, and to fit the best model based on these.

Note that we're doing a four-term truncated Fourier fit:



In [4]:

    
from multiband_LS import LombScargleMultibandFast
model = LombScargleMultibandFast()
model.optimizer.period_range = (0.4, 0.8)
model.fit(t, y, dy, filts)









    Out[4]:





<multiband_LS.lomb_scargle_multiband.LombScargleMultibandFast at 0x103790f10>

Let's find the best period:



In [5]:

    
print("best period: {0:.4f}".format(model.best_period))









    



Finding optimal frequency:
 - Using omega_step = 0.00038
 - Computing periods at 20758 steps from 0.40 to 0.80
Zooming-in on 5 candidate peaks:
 - Computing periods at 1005 steps
best period: 0.6143

Let's zoom-in around this best value:



In [6]:

    
periods = np.linspace(model.best_period - 0.01,
                      model.best_period + 0.01, 1000)

P = model.periodogram(periods)
plt.plot(periods, P);

Now given this optimal period, we can plot the folded light curves along with the model fit:



In [7]:

    
def plot_folded(model):
    tfit = np.linspace(0, model.best_period, 100)[:, None]
    yfits = model.predict(tfit, filts=['u', 'g', 'r', 'i', 'z'])

    colors = plt.rcParams['axes.color_cycle']

    for i, (color, filt) in enumerate(zip(colors[:5], 'ugriz')):
        mask = (filts == filt)
        plt.errorbar(t[mask] % model.best_period, y[mask], dy[mask], fmt='o', c=color)
        plt.plot(tfit, yfits[:, i], '-', c=color, alpha=0.6, label=filts[i])

    plt.gca().invert_yaxis()
    plt.legend()
    
plot_folded(model)

From this model, we can plot the trend in color-color space as a function of phase:



In [8]:

    
def plot_color_mag(model):
    tfit = np.linspace(0, model.best_period, 100)[:, None]
    yfits = model.predict(tfit, filts=['u', 'g', 'r', 'i', 'z'])
    
    g, dg = y[filts == 'g'], dy[filts == 'g']
    r, dr = y[filts == 'r'], dy[filts == 'r']
    dgr = np.sqrt(dg ** 2 + dr ** 2)

    gfit = yfits[:, 1]
    rfit = yfits[:, 2]

    plt.errorbar(g - r, r, yerr=dg, xerr=dgr, fmt='o')
    plt.plot(gfit - rfit, rfit, '-k', alpha=0.7)
    plt.xlabel('g - r'); plt.ylabel('r');
    
plot_color_mag(model)

This can also be repeated with a multi-term truncated Fourier series. Let's do four terms:



In [9]:

    
model = LombScargleMultibandFast(Nterms=4)
model.optimizer.period_range = (0.6, 0.65)
model.fit(t, y, dy, filts)









    Out[9]:





<multiband_LS.lomb_scargle_multiband.LombScargleMultibandFast at 0x10adbb690>



In [10]:

    
print("best period: {0:.4f} days".format(model.best_period))









    



Finding optimal frequency:
 - Using omega_step = 0.00038
 - Computing periods at 2130 steps from 0.60 to 0.65
Zooming-in on 5 candidate peaks:
 - Computing periods at 1005 steps
best period: 0.6143 days



In [11]:

    
plot_folded(model)



In [12]:

    
plot_color_mag(model)

Other Multiband Method...

Here's the more flexible multiband method. Here we can fit a certain number of terms for the base model, and then a certain number of terms for the color offsets. A very light regularization is used to make sure that the models are well-behaved, and to drive the bulk of the variation into the base model.



In [13]:

    
from multiband_LS import LombScargleMultiband
model = LombScargleMultiband(Nterms_base=4, Nterms_band=2)
model.optimizer.period_range = (0.6, 0.65)
model.fit(t, y, dy, filts)









    Out[13]:





<multiband_LS.lomb_scargle_multiband.LombScargleMultiband at 0x10b09b390>



In [14]:

    
print("best period: {0:.4f}".format(model.best_period))









    



Finding optimal frequency:
 - Using omega_step = 0.00038
 - Computing periods at 2130 steps from 0.60 to 0.65
Zooming-in on 5 candidate peaks:
 - Computing periods at 1005 steps
best period: 0.6143



In [15]:

    
plot_folded(model)



In [16]:

    
plot_color_mag(model)