كيف يمكنني استخدام numpy.ربط للقيام الارتباط الذاتي?

StackOverflow https://stackoverflow.com/questions/643699

سؤال

أنا بحاجة إلى القيام السيارات-ارتباط مجموعة من الأرقام ، كما أفهمها هي علاقة المجموعة مع نفسها.

لقد حاولت ذلك باستخدام numpy هو ربط وظيفة, ولكن لا أعتقد ونتيجة لذلك ، كما انها دائما تقريبا يعطي ناقلات حيث أن الرقم الأول هو لا أكبر, كما يجب أن يكون.

لذا, هذا السؤال هو في الحقيقة سؤالين:

  1. ما هو بالضبط numpy.correlate تفعل ؟
  2. كيف يمكنني استخدامه (أو شيء آخر) للقيام السيارات-الارتباط ؟
هل كانت مفيدة؟

المحلول

للإجابة على السؤال الأول ، numpy.correlate(a, v, mode) أداء الالتواء a مع عكس v وإعطاء نتائج قص المحدد الوضعية.على تعريف الالتواء, C(t)=∑ -∞ < أنا < ∞ aأناvt+i حيث∞ < t < ∞, يسمح النتائج من -∞ إلى ∞, ولكن من الواضح أنك لا يمكن تخزين بلا حدود مجموعة طويلة.لذا يجب أن يكون قص و هذا هو الوضع.هناك 3 طرق مختلفة:كامل, نفس, & صالحة:

  • "الكامل" يعود وضع النتائج لكل t حيث كل a و v بعض التداخل.
  • "نفس" وضع بإرجاع نتيجة مع نفس طول أقصر ناقلات (a أو v).
  • "صالح" وضع بإرجاع نتائج فقط عندما a و v تماما تتداخل مع بعضها البعض.على الوثائق بالنسبة numpy.convolve يعطي المزيد من التفاصيل على وسائط.

بالنسبة لسؤالك الثاني اعتقد numpy.correlate هو مما يتيح لك الارتباط الذاتي ، هو فقط مما يتيح لك أكثر قليلا كذلك.إن الارتباط الذاتي يستخدم لمعرفة كيف مماثلة إشارة ، أو وظيفة ، هو نفسه في وقت معين الفرق.في وقت فرق 0, السيارات-الارتباط يجب أن تكون أعلى لأن الإشارة هي متطابقة إلى نفسه ، لذا من المتوقع أن العنصر الأول في الارتباط الذاتي نتيجة مجموعة سيكون أعظم.ولكن العلاقة لم تبدأ في وقت فرق من 0.ويبدأ في سلبية فارق التوقيت ، يغلق إلى 0 ثم يذهب إيجابية.التي كنت تتوقع:

الارتباط الذاتي(أ) = ∑ -∞ < أنا < ∞ aأناvt+i حيث 0 <= t < ∞

ولكن ما حصلت عليه هو:

الارتباط الذاتي(أ) = ∑ -∞ < أنا < ∞ aأناvt+i حيث∞ < t < ∞

ما تحتاج إلى القيام به هو أن تأخذ النصف الأخير من ارتباط النتيجة ، وأنه ينبغي أن يكون الارتباط الذاتي كنت تبحث عن.الثعبان بسيط الدالة على ذلك:

def autocorr(x):
    result = numpy.correlate(x, x, mode='full')
    return result[result.size/2:]

سوف تحتاج بطبيعة الحال خطأ التحقق للتأكد من أن x هو في الواقع 1-د الصفيف.كما أن هذا التفسير ربما ليس رياضيا دقيقا.لقد تم رمي حول علامات ما لا نهاية لأن تعريف الالتواء يستخدم لهم, ولكن هذا لا ينطبق بالضرورة على الارتباط الذاتي.لذا الجزء النظري من هذا التفسير قد يكون قليلا متزعزع ، ولكن نأمل أن نتائج عملية مفيدة. هذه صفحات على الارتباط الذاتي هي مفيدة جدا و يمكن أن تعطيك أفضل بكثير الخلفية النظرية إذا كنت لا تمانع في الخوض في التدوين و الثقيلة المفاهيم.

نصائح أخرى

السيارات-الارتباط يأتي في نسختين:الإحصائية و الإلتواء.كلاهما تفعل نفس الشيء, ما عدا بعض التفاصيل:الإحصائية الإصدار هو تطبيع على الفترة [-1,1].هنا هو مثال على كيفية عمل إحصائية واحدة:

def acf(x, length=20):
    return numpy.array([1]+[numpy.corrcoef(x[:-i], x[i:])[0,1]  \
        for i in range(1, length)])

استخدام numpy.corrcoef وظيفة بدلا من numpy.correlate لحساب علاقة إحصائية عن تأخر t:

def autocorr(x, t=1):
    return numpy.corrcoef(numpy.array([x[:-t], x[t:]]))

كما أنا فقط وقعت في نفس المشكلة ، أود أن حصة بضعة أسطر من التعليمات البرمجية معك.في الواقع هناك عدة بدلا وظائف مماثلة عن الارتباط الذاتي في ستاكوفيرفلوو الآن.إذا كان تعريف الارتباط الذاتي كما a(x, L) = sum(k=0,N-L-1)((xk-xbar)*(x(k+L)-xbar))/sum(k=0,N-1)((xk-xbar)**2) [هذا هو التعريف الوارد في IDL a_correlate وظيفة وهو يتفق مع ما أرى في الجواب 2 من السؤال #12269834] ، ثم يلي ويبدو أن تعطي النتائج الصحيحة:

import numpy as np
import matplotlib.pyplot as plt

# generate some data
x = np.arange(0.,6.12,0.01)
y = np.sin(x)
# y = np.random.uniform(size=300)
yunbiased = y-np.mean(y)
ynorm = np.sum(yunbiased**2)
acor = np.correlate(yunbiased, yunbiased, "same")/ynorm
# use only second half
acor = acor[len(acor)/2:]

plt.plot(acor)
plt.show()

كما ترون لقد اختبرت هذا مع الخطيئة منحنى موحد توزيع عشوائي كل النتائج تبدو وكأنها لا يمكن أن نتوقع منهم.علما بأنني استخدمت mode="same" بدلا من mode="full" كما فعل غيره.

أعتقد أن هناك 2 من الأشياء التي تضيف الارتباك على هذا الموضوع:

  1. إحصائية v. s.معالجة الإشارات التعريف:كما أشار آخرون ، في الإحصاءات نحن تطبيع السيارات-الارتباط إلى [-1,1].
  2. جزئية v. s.غير جزئية يعني/الفرق:عندما timeseries التحولات في تأخر>0 ، التداخل الحجم دائما < الطول الأصلي.لا يعني std الأصلي (غير جزئية) أو دائما حساب جديد يعني و الأمراض المنقولة جنسيا باستخدام المتغيرة التداخل (جزئي) فرقا.(ربما هناك الرسمية المصطلح هذا ، ولكن سأستخدم "جزئية" في الوقت الحالي).

لقد خلق 5 وظائف حساب السيارات-العلاقة 1d مجموعة جزئية v. s.غير الفروق الجزئية.بعض استخدام صيغة من إحصاءات استخدام بعض ترتبط في معالجة الإشارات المعنى ، والتي يمكن أيضا أن يتم ذلك عن طريق الاتحاد الفرنسي للتنس.ولكن كل النتائج السيارات-الارتباطات في الإحصاءات الوضوح ، بحيث توضح كيف ترتبط مع بعضها البعض.رمز أدناه:

import numpy
import matplotlib.pyplot as plt

def autocorr1(x,lags):
    '''numpy.corrcoef, partial'''

    corr=[1. if l==0 else numpy.corrcoef(x[l:],x[:-l])[0][1] for l in lags]
    return numpy.array(corr)

def autocorr2(x,lags):
    '''manualy compute, non partial'''

    mean=numpy.mean(x)
    var=numpy.var(x)
    xp=x-mean
    corr=[1. if l==0 else numpy.sum(xp[l:]*xp[:-l])/len(x)/var for l in lags]

    return numpy.array(corr)

def autocorr3(x,lags):
    '''fft, pad 0s, non partial'''

    n=len(x)
    # pad 0s to 2n-1
    ext_size=2*n-1
    # nearest power of 2
    fsize=2**numpy.ceil(numpy.log2(ext_size)).astype('int')

    xp=x-numpy.mean(x)
    var=numpy.var(x)

    # do fft and ifft
    cf=numpy.fft.fft(xp,fsize)
    sf=cf.conjugate()*cf
    corr=numpy.fft.ifft(sf).real
    corr=corr/var/n

    return corr[:len(lags)]

def autocorr4(x,lags):
    '''fft, don't pad 0s, non partial'''
    mean=x.mean()
    var=numpy.var(x)
    xp=x-mean

    cf=numpy.fft.fft(xp)
    sf=cf.conjugate()*cf
    corr=numpy.fft.ifft(sf).real/var/len(x)

    return corr[:len(lags)]

def autocorr5(x,lags):
    '''numpy.correlate, non partial'''
    mean=x.mean()
    var=numpy.var(x)
    xp=x-mean
    corr=numpy.correlate(xp,xp,'full')[len(x)-1:]/var/len(x)

    return corr[:len(lags)]


if __name__=='__main__':

    y=[28,28,26,19,16,24,26,24,24,29,29,27,31,26,38,23,13,14,28,19,19,\
            17,22,2,4,5,7,8,14,14,23]
    y=numpy.array(y).astype('float')

    lags=range(15)
    fig,ax=plt.subplots()

    for funcii, labelii in zip([autocorr1, autocorr2, autocorr3, autocorr4,
        autocorr5], ['np.corrcoef, partial', 'manual, non-partial',
            'fft, pad 0s, non-partial', 'fft, no padding, non-partial',
            'np.correlate, non-partial']):

        cii=funcii(y,lags)
        print(labelii)
        print(cii)
        ax.plot(lags,cii,label=labelii)

    ax.set_xlabel('lag')
    ax.set_ylabel('correlation coefficient')
    ax.legend()
    plt.show()

هنا هو إخراج الشكل:

enter image description here

نحن لا نرى كل 5 خطوط لأن 3 منهم التداخل (في الأرجواني).التداخل كلها غير جزئية السيارات-الارتباطات.وذلك لأن حسابات من إشارة أساليب المعالجة (np.correlate, FFT) لا حساب مختلف يعني/std لكل التداخل.

نلاحظ أيضا أن fft, no padding, non-partial (الخط الأحمر) نتيجة مختلفة ، لأنه لم وسادة timeseries مع 0s قبل الاتحاد الفرنسي للتنس ، لذلك التعميم الاتحاد الفرنسي للتنس.لا أستطيع أن أشرح بالتفصيل لماذا, هذا ما تعلمته من أماكن أخرى.

السؤال 1 وقد تم بالفعل مناقشة مستفيضة في عدة إجابات ممتازة هنا.

ظننت أن أطلعكم على بضعة أسطر من التعليمات البرمجية التي تسمح لك لحساب الارتباط الذاتي من إشارة تستند فقط على الخصائص الرياضية من الارتباط الذاتي.أي أن الارتباط الذاتي قد يكون حسابها بالطريقة التالية:

  1. طرح يعني من إشارة على غير متحيزة إشارة

  2. حساب تحويل فورييه منحازة إشارة

  3. حساب الكثافة الطيفية من إشارة طريق اتخاذ مربع القاعدة من كل قيمة تحويل فورييه منحازة إشارة

  4. حساب معكوس تحويل فورييه الكثافة الطيفية

  5. تطبيع معكوس تحويل فورييه الكثافة الطيفية من مجموع المربعات من غير منحازة إشارة تأخذ فقط نصف الناتج ناقلات

رمز للقيام بذلك هي التالية:

def autocorrelation (x) :
    """
    Compute the autocorrelation of the signal, based on the properties of the
    power spectral density of the signal.
    """
    xp = x-np.mean(x)
    f = np.fft.fft(xp)
    p = np.array([np.real(v)**2+np.imag(v)**2 for v in f])
    pi = np.fft.ifft(p)
    return np.real(pi)[:x.size/2]/np.sum(xp**2)

يمكنني استخدام طالب.CORREL عن الارتباط الذاتي مثل هذا, وأظن أنك يمكن أن تفعل الشيء نفسه مع الحزم الأخرى:

def autocorrelate(x, period):

    # x is a deep indicator array 
    # period of sample and slices of comparison

    # oldest data (period of input array) may be nan; remove it
    x = x[-np.count_nonzero(~np.isnan(x)):]
    # subtract mean to normalize indicator
    x -= np.mean(x)
    # isolate the recent sample to be autocorrelated
    sample = x[-period:]
    # create slices of indicator data
    correls = []
    for n in range((len(x)-1), period, -1):
        alpha = period + n
        slices = (x[-alpha:])[:period]
        # compare each slice to the recent sample
        correls.append(ta.CORREL(slices, sample, period)[-1])
    # fill in zeros for sample overlap period of recent correlations    
    for n in range(period,0,-1):
        correls.append(0)
    # oldest data (autocorrelation period) will be nan; remove it
    correls = np.array(correls[-np.count_nonzero(~np.isnan(correls)):])      

    return correls

# CORRELATION OF BEST FIT
# the highest value correlation    
max_value = np.max(correls)
# index of the best correlation
max_index = np.argmax(correls)

باستخدام تحويل فورييه و الإلتواء نظرية

الوقت هو complexicity N*log(N)

def autocorr1(x):
    r2=np.fft.ifft(np.abs(np.fft.fft(x))**2).real
    return r2[:len(x)//2]

هنا هو تطبيع وغير متحيزة الإصدار هو أيضا N*log(N)

def autocorr2(x):
    r2=np.fft.ifft(np.abs(np.fft.fft(x))**2).real
    c=(r2/x.shape-np.mean(x)**2)/np.std(x)**2
    return c[:len(x)//2]

الطريقة التي تقدمها A.ليفي يعمل ، ولكن اختبرته في جهاز الكمبيوتر الخاص بي ، وقتها complexicity يبدو أن N*N

def autocorr(x):
    result = numpy.correlate(x, x, mode='full')
    return result[result.size/2:]

أعتقد أن الجواب الحقيقي إلى OP السؤال بإيجاز الواردة في هذا مقتطف من Numpy.ربط الوثائق:

mode : {'valid', 'same', 'full'}, optional
    Refer to the `convolve` docstring.  Note that the default
    is `valid`, unlike `convolve`, which uses `full`.

وهذا يعني أن, عندما تستخدم مع أي وضع تعريف ، Numpy.ترتبط وظيفة سيعود العددية ، عندما تعطى نفس متجه اثنين من وسائط الإدخال (أي- عندما تستخدم لأداء الارتباط الذاتي).

حل بسيط دون الباندا:

import numpy as np

def auto_corrcoef(x):
   return np.corrcoef(x[1:-1], x[2:])[0,1]

مؤامرة الإحصائية الارتباط الذاتي نظرا الباندا datatime سلسلة من العوائد:

import matplotlib.pyplot as plt

def plot_autocorr(returns, lags):
    autocorrelation = []
    for lag in range(lags+1):
        corr_lag = returns.corr(returns.shift(-lag)) 
        autocorrelation.append(corr_lag)
    plt.plot(range(lags+1), autocorrelation, '--o')
    plt.xticks(range(lags+1))
    return np.array(autocorrelation)

أنا في علم الأحياء الحسابي ، وعندما كان علي حساب السيارات/عبر الارتباطات بين الأزواج من الوقت سلسلة من العمليات العشوائية أدركت أن np.correlate كان لا يقوم بهذه المهمة احتاجه.

والواقع أن ما يبدو أن يكون في عداد المفقودين من np.correlate هو حيث بلغ متوسطها على جميع الأزواج من وقت نقاط في المسافة 𝜏.

هنا هو كيف يمكنني تعريف وظيفة تفعل ما يلزم:

def autocross(x, y):
    c = np.correlate(x, y, "same")
    v = [c[i]/( len(x)-abs( i - (len(x)/2)  ) ) for i in range(len(c))]
    return v

يبدو لي أيا من الإجابات السابقة تغطي هذه الحالة السيارات/عبر الارتباط:نأمل أن يكون هذا الجواب قد يكون من المفيد أن يعمل على العمليات العشوائية مثلي.

بديل numpy.ربط متاح في statsmodels.tsa.stattools.acf().هذا ينتج باستمرار تناقص الارتباط الذاتي وظيفة مثل الموصوفة من قبل المرجع.تنفيذ الأمر بسيط إلى حد ما:

from statsmodels.tsa import stattools
# x = 1-D array
# Yield normalized autocorrelation function of number lags
autocorr = stattools.acf( x )

# Get autocorrelation coefficient at lag = 1
autocorr_coeff = autocorr[1]

السلوك الافتراضي هو التوقف عند 40 nlags ، ولكن هذا يمكن تعديلها مع nlag= الخيار للتطبيق الخاص بك محددة.هناك الاقتباس في أسفل الصفحة إحصاءات وراء وظيفة.

مرخصة بموجب: CC-BY-SA مع الإسناد
لا تنتمي إلى StackOverflow
scroll top