Making sense of the lomb-scargle periodogram

Question

I am trying to use the periodogram to tell when a signal is periodic or not by following the tutorial for the astropy Lomb-scargle periodogram here.

http://docs.astropy.org/en/stable/stats/lombscargle.html

I simulated some data, one that is a sinusoid (period = 200) and one which is a skew gaussian (i.e, a single transient event). The hope was that periodogram would pick out the period for the periodic object and and give a period for the transient that would imply a single event occured in the window.

Unfortunately, the results don't make sense at all. I have attached code at the end and the figures generated below. Due to random noise in my simulation, each result is different and I provide two examples of the results below. I use the same method outlined in the link above where we use the function LombScargle.model() on the best fit frequency from frequency[argmax[power]]

The red line is the true function I simulated the data from. The green is the best fit from the periodogram. The right hand plots are the PSD from the periodogram.

Example 1

Here we can see the best fit frequency for the sinusoid (top right plot) picked out is 0.105 (i.e. a period of 9-10 days) which is not near a frequency of 0.05 I'd expect for something with a period of 200 days, yet when I feed this best fit frequency of 0.105 to the lomb-scargle model fitter, a nicely matching periodic curve comes out with a period of 200 days?

This does not make sense.

Example 2

Here I ran the code again and this time the results are switched around? It fit a very large period to the transient so that I can confidently say it is a single transient event yet the sinusoidal fit is terrible. The best frequency is still 0.105 (period=10) yet the lomb-scargle model fitter overlays something that appears to have a period of 60 days which is wrong?

Could I get clarification on if i'm doing something wrong? I've been told the periodogram is the de facto tool for unenvenly sampled data like this yet... the results seem awful half the time.

To clarify my questions

1. Can it be explained how in the first plot the best-fit frequncy of 0.105 that I feed into the astropy lomb-scargle model fitter somehow creates a correctly matching sinusoid with frequency 0.05? What is the explanation?

2. Why are there 5 strong peaks in the top right plot for the periodogram in example 1 when I only expect 1? The middle two are close to the real value of 0.05 (at 0.045 and 0.055)

Here is the short code used to simulate and plot the data

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as sc
from astropy.stats import LombScargle
import math


#simulate parameters
data_range = [i for i in range(1,1001)]
number_of_samples = 50
gauss_skew = sc.skewnorm.pdf
skew = -10
period = 200
location = data_range[int(len(data_range)/2)]


y1= [(2 * (1. + np.sin(2. * np. pi * x/period)) + np.random.normal(loc =0.0, scale = 0.5)) for x in data_range]
errors1 = [np.random.normal(loc = 0.0, scale = 1) for x in data_range]
y2 = [(1000* (gauss_skew(x,skew,loc=location ,scale = 50)) + np.random.normal(loc =0.0, scale = 1)) for x in data_range]
errors2 = [np.random.normal(loc = 0.0, scale = 1) for x in data_range]

sample_rate = int(len(data_range)/number_of_samples)# To thin the data a bit
y1 = y1[::sample_rate]
y2 = y2[::sample_rate]
errors1 = errors1[::sample_rate]
errors2 = errors2[::sample_rate]
x1= x2 =data_range[::sample_rate]

truth1 = [2* (1. + np.sin(2. * np. pi * x/period)) for x in data_range]
truth2 = [1000 * (gauss_skew(x,-10,loc=location ,scale = 50)) for x in data_range]

truths = [truth1,truth2]
x= [x1,x2]
y=[y1,y2]
errors = [errors1,errors2]

fig,ax  = plt.subplots(nrows=2,ncols=2)

ax[0][0].errorbar(x1,y1,yerr=errors1,fmt='o')
ax[0][0].set_xlabel('time')
ax[1][0].errorbar(x2,y2,yerr=errors2,fmt='o')
ax[1][0].set_xlabel('time')
ax[0][1].set_xlabel('frequency')
ax[1][1].set_xlabel('frequency')


for i in range(0,2):
    frequency, power = LombScargle(x[i],y[i],errors[i]).autopower()
    #Get the best fit frequency as in the tutorial
    best_frequency = frequency[np.argmax(power)]
    print('best frequency:',best_frequency)
    t_fit = np.linspace(x[i][0], math.floor(x[i][-1]),num =50)

    #Fit the best fit frequency
    #plot the best best model based on the best fit
    y_fit = LombScargle(x[i], y[i], errors[i]).model(t_fit, best_frequency)
    ax[i][0].plot(t_fit,y_fit,'g')
    ax[i][0].plot(data_range,truths[i],'r')

    #Plot the PSD
    ax[i][1].plot(frequency,power)
    ax[i][1].axvline(x=best_frequency,color='black', ls='--')
plt.show()

Answer 1

I think I see what is going on. Your best fit is actually a beat frequency between the true frequency (0.005) and double the sampling frequency (0.05). This does indeed produce a model that goes through your data points but because you have only plotted 50 points in the model you haven't been able to see it.

If you change this line t_fit = np.linspace(x[i][0], math.floor(x[i][-1]),num =1000)

so that it shows the model at more points you can see that it is indeed at a much higher frequency.

The trick is to know (roughly) the answer before you start.

If you limit your frequency range in the Lomb-Scargle command to below the Nyquist frequency:

frequency, power = LombScargle(x[i],y[i],errors[i]).autopower(nyquist_factor=1.5)

You will find the behaviour you expected. e.g. (NB I reduced the error bar size to focus in on your problem, but also fixed the error bar size to a single positive value. Your code was introducing very small and even negative errors and I do not know how the code would deal with those).

In unevenly spaced data, things are not so simple because the Nyquist frequency is not well defined.