Making sense of the 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 cross-posted this to the astronomy stack exchange and I received a satisfactory answer. Those interested can read the below.

https://astronomy.stackexchange.com/questions/29956/making-sense-of-the-lomb-scargle-periodogram

Answer 2

I think the problem lies with you downsampling the data you generate, what you call "Thin the data". By not low pass filtering it first, you generate lots of alias, all the more, since you are downsampling by factor 20!

In your first example, the alias does not seem to be the problem, as the fitter does still yield the correct result, but your second example seems to be messed up by the aliasing.

The factor 20 also explains the apparent difference between the frequency and the resulting sine. $10d \cdot 20=200d$