Peak-Erkennung mit der Lomb-Scargle-Methode

Question

Ich versuche, einen Python-Code arbeiten, der Spitzen in Daten mit der Lomb-Scargle-Methode findet.

http://www.astropython.org/snippets/fast-lomb-scargle-algorithm32/

Mit dieser Methode, wie unten,

import lomb 
x = np.arange(10) 
y = np.sin(x) 
fx,fy, nout, jmax, prob = lomb.fasper(x,y, 6., 6.) 
print jmax 

arbeitet, fein ohne Probleme. Er druckt 8. jedoch auf einem anderen Stück von Daten (Daten-Dump unten),

df = pd.read_csv('extinct.csv',header=None) 
Y = pd.rolling_mean(df[0],window=5) 
fx,fy, nout, jmax, prob = lomb.fasper(np.array(Y.index),np.array(Y),6.,6.) 
print jmax 

zeigt nur 0. Ich habe versucht, verschiedene ofac vorbei, hifac Werte, keine gibt mir sinnvolle Werte.

Hauptfunktion

""" 
from numpy import * 
from numpy.fft import * 

def __spread__(y, yy, n, x, m): 
    """ 
    Given an array yy(0:n-1), extirpolate (spread) a value y into 
    m actual array elements that best approximate the "fictional" 
    (i.e., possible noninteger) array element number x. The weights 
    used are coefficients of the Lagrange interpolating polynomial 
    Arguments: 
    y : 
    yy : 
    n : 
    x : 
    m : 
    Returns: 

    """ 
    nfac=[0,1,1,2,6,24,120,720,5040,40320,362880] 
    if m > 10. : 
    print 'factorial table too small in spread' 
    return 

    ix=long(x) 
    if x == float(ix): 
    yy[ix]=yy[ix]+y 
    else: 
    ilo = long(x-0.5*float(m)+1.0) 
    ilo = min(max(ilo , 1), n-m+1) 
    ihi = ilo+m-1 
    nden = nfac[m] 
    fac=x-ilo 
    for j in range(ilo+1,ihi+1): fac = fac*(x-j) 
    yy[ihi] = yy[ihi] + y*fac/(nden*(x-ihi)) 
    for j in range(ihi-1,ilo-1,-1): 
     nden=(nden/(j+1-ilo))*(j-ihi) 
     yy[j] = yy[j] + y*fac/(nden*(x-j)) 

def fasper(x,y,ofac,hifac, MACC=4): 
    """ function fasper 
    Given abscissas x (which need not be equally spaced) and ordinates 
    y, and given a desired oversampling factor ofac (a typical value 
    being 4 or larger). this routine creates an array wk1 with a 
    sequence of nout increasing frequencies (not angular frequencies) 
    up to hifac times the "average" Nyquist frequency, and creates 
    an array wk2 with the values of the Lomb normalized periodogram at 
    those frequencies. The arrays x and y are not altered. This 
    routine also returns jmax such that wk2(jmax) is the maximum 
    element in wk2, and prob, an estimate of the significance of that 
    maximum against the hypothesis of random noise. A small value of prob 
    indicates that a significant periodic signal is present. 

    Reference: 
    Press, W. H. & Rybicki, G. B. 1989 
    ApJ vol. 338, p. 277-280. 
    Fast algorithm for spectral analysis of unevenly sampled data 
    (1989ApJ...338..277P) 

    Arguments: 
     X : Abscissas array, (e.g. an array of times). 
     Y : Ordinates array, (e.g. corresponding counts). 
     Ofac : Oversampling factor. 
     Hifac : Hifac * "average" Nyquist frequency = highest frequency 
      for which values of the Lomb normalized periodogram will 
      be calculated. 

    Returns: 
     Wk1 : An array of Lomb periodogram frequencies. 
     Wk2 : An array of corresponding values of the Lomb periodogram. 
     Nout : Wk1 & Wk2 dimensions (number of calculated frequencies) 
     Jmax : The array index corresponding to the MAX(Wk2). 
     Prob : False Alarm Probability of the largest Periodogram value 
     MACC : Number of interpolation points per 1/4 cycle 
      of highest frequency 

    History: 
    02/23/2009, v1.0, MF 
     Translation of IDL code (orig. Numerical recipies) 
    """ 
    #Check dimensions of input arrays 
    n = long(len(x)) 
    if n != len(y): 
    print 'Incompatible arrays.' 
    return 

    nout = 0.5*ofac*hifac*n 
    nfreqt = long(ofac*hifac*n*MACC) #Size the FFT as next power 
    nfreq = 64L    # of 2 above nfreqt. 

    while nfreq < nfreqt: 
    nfreq = 2*nfreq 

    ndim = long(2*nfreq) 

    #Compute the mean, variance 
    ave = y.mean() 
    ##sample variance because the divisor is N-1 
    var = ((y-y.mean())**2).sum()/(len(y)-1) 
    # and range of the data. 
    xmin = x.min() 
    xmax = x.max() 
    xdif = xmax-xmin 

    #extirpolate the data into the workspaces 
    wk1 = zeros(ndim, dtype='complex') 
    wk2 = zeros(ndim, dtype='complex') 

    fac = ndim/(xdif*ofac) 
    fndim = ndim 
    ck = ((x-xmin)*fac) % fndim 
    ckk = (2.0*ck) % fndim 

    for j in range(0L, n): 
    __spread__(y[j]-ave,wk1,ndim,ck[j],MACC) 
    __spread__(1.0,wk2,ndim,ckk[j],MACC) 

    #Take the Fast Fourier Transforms 
    wk1 = ifft(wk1)*len(wk1) 
    wk2 = ifft(wk2)*len(wk1) 

    wk1 = wk1[1:nout+1] 
    wk2 = wk2[1:nout+1] 
    rwk1 = wk1.real 
    iwk1 = wk1.imag 
    rwk2 = wk2.real 
    iwk2 = wk2.imag 

    df = 1.0/(xdif*ofac) 

    #Compute the Lomb value for each frequency 
    hypo2 = 2.0 * abs(wk2) 
    hc2wt = rwk2/hypo2 
    hs2wt = iwk2/hypo2 

    cwt = sqrt(0.5+hc2wt) 
    swt = sign(hs2wt)*(sqrt(0.5-hc2wt)) 
    den = 0.5*n+hc2wt*rwk2+hs2wt*iwk2 
    cterm = (cwt*rwk1+swt*iwk1)**2./den 
    sterm = (cwt*iwk1-swt*rwk1)**2./(n-den) 

    wk1 = df*(arange(nout, dtype='float')+1.) 
    wk2 = (cterm+sterm)/(2.0*var) 
    pmax = wk2.max() 
    jmax = wk2.argmax() 


    #Significance estimation 
    #expy = exp(-wk2)   
    #effm = 2.0*(nout)/ofac  
    #sig = effm*expy 
    #ind = (sig > 0.01).nonzero() 
    #sig[ind] = 1.0-(1.0-expy[ind])**effm 

    #Estimate significance of largest peak value 
    expy = exp(-pmax)   
    effm = 2.0*(nout)/ofac  
    prob = effm*expy 

    if prob > 0.01: 
    prob = 1.0-(1.0-expy)**effm 

    return wk1,wk2,nout,jmax,prob 

def getSignificance(wk1, wk2, nout, ofac): 
    """ returns the peak false alarm probabilities 
    Hence the lower is the probability and the more significant is the peak 
    """ 
    expy = exp(-wk2)   
    effm = 2.0*(nout)/ofac  
    sig = effm*expy 
    ind = (sig > 0.01).nonzero() 
    sig[ind] = 1.0-(1.0-expy[ind])**effm 
    return sig 

Daten,

13.5945121951 
13.5945121951 
12.6615853659 
12.6615853659 
12.6615853659 
4.10975609756 
4.10975609756 
4.10975609756 
7.99695121951 
7.99695121951 
16.237804878 
16.237804878 
16.237804878 
16.0823170732 
16.237804878 
16.237804878 
8.92987804878 
8.92987804878 
10.6402439024 
10.6402439024 
28.0548780488 
28.0548780488 
28.0548780488 
27.8993902439 
27.8993902439 
41.5823170732 
41.5823170732 
41.5823170732 
41.5823170732 
41.5823170732 
41.5823170732 
18.7256097561 
15.9268292683 
15.9268292683 
15.9268292683 
15.9268292683 
15.9268292683 
15.9268292683 
14.0609756098 
14.0609756098 
14.0609756098 
14.0609756098 
14.0609756098 
23.8567073171 
23.8567073171 
23.8567073171 
23.8567073171 
25.4115853659 
25.4115853659 
28.0548780488 
40.0274390244 
40.0274390244 
40.0274390244 
40.0274390244 
40.0274390244 
40.0274390244 
20.5914634146 
20.5914634146 
20.4359756098 
19.6585365854 
18.2591463415 
19.3475609756 
18.2591463415 
10.3292682927 
27.743902439 
27.743902439 
27.743902439 
27.743902439 
27.743902439 
27.743902439 
22.3018292683 
22.3018292683 
21.368902439 
21.368902439 
21.368902439 
21.5243902439 
20.4359756098 
20.4359756098 
20.4359756098 
20.4359756098 
20.4359756098 
20.4359756098 
20.4359756098 
11.8841463415 
11.8841463415 
1.0 
11.1067073171 
10.1737804878 
14.5274390244 
14.5274390244 
14.5274390244 
14.5274390244 
14.5274390244 
14.5274390244 
11.7286585366 
11.7286585366 
12.6615853659 
11.7286585366 
8.15243902439 
1.0 
7.84146341463 
6.90853658537 
12.6615853659 
12.6615853659 
12.6615853659 
12.6615853659 
12.6615853659 
12.6615853659 
12.6615853659 
12.6615853659 
12.6615853659 
13.1280487805 
12.9725609756 
12.9725609756 
12.9725609756 
10.3292682927 
10.3292682927 
10.3292682927 
10.3292682927 
9.55182926829 
10.4847560976 
29.9207317073 
29.9207317073 
29.9207317073 
29.9207317073 
30.0762195122 
30.0762195122 
26.1890243902 
7.99695121951 
25.256097561 
7.99695121951 
7.99695121951 
7.99695121951 
6.59756097561 
6.59756097561 
6.59756097561 
6.59756097561 
7.53048780488 
7.53048780488 
7.53048780488 
7.53048780488 
7.53048780488 
7.53048780488 
7.53048780488 
7.53048780488 
10.0182926829 
10.0182926829 
10.0182926829 
10.0182926829 
10.0182926829 
10.0182926829 
10.4847560976 
15.9268292683 
15.9268292683 
15.9268292683 
15.9268292683 
15.9268292683 
16.8597560976 
15.9268292683 
15.9268292683 
16.8597560976 
16.7042682927 
16.7042682927 
16.7042682927 
9.08536585366 
8.46341463415 
8.46341463415 
8.46341463415 
8.46341463415 
6.90853658537 
7.84146341463 
6.90853658537 
4.26524390244 
12.3506097561 
12.3506097561 
12.3506097561 
12.3506097561 
12.3506097561 
12.3506097561 
12.3506097561 
12.3506097561 
12.3506097561 
12.3506097561 
12.3506097561 
14.2164634146 
14.2164634146 
14.2164634146 
14.0609756098 
14.0609756098 
14.0609756098 
14.0609756098 
16.8597560976 
16.8597560976 
16.7042682927 
16.7042682927 
16.7042682927 
16.7042682927 
17.9481707317 
17.9481707317 
19.6585365854 
19.6585365854 
19.6585365854 
19.6585365854 
10.7957317073 
10.7957317073 
10.7957317073 
10.7957317073 
10.7957317073 
12.1951219512 
12.1951219512 
22.9237804878 
22.9237804878 
22.9237804878 
22.9237804878 
22.9237804878 
22.9237804878 
22.9237804878 
7.84146341463 
7.84146341463 
7.84146341463 
7.84146341463 
8.7743902439 
8.7743902439 
7.84146341463 
8.61890243902 
8.61890243902 
8.61890243902 
8.61890243902 
18.2591463415 
18.2591463415 
18.2591463415 
18.2591463415 
18.2591463415 
18.2591463415 
18.2591463415 
18.2591463415 
18.2591463415 
9.39634146341 
9.39634146341 
9.24085365854 
9.24085365854 
9.24085365854 
9.24085365854 
9.08536585366 
9.08536585366 
9.08536585366 
9.08536585366 
9.55182926829 
9.55182926829 
9.55182926829 
9.55182926829 
9.55182926829 
16.5487804878 
16.5487804878 
16.5487804878 
16.5487804878 
16.5487804878 
16.5487804878 
16.5487804878 
16.5487804878 
16.5487804878 
16.5487804878 
16.5487804878 
16.5487804878 
16.5487804878 
16.5487804878 
1.0 
16.0823170732 
16.0823170732 
16.0823170732 
16.0823170732 
16.0823170732 
16.0823170732 
16.0823170732 
16.0823170732 
16.0823170732 
17.1707317073 
17.0152439024 
21.9908536585 
21.9908536585 
21.9908536585 
21.9908536585 
21.9908536585 
21.9908536585 
21.9908536585 
7.84146341463 
8.7743902439 
7.84146341463 
6.75304878049 
5.9756097561 
5.9756097561 
5.9756097561 
5.9756097561 
5.9756097561 
5.9756097561 
3.95426829268 
7.06402439024 
7.06402439024 
7.06402439024 
11.262195122 
11.262195122 
11.262195122 
11.262195122 
11.262195122 
11.262195122 
9.08536585366 
9.86280487805 
7.99695121951 
7.99695121951 
14.2164634146 
14.0609756098 
14.0609756098 
14.0609756098 
14.0609756098 
14.0609756098 
2.24390243902 
2.08841463415 
3.02134146341 
3.02134146341 
2.08841463415 
4.73170731707 
4.73170731707 
4.73170731707 
4.73170731707 
6.44207317073 
6.44207317073 
6.44207317073 
6.44207317073 
6.44207317073 
6.44207317073 
6.44207317073 
6.44207317073 
6.44207317073 
6.44207317073 
6.59756097561 
6.59756097561 
6.59756097561 
6.75304878049 
1.0 
6.28658536585 
6.28658536585 
7.21951219512 
6.28658536585 
10.6402439024 
10.6402439024 
10.6402439024 
10.6402439024 
10.6402439024 
10.6402439024 
10.6402439024 
14.3719512195 
14.3719512195 
15.6158536585 
15.6158536585 
15.6158536585 
35.6737804878 
35.6737804878 
35.6737804878 
35.6737804878 
35.6737804878 
35.6737804878 
35.6737804878 
35.6737804878 
35.6737804878 
35.6737804878 
35.6737804878 
28.6768292683 
28.6768292683 
28.6768292683 
28.6768292683 
28.6768292683 
51.8445121951 
51.8445121951 
51.8445121951 
51.8445121951 
51.8445121951 
52.0 
52.0 
4.42073170732 
4.42073170732 
5.9756097561 
5.9756097561 
5.9756097561 
5.9756097561 
5.9756097561 
5.9756097561 
4.10975609756 
3.95426829268 
3.64329268293 
3.64329268293 
4.73170731707 
4.73170731707 
6.28658536585 
6.28658536585 
6.28658536585 
6.28658536585 
6.28658536585 
6.28658536585 
6.28658536585 
5.9756097561 
5.82012195122 
5.82012195122 
5.82012195122 
5.82012195122 
5.82012195122 
12.1951219512 
12.1951219512 
12.1951219512 
12.1951219512 
12.1951219512 
12.1951219512 
12.1951219512 
12.1951219512 
1.0 
11.7286585366 
11.7286585366 
11.7286585366 
11.7286585366 
11.7286585366 
11.7286585366 
11.1067073171 
11.1067073171 
11.1067073171 
11.1067073171 
11.1067073171 
11.1067073171 
11.1067073171 
11.1067073171 
10.0182926829 
10.0182926829 
16.7042682927 
16.7042682927 
16.7042682927 
16.7042682927 
16.7042682927 
16.7042682927 
29.1432926829 
29.1432926829 
29.1432926829 
29.1432926829 
29.1432926829 
29.1432926829 
29.1432926829 
29.1432926829 
29.1432926829 
1.15548780488 
2.71036585366 
2.71036585366 
2.71036585366 
2.71036585366 
2.71036585366 
2.71036585366 
2.71036585366 
3.17682926829 
4.10975609756 
4.10975609756 
5.9756097561 
5.9756097561 
5.9756097561 
6.90853658537 
5.9756097561 
10.1737804878 
10.1737804878 
10.1737804878 
8.61890243902 
8.46341463415 
8.46341463415 
9.39634146341 
8.46341463415 
8.46341463415 
5.35365853659 
5.35365853659 
5.35365853659 
5.35365853659 
5.35365853659 
5.35365853659 
3.33231707317 
4.42073170732 
3.33231707317 
6.59756097561 
6.44207317073 
5.82012195122 
6.75304878049 
5.82012195122 
5.82012195122 
5.82012195122 
4.73170731707 
5.66463414634 
5.66463414634 
4.73170731707 
4.73170731707 
5.66463414634 
5.66463414634 
5.50914634146 
2.71036585366 
5.50914634146 
2.71036585366 
2.71036585366 
5.50914634146 
5.50914634146 
5.50914634146 
6.28658536585 
6.28658536585 
5.9756097561 
5.9756097561 
7.06402439024 
5.9756097561 
7.53048780488 
8.46341463415 
8.46341463415 
13.2835365854 
13.2835365854 
13.2835365854 
13.2835365854 
2.55487804878 
2.55487804878 
2.55487804878 
2.55487804878 
4.10975609756 
3.17682926829 
3.17682926829 
4.26524390244 
3.64329268293 
3.64329268293 
3.64329268293 
3.33231707317 
3.33231707317 
3.33231707317 
2.24390243902 
3.33231707317 
2.24390243902 
2.24390243902 
3.64329268293 
3.64329268293 
3.64329268293 
3.64329268293 
3.64329268293 
3.64329268293 
7.53048780488 
7.53048780488 
7.53048780488 
7.53048780488 
7.53048780488 
7.53048780488 
7.53048780488 
7.53048780488 
7.53048780488 
6.28658536585 
6.28658536585 
7.21951219512 
6.28658536585 
6.28658536585 
6.28658536585 
6.28658536585 
6.28658536585 
6.28658536585 
3.7987804878 
4.73170731707 
3.7987804878 
3.7987804878 
3.7987804878 
3.7987804878 
3.7987804878 
3.7987804878 
4.26524390244 
4.26524390244 
5.19817073171 
5.19817073171 
6.28658536585 
6.28658536585 
6.28658536585 
6.28658536585 
6.28658536585 
6.28658536585 
6.28658536585 
6.28658536585 
7.53048780488 
7.53048780488 
7.53048780488 
7.53048780488 
7.53048780488 
7.53048780488 
3.7987804878 
3.7987804878 
3.95426829268 
3.02134146341 
3.02134146341 
3.02134146341 
1.0 
1.93292682927 
2.55487804878 
2.55487804878 
5.9756097561 
5.9756097561 
5.9756097561 
5.9756097561 
5.9756097561 
5.9756097561 
5.9756097561 
5.9756097561 
5.9756097561 
5.9756097561 
5.9756097561 
6.28658536585 
6.28658536585 
6.28658536585 
6.28658536585 
6.28658536585 
6.28658536585 
16.0823170732 
16.0823170732 
31.3201219512 
31.3201219512 
31.3201219512 
31.3201219512 
31.3201219512 
31.3201219512 
31.3201219512 
31.3201219512 
3.64329268293 
3.64329268293 
4.26524390244 
4.26524390244 
3.7987804878 
4.73170731707 
3.7987804878 
3.7987804878 
2.55487804878 
3.48780487805 
2.55487804878 
2.55487804878 
3.17682926829 
3.17682926829 
3.17682926829 
3.17682926829 
3.17682926829 
3.17682926829 
3.17682926829 
3.17682926829 
3.17682926829 
3.17682926829 
3.17682926829 
3.17682926829 
3.17682926829 
3.17682926829 
3.17682926829 
3.17682926829 
3.17682926829 
3.17682926829 
3.17682926829 
3.17682926829 
3.33231707317 
12.3506097561 
12.3506097561 
12.3506097561 
12.3506097561 
12.3506097561 
12.3506097561 
4.73170731707 
4.73170731707 
4.73170731707 
4.73170731707 
4.73170731707 
4.73170731707 
4.73170731707 
4.73170731707 
2.86585365854 
2.86585365854 
1.46646341463 
1.46646341463 
1.46646341463 
1.46646341463 
1.46646341463 
1.46646341463 
1.62195121951 
1.62195121951 
1.62195121951 
1.77743902439 
1.77743902439 
4.42073170732 
4.42073170732 
4.42073170732 
4.42073170732 
4.42073170732 
4.42073170732 
4.42073170732 
3.95426829268 
3.95426829268 
2.71036585366 
2.71036585366 
2.71036585366 
2.71036585366 
2.71036585366 
1.77743902439 
2.86585365854 
3.02134146341 
2.86585365854 
2.86585365854 
3.17682926829 
3.17682926829 

Plot

Jede Hilfe dankbar,

Answer 1

Afte würde Wenn man etwas graben will, sieht es so aus, als ob die AstroML-Methode die beste ist.

import numpy as np 
from matplotlib import pyplot as plt 
from astroML.time_series import lomb_scargle, search_frequencies 
import pandas as pd 

df = pd.read_csv('extinct.csv',header=None) 
Y = df[0] 

dy = 0.5 + 0.5 * np.random.random(len(df)) 
omega = np.linspace(10, 100, 1000) 

sig = np.array([0.1, 0.01, 0.001]) 
PS, z = lomb_scargle(df.index, Y, dy, omega, generalized=True, significance=sig) 

plt.plot(omega,PS) 
plt.hold(True) 

xlim = (omega[0], omega[-1]) 
for zi, pi in zip(z, sig): 
    plt.plot(xlim, (zi, zi), ':k', lw=1) 
    plt.text(xlim[-1] - 0.001, zi - 0.02, "$%.1g$" % pi, ha='right', va='top') 
    plt.hold(True) 
plt.show() 

die

Signifikanzniveaus sind in dem Diagramm gibt auch gezeigt. Ich verallgemeinerte LS und benutzte keine Glättung.