Fig 10.14 : spectral window function: not working....

[Krzysztof Suberlak]

Hi!


Using the code to reproduce Astro ML book  Fig. 10.14 http://www.astroml.org/book_figures/chapter10/fig_sampling.html the PSD of sampling times returns an array of inf's.


In particular, the minimal reproducible code for the error: 

import numpy as np
from astroML.time_series import lomb_scargle

t_obs = 100 * np.random.random(40)
y_obs1 = np.sin(np.pi * t_obs / 3)
y_window = np.ones_like(y_obs1)
omega = np.linspace(0, 5, 1001)[1:]
P_window = lomb_scargle(t_obs, y_window, 1, omega, generalized=False, subtract_mean=False)


results in P_window = array([ inf, inf, inf, inf, inf,......])


Why? I expected to be able to reproduce Fig. 10.14 with the online code. Is this a bug ?


Thank you ,

Best wishes
Chris

[Jake Vanderplas]

Hi Chris,

There are some things that have changed in astroML's dependency libraries, and I think this is the result of one of them. I'll take a look at it today,

   Jake

--
You received this message because you are subscribed to the Google Groups "astroML-general" group.
To unsubscribe from this group and stop receiving emails from it, send an email to astroml-general+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Jake Vanderplas]

Yes, there are some changes in how numpy handles corner cases that affect the behavior of the periodogram. Since astroML was released, there's been a really nice implementation of the Lomb-Scargle periodogram added to astropy; you can read the documentation here: http://docs.astropy.org/en/stable/stats/lombscargle.html

I'm planning to remove all the lomb-scargle code from astroML and use this instead. I think the translation from astroML's syntax to astropy's syntax should be pretty clear in the mean-time.

Thanks for pointing out the bug!

   Jake

[Jake Vanderplas]

I just pushed an update to master that removes the astropy Lomb-Scargle implementation and just wraps astropy instead. I think that should fix the bug.

   Jake

[Krzysztof Suberlak]

Hi Jake!

Thank you for looking into it. This fixes the PSD of sampling times issue. 

However, I see that when using the .autopower() , which seems recommended,  when reproducing the Fig. 10.14 I run into a big difference between employing the  .autopower()  vs .power(), where I   specify  the frequency grid as in the figure  (which, from  reading your blog post, may lead to either over- or undersampling, so that we may miss the true power peak). 

The minimum code to reproduce that : 

# Problem :  why autopower makes such terrible mistake  ?

f0 = 1 / 6.0

t = np.linspace(0, 100, 10000)

y = np.sin(2.0 * np.pi *f0 *  t )

# set grid 

f_grid = np.linspace(0, 1, 1001)[1:] 

P_true_grid = LombScargle(t, y, 1).power(f_grid)

# auto grid 

f_auto, P_true_auto = LombScargle(t, y).autopower()   

    

fig, ax= plt.subplots(3,1,figsize=(5,7))

ax[0].plot(f_grid, P_true_grid)

ax[0].text(0.96, 0.92, "True PSD, \nSet frequency grid", ha='right', va='top', transform=ax[0].transAxes)

ax[0].set_ylim(-0.1, 1.1)

ax[0].set_xlabel('$f$')

ax[0].set_ylabel(r'$P_{\rm LS}(f)$')

ax[0].axvline(f0, ls='--', c='red', lw=2)

ax[1].plot(f_auto, P_true_auto)

ax[1].text(0.96, 0.92, "True PSD, \nAuto frequency grid", ha='right', va='top', transform=ax[1].transAxes)

ax[1].set_ylim(-0.1, 1.1)

ax[1].set_xlim(-0.4,)

ax[1].set_xlabel('$f$')

ax[1].set_ylabel(r'$P_{\rm LS}(f)$')

ax[1].axvline(f0, ls='--', c='red', lw=2)

 

ax[2].scatter(t, y,  c='gray', lw=0, s=1, alpha=0.6)

ax[2].text(0.96, 0.92, "Data", ha='right', va='top', transform=ax[2].transAxes)

ax[2].set_ylim(-1.5, 1.8)

ax[2].set_xlabel('$t$')

ax[2].set_ylabel('$y(t)$')

ax[2].set_ylim(-1.5, 1.8)

plt.show()

    

Is there something wrong with using .autopower() in this case ?  Are there unspoken caveats about it? I looked into the documentation and it seemed fine to use it...

Many thanks,

Chris 

To unsubscribe from this group and stop receiving emails from it, send an email to astroml-gener...@googlegroups.com.

[Jake Vanderplas]

More clarification on that: the autofrequency makes some assumptions about the range of frequencies the data is sensitive to, and the user needs to make sure that range is appropriate for the input data. In this case, the input data is on a regular grid, so any frequency that is a multiple of 100 will produce a perfect fit.

If you set the x limit in the auto frequency case to [0,1] like in the hand-chosen grid, you'll see that the results are the same.

   Jake

To unsubscribe from this group and stop receiving emails from it, send an email to astroml-general+unsubscribe@googlegroups.com.