<p>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.</p>

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

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

<p>The trick is to know (roughly) the answer before you start. </p>

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

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

<p>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 <strong>positive</strong> value. Your code was introducing very small and even negative errors and I do not know how the code would deal with those).</p>

<p><a href="https://i.sstatic.net/k54Xn.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/k54Xn.png" alt="Periodograms"></a> </p>

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

CC BY-SA 4.0
2019-03-12T13:00:12.740Z
Yes, thank you! This is pretty much it. I have been reading more about this in a brilliant periodogram review below

https://arxiv.org/pdf/1703.09824.pdf

and it detailed the effects of the sampling frequency being convolved with the function frequency which I hadn't considered.
wrahman
CC BY-SA 4.0
2019-07-16T16:12:27.830Z
Hey Rob, do you understand what the significance of Lomb-Scargle power is? It's different from the power-spectrum obtained with FFT, and I'm trying to understand what the definition is. I've been asking this question multiple times but got no helpful answer. Please check [here](https://dsp.stackexchange.com/questions/59546/what-is-the-significance-of-y-axis-in-lomb-scargle-periodogram?noredirect=1#comment118913_59546)
Eric Kim
CC BY-SA 4.0
2019-03-12T11:28:52.130Z
2019-03-12T11:28:52.130Z
ProfRob
<p>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.</p>

<p><a href="http://docs.astropy.org/en/stable/stats/lombscargle.html" rel="nofollow noreferrer">http://docs.astropy.org/en/stable/stats/lombscargle.html</a></p>

<p>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. </p>

<p>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 <code>LombScargle.model()</code> on the best fit frequency from <code>frequency[argmax[power]]</code></p>

<p>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.</p>

<p>Example 1
<a href="https://i.sstatic.net/06Bmn.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/06Bmn.png" alt="enter image description here"></a></p>

<p>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? </p>

<p>This does not make sense.</p>

<p>Example 2</p>

<p><a href="https://i.sstatic.net/DJbvQ.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/DJbvQ.png" alt="enter image description here"></a></p>

<p>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?</p>

<p>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.</p>

<p>To clarify my questions</p>

<ol>
<li><p>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?</p></li>
<li><p>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)</p></li>
</ol>

<p>Here is the short code used to simulate and plot the data</p>

<pre><code>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()
</code></pre>

CC BY-SA 4.0
2019-03-12T07:17:05.487Z
Why are you making your error bar have a random size and applying a different error to the actual data? I noted on the plot that some error bars were tiny compared to the others. The data are weighted by the square of the error, so this could have unpredictable consequences.
ProfRob
CC BY-SA 4.0
2019-03-12T09:38:39.340Z
Hi, I am trying to simulate telescope observational data which does not have uniform percentage errors for each data point for observational reasons. At least, the feature of tiny error bars on some data points and large ones on others is something i've seen before in astronomical data, such as the data from the Kaggle LSST Plasticc challenge. 

This is just an initial approximate simulation and I understand it won't be perfect, but by eye, periodicity is clear to me and I thought a well-known tool like the Lomb-Scargle would cope fine with this.
wrahman
CC BY-SA 4.0
2019-03-12T09:41:15.873Z
I'm trying to determine here, if the tool is just failing here or if there is something wrong with my code or that I have a fundamental misunderstanding of the how the tool works. Any help appreciated, thanks!
wrahman
<observational-astronomy><fundamental-astronomy><astropy>
Making sense of the lomb-scargle periodogram
wrahman
# Making sense of the lomb-scargle periodogram

(Asked by: username_0 on 2019-03-12)

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

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?
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()

**Tags:** observational-astronomy, fundamental-astronomy, astropy
> **Comments:**
> **username_1**: Why are you making your error bar have a random size and applying a different error to the actual data? I noted on the plot that some error bars were tiny compared to the others. The data are weighted by the square of the error, so this could have unpredictable consequences.  
> **username_0**: Hi, I am trying to simulate telescope observational data which does not have uniform percentage errors for each data point for observational reasons. At least, the feature of tiny error bars on some data points and large ones on others is something i've seen before in astronomical data, such as the data from the Kaggle LSST Plasticc challenge.   This is just an initial approximate simulation and I understand it won't be perfect, but by eye, periodicity is clear to me and I thought a well-known tool like the Lomb-Scargle would cope fine with this.  
> **username_0**: I'm trying to determine here, if the tool is just failing here or if there is something wrong with my code or that I have a fundamental misunderstanding of the how the tool works. Any help appreciated, thanks!  

---

## 1 answers

### Answer by username_1 on 2019-03-12. Score: 2 | Accepted answer

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).

[![Periodograms](https://i.sstatic.net/k54Xn.png)](https://i.sstatic.net/k54Xn.png)

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

> **Comments:**
> **username_0**: Yes, thank you! This is pretty much it. I have been reading more about this in a brilliant periodogram review below  https://arxiv.org/pdf/1703.09824.pdf  and it detailed the effects of the sampling frequency being convolved with the function frequency which I hadn't considered.  
> **username_2**: Hey Rob, do you understand what the significance of Lomb-Scargle power is? It's different from the power-spectrum obtained with FFT, and I'm trying to understand what the definition is. I've been asking this question multiple times but got no helpful answer. Please check [here](https://dsp.stackexchange.com/questions/59546/what-is-the-significance-of-y-axis-in-lomb-scargle-periodogram?noredirect=1#comment118913_59546)  

---

