import numpy as np
from astroML.time_series import lomb_scargle_bootstrap
from gatspy.periodic import LombScargleFast
print np.__version__ #1.9.1
print gatspy.__version__ #0.2.1
print astroML.__version__ #0.3
time=np.array([ 2456627.79138722, 2456706.70756269, 2456707.71879721,
2456710.7089872 , 2456714.6781531 , 2457003.78622566,
2457004.78832308, 2457008.79132022, 2457011.77312267,
2457013.77502331])
rv=np.array([ 30356.93030223, 6381.7927134 , 6121.9466685 , 5884.45430149,
5641.85369914, 11957.91925097, 12068.74753209, 10966.23226051,
11870.72687516, 9647.03955393])
e=np.array([ 80.0980572 , 80.33726674, 80.98308236, 82.29154533,
91.65632576, 87.24849974, 81.73203085, 82.35611381,
82.30579078, 79.94037756])
print 'N={}'.format(time.size) #10
print 'T.max()-T.min(): {:.0f}'.format(time.max()-time.min()) #386
print 'stdev(RV): {:.0f}'.format(rv.std()) #6.9 km/s
print 'mean error: {:.0f}'.format(e.mean()) #83 m/s
model = LombScargleFast().fit(time,rv,e)
model.optimizer.period_range=(1.2,3000)
fbpP,fbpS=model.find_best_periods(n_periods=6, return_scores=True)
periods=np.arange(1.2,3000,.1)
scores=model.score(periods)
autoP,autoS=model.periodogram_auto(nyquist_factor=100)
print ('Do scores returned by model.find_best_periods '
'match those from model.score?')
try:
assert (model.score(fbpP)==fbpS).all()
print ' Yes'
except AssertionError:
print ' No'
print ('Do scores returned by model.periodogram '
'match those from model.score?')
try:
assert (model.score(fbpP)==model.periodogram(fbpP)).all()
print ' Yes'
except AssertionError:
print ' No'
print ('Do the periods returned by model.find_best_periods have ranking as '
'those from model.score?')
try:
assert (model.score(fbpP).argsort()==fbpS.argsort()).all()
print ' Yes'
except AssertionError:
print ' No'
D=lomb_scargle_bootstrap(time, rv, e, 1/periods, generalized=True,
N_bootstraps=1000, random_state=0)
sig1 = np.percentile(D,[99])
D=lomb_scargle_bootstrap(time, rv, e, 1/autoP, generalized=True,
N_bootstraps=1000, random_state=0)
sig1_auto = np.percentile(D,[99])
D=lomb_scargle_bootstrap(time, rv, e, 1/fbpP, generalized=True,
N_bootstraps=1000, random_state=0)
sig1_fbpP = np.percentile(D,[99])
plt.plot(autoP,autoS,'k')
plt.axhline(sig1_auto,color='k',linestyle='--')
plt.plot(periods,scores,'r')
plt.axhline(sig1,color='r',linestyle='--')
plt.axhline(sig1_fbpP,color='g',linestyle='--')
plt.ylim(0,2)
plt.xlim(.9,3000)
plt.semilogx()
for i,x,y in zip(fbpS.size-fbpS.argsort(),fbpP, model.score(fbpP)):
plt.annotate(str(i),(x,y),(x,1.5),
arrowprops=dict(facecolor='black', shrink=0.05,
width=.5,headwidth=.5,frac=.1))
plt.axvline(1.2,color='k',linestyle=':')
plt.xlabel('P (days)')
plt.ylabel('LS Power')
handles=[plt.plot([-1,-1],[-1,-2],'k--',label='99% CI')[0],
plt.plot([-1,-1],[-1,-2],'k',label='Auto P')[0],
plt.plot([-1,-1],[-1,-2],'k',label='Manual P')[0],
plt.plot([-1,-1],[-1,-2],'k:',label='Min. physical P')[0]]
plt.legend(loc='upper left',frameon=0)
plt.subplots_adjust(.11,.12,.96,.96)
"""
I would think that at a bare minimum the relative ordering of the periods
from find_best_period would match those from .score and .periodogram.
Thought not demonstrated here, I'm also somewhat unnerved that the scores
returned by.score are a function of the input periods at least at the 3% level.
Which is higher than I'd expect from something like variations due to machine
precision alone.
This, combined with the issue above, makes be uncomfortable quiting
significance of peaks deriving from both find_best_periods and a manual search
of scores. My initial interest in find_best_periods was that it seemed to offer
an efficient way to account for the frequency gridding seen in peaks located
using say autoS.argsorted()[::-1][:6] alone. In a large sample of stars this
simplistic approach will see long (~300-2000 day) periods pile up at
~950,620, 500, etc. days.
""
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent