Using Supersmoother for Periodic Fits

This is an example of using the supersmoother algorithm of Riemann 1995 to find periodic signals in time series data.



In [1]:

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

# use seaborn for matplotlib settings
import seaborn; seaborn.set()

The Data

We'll use astroML to fetch some data from the LINEAR survey. AstroML can be installed using pip install astroML



In [2]:

    
from astroML.datasets import fetch_LINEAR_sample
data = fetch_LINEAR_sample()









    



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



In [3]:

    
lcid = 10022663
t, y, dy = data[lcid].T
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));

Fitting the Period

Now we use the supersmoother period fitter as follows:



In [4]:

    
# add up-directory to Python path
import sys, os; sys.path.append(os.path.abspath('..'))
from periodogram.supersmoother import SuperSmoother



In [5]:

    
model = SuperSmoother().fit(t, y, dy)
periods = np.linspace(0.2, 0.8, 10000)
scores = model.score(periods)

plt.plot(periods, scores);

There is a clear peak near 0.61! Let's zoom-in and take a closer look:



In [6]:

    
periods = np.linspace(0.61, 0.62, 1000)
scores = model.score(periods)

plt.plot(periods, scores);

For good-measure, let's do another zoom and make sure we have a really precise maximum:



In [7]:

    
periods = np.linspace(0.615, 0.617, 1000)
scores = model.score(periods)
best_period = periods[np.nanargmax(scores)]
print(best_period)









    



0.615958958959

Now we can plot the folded light-curve and the best-fit supersmoother model:



In [8]:

    
tfit = np.linspace(0, best_period, 1000)
plt.errorbar(t % best_period, y, dy, fmt='o')
plt.plot(tfit, model.predict(tfit, period=best_period), '-k');

Ta daa!! Now it would be useful to compare this to the Lomb Scargle and see what the results are.