الگوریتم محاسبه واریانس

از ویکی‌پدیا، دانشنامهٔ آزاد
پرش به: ناوبری، جستجو

الگوریتم Naïve[ویرایش]

فرمول محاسبه واریانس جامعه آماری با اندازه N به این صورت است:

\sigma^2 = \displaystyle\frac {\sum_{i=1}^N x_i^2 - (\sum_{i=1}^N x_i)^2/N}{N}. \!

فرمول محاسبه تخمین unbiased از واریانس جامعه آماری از یک نمونه محدود از n مشاهده به این صورت میباشد:

s^2 = \displaystyle\frac {\sum_{i=1}^n x_i^2 - (\sum_{i=1}^n x_i)^2/n}{n-1}. \!

بنابراین الگوریتم naive برای محاسبه واریانس تخمینی به صورت زیر میباشد:

def naive_variance(data):
    n = 0
    Sum = 0
    Sum_sqr = 0
 
    for x in data:
        n = n + 1
        Sum = Sum + x
        Sum_sqr = Sum_sqr + x*x
 
    variance = (Sum_sqr - (Sum*Sum)/n)/(n - 1)
    return variance

این الگوریتم به راحتی میتواند برای محاسبه واریانس یک جامعه محدود به کار رود: فقط کافیست که در خط آخر به جای تقسیم بر N، بر n − 1 تقسیم گردد. چون Sum_sqr و (Sum*Sum)/n میتوانند مقدار بسیار نزدیکی داشته باشند، تقریب زدن میتواند منجر شود که دقت نتیجه بسیار کمتر از دقت ذاتی موجود در محاسبات ممیز شناور شود. بنابراین این الگوریتم نباید در عمل مورد استفاده قرار گیرد.[۱][۲] این موضوع بخصوص در مواردی که انحراف معیار نسبت به میانگین بسیار کوچک است، تشدید میشود. اگرچه، این الگوریتم میتواند با استفاده از روش میانگین مفروض بهبود یابد.

الگوریتم دومرحله ای[ویرایش]

یک روش جایگزین، استفاده از فرمول دیگری برای محاسبه واریانس میباشد، به این صورت که ابتدا میانگیم نمونه را محاسبه میکند:

\bar x = \displaystyle \frac {\sum_{j=1}^n x_j}{n},

و سپس مجموع مربعات فاصله از میانگین را بدست می آورد:

\mathrm{variance} = s^2 = \displaystyle\frac {\sum_{i=1}^n (x_i - \bar x)^2}{n-1} \!,

S برابر با مقدار انحراف معیار است، که با استفاده از شبه کد زیر به دست می آید:

def two_pass_variance(data):
    n    = 0
    sum1 = 0
    sum2 = 0
 
    for x in data:
        n    = n + 1
        sum1 = sum1 + x
 
    mean = sum1/n
 
    for x in data:
        sum2 = sum2 + (x - mean)*(x - mean)
 
    variance = sum2/(n - 1)
    return variance

این الگوریتم همواره از لحاظ عددی، سازگار است، مگر برای nهای بزرگ.[۱][۳] اگرچه میتواند برای نمونه هایی که در شرایط زیر صدق میکنند، دارای خطا باشد: 1. اکثر داده های آن، بسیار نزدیک به میانگین(و نه برابر با آن) باشند. 2. فقط تعداد اندکی از داده ها، با فاصله ی بسیار زیاد از میانگین باشند.[نیازمند منبع]. نتایج هر دو الگوریتم میتواند بی اندازه، به ترتیب داده ها وابسته باشد و میتواند نتایج ضعیفی برای مجموعه داده های بسیار بزرگ (بخاطر تکرر خطای گرد کردن) داشته باشد. تکنیکهایی همچون مجموعِ جبران شده میتواند تا درجه ای برای از بین بردن این خطا مورد استفاده قرار گیرند.

مغایرت جبران شده[ویرایش]

نسخه ی مجموعِ جبران شده از الگوریتم فوق بدین صورت است:[۴]

def compensated_variance(data):
    n = 0
    sum1 = 0
    for x in data:
        n = n + 1
        sum1 = sum1 + x
    mean = sum1/n
 
    sum2 = 0
    sum3 = 0
    for x in data:
        sum2 = sum2 + (x - mean)**2
        sum3 = sum3 + (x - mean)
    variance = (sum2 - sum3**2/n)/(n - 1)
    return variance

الگوریتم برخط[ویرایش]

محاسبه ی واریانس در یک مرحله، اغلب میتواند بسیار مفید باشد، به عنوان مثال، وقتی داده ها در حال جمع آوری هستند و فضای کافی برای ذخیره همه داده ها نداریم، و یا موقعی که هزینه دسترسی به حافظه، بر هزینه محاسبات تاثیر زیادی دارد. برای چنین الگوریتم برخطی، رابطه رخداد مجدد بین کمیتهایی مورد استفاده در محاسبات، مورد نیاز میباشد تا محاسبات ما به صورت به روشی سازگار انجام شود. فرمولهای زیر برای به روزرسانی میانگین و واریانس(تخمینی) داده ها، در زمان افزودن عنصر جدید x_{\mathrm{new}} استفاده میشود. در اینجا، xn نشان دهنده میانگین n نمونه ی اول (x1, ..., xn) و s2n واریانس نمونه و σ2N واریانس جامعه ی آماری آنهاست.

\bar x_n = \frac{(n-1) \, \bar x_{n-1} + x_n}{n} = \bar x_{n-1} + \frac{x_n - \bar x_{n-1}}{n} \!
s^2_n = \frac{(n-2)}{(n-1)} \, s^2_{n-1} + \frac{(x_n - \bar x_{n-1})^2}{n}, \quad n>1
\sigma^2_N = \frac{(N-1) \, \sigma^2_{N-1} + (x_N - \bar x_{N-1})(x_N - \bar x_{N})}{N}.

همچنین میتوان نشان داد که کمیت مناسبتر برای به روزرسانی مجموع مربعات فاصله از میانگین کنونی (\textstyle\sum_{i=1}^n (x_i - \bar x_n)^2) ، با عنوان M_{2,n} دراینجا نشان داده شده است:

M_{2,n}\! = M_{2,n-1} + (x_n - \bar x_{n-1})(x_n - \bar x_n)
s^2_n = \frac{M_{2,n}}{n-1}
\sigma^2_N = \frac{M_{2,N}}{N}

یک الگوریتم عددی سازگار به صورت زیر داده شده است، که البته میانگین را نیز محاسبه می نماید. این الگوریتم، منتسب به Knuth[۵] است که او نیز از الگوریتمی از Welford[۶] استفاده کرده و به صورت کامل تحلیل نموده است.[۷][۸] همچنین معمول است که متغیرهای زیر به جای همدیگر استفاده شوند: M_k = \bar x_k and S_k = M_{2,k}[۹]

def online_variance(data):
    n = 0
    mean = 0
    M2 = 0
 
    for x in data:
        n = n + 1
        delta = x - mean
        mean = mean + delta/n
        M2 = M2 + delta*(x - mean)
 
    variance = M2/(n - 1)
    return variance

این الگوریتم بخاطر کاهش میزان گرد کردن، میتواند کاهشِ دقتِ کمتری داشته باشد، اما احتمالاً نمی‌تواند خیلی بهینه باشد، چرا که عمل تقسیم درون حلقه انجام میشود. برای الگوریتم عظیمتر دومرحله ای برای محاسبه واریانس، بهتر است ابتدا تخمینی از میانگین را محاسبه کنیم و سپس در ادامه از این الگوریتم استفاده نماییم. الگوریتم موازی که در ادامه خواهد آمد، نشان میدهد که چگونه چندین مجموعه آماری را که به صورت برخط محاسبه میشوند، ادغام نمود.

الگوریتم افزایشی وزن دار[ویرایش]

این الگوریتم میتواند برای کنترل نمونه های دارای وزن تعمیم یابد، بدین صورت که میتوان شمارنده n را برای مجموع وزنهایی مشاهده شده تاکنون استفاده کرد. West (1979) [۱۰] این الگوریتم افزایشی را پیشنهاد میدهد:

def weighted_incremental_variance(dataWeightPairs):
    sumweight = 0
    mean = 0
    M2 = 0
 
    for x, weight in dataWeightPairs:  # Alternatively "for x, weight in zip(data, weights):"
        temp = weight + sumweight
        delta = x - mean
        R = delta * weight / temp
        mean = mean + R
        M2 = M2 + sumweight * delta * R  # Alternatively, "M2 = M2 + weight * delta * (x−mean)"
        sumweight = temp
 
    variance_n = M2/sumweight
    variance = variance_n * len(dataWeightPairs)/(len(dataWeightPairs) - 1)

الگوریتم موازی[ویرایش]

Chan et al.[۴] میگوید که سومین الگوریتم برخط فوق، نمونه ای خاص از الگوریتم است که برای هر تقسیم بندی از نمونه X به مجموعه های X_A, X_B درست عمل میکند:

\delta\! = \bar x_B - \bar x_A
\bar x_X = \bar x_A + \delta\cdot\frac{n_B}{n_X}
M_{2,X} = M_{2,A} + M_{2,B} + \delta^2\cdot\frac{n_A n_B}{n_X}.

این موضوع میتواند زمانی مفید باشد که برای مثال، چندین واحد پرداش برای قسمتهای مختلف از ورودی در نظر گرفته شده باشد. شیوه Chan روشی برای تخمین میانگین، زمانی که n_A \approx n_B و هر دو بسیار بزرگ باشند، از لحاظ عددی نامطمئن است، چراکه خطای عددی موجود در \bar x_B - \bar x_A متناسب با زمانی نیست که n_B = 1. در چنین مواردی، \bar x_X = \frac{n_A \bar x_A + n_B \bar x_B}{n_A + n_B} ترجیح داده میشود.

مثال[ویرایش]

فرض کنید که عملیات ممیز شناور از استاندارد محاسباتی IEEE 754 double-precision استفاده میکند. نمونه ی {4، 7، 13، 16} را از یک جامعه آماری نامحدود فرض کنید. براساس این نمونه، میانگین تخمینی جامعه برابر با 10 است و واریانس تخمینی بدون پیشقدرِ جامعه برابر با 30 است. هر دو الگورتیم اول و دوم این مقادیر را به درستی محاسبه میکند. سپس نمونه { 108 + 4, 108 + 7, 108 + 13, 108 + 16 } را در نظر بگیرید، که مقدار واریانسی برابر با نمونه ی قبلی به ما میدهد. الگوریتم دوم این واریانس تخمینی را به درستی محاسبه میکند، اما الگوریتم اول مقدار 29.33333333 را به جای 30 برمیگرداند. اگرچه این کمبود دقت قابل تحمل است و به عنوان نقص کوچکی در الگوریتم اول مشاهده میشود، اما به سادگی داده هایی یافت که نقص اساسی در الگوریتم naive را آشکار میسازد: نمونه ی { 109 + 4, 109 + 7, 109 + 13, 109 + 16 } را در نظر بگیرید. دوباره واریانس تخمینی برابر با 30 میباشد که توسط الگوریتم دوم به سادگی قابل محاسبه است، اما این دفعه الگوریتم naive مقدار −170.66666666666666 را به عنوان جواب برمیگرداند. این یک مشکل اساسی در الگوریتم اول میباشد و بخاطر تقریب زدن فاجعه آمیزی رخ داده است که در مرحله آخر الگوریتم، هنگام تفریق دو عدد مشابه صورت گرفته است.

آمارهای رده بالاتر[ویرایش]

Terriberry[۱۱] فرمول Chan را برای محاسبه سومین و چهارمین گشتاور مرکزی تعمیم میدهد. برای مثال برای تخمین چولگی و کشیدگی مورد نیاز است:

M_{3,X} = M_{3,A} + M_{3,B} + \delta^3\frac{n_A n_B (n_A - n_B)}{n_X^2} + 3\delta\frac{n_AM_{2,B} - n_BM_{2,A}}{n_X}
\begin{align}
M_{4,X} = M_{4,A} + M_{4,B} & + \delta^4\frac{n_A n_B \left(n_A^2 - n_A n_B + n_B^2\right)}{n_X^3} \\
                    & + 6\delta^2\frac{n_A^2 M_{2,B} + n_B^2 M_{2,A}}{n_X^2} + 4\delta\frac{n_AM_{3,B} - n_BM_{3,A}}{n_X} \\
\end{align}

اینجا M_k برابر با توانهای فاصله با میانه \Sigma(x - \overline{x})^k است که نتیجه میدهد:

چولگی: g_1 = \frac{\sqrt{n} M_3}{M_2^{3/2}},
کشیدگی: g_2 = \frac{n M_4}{M_2^2}-3.

برای نسخه افزایشی (مثلاً B = \{x\})، به صورت زیر ساده میشود:

\delta\! = x - m
m' = m + \frac{\delta}{n}
M_2' = M_2 + \delta^2 \frac{ n-1}{n}
M_3' = M_3 + \delta^3 \frac{ (n - 1) (n - 2)}{n^2} - \frac{3\delta M_2}{n}
M_4' = M_4 + \frac{\delta^4 (n - 1) (n^2 - 3n + 3)}{n^3} + \frac{6\delta^2 M_2}{n^2} - \frac{4\delta M_3}{n}

با نگهداشتن مقدار \delta / n ، فقط یک عمل تقسیم مورد نیاز است و آمار رده بالاتر میتواند بنابراین برای هزینه افزایشی محاسبه شود. مثال الگوریتم برخط برای کشیدگی به صورت زیر پیاده سازی شده است:

def online_kurtosis(data):
    n = 0
    mean = 0
    M2 = 0
    M3 = 0
    M4 = 0
 
    for x in data:
        n1 = n
        n = n + 1
        delta = x - mean
        delta_n = delta / n
        delta_n2 = delta_n * delta_n
        term1 = delta * delta_n * n1
        mean = mean + delta_n
        M4 = M4 + term1 * delta_n2 * (n*n - 3*n + 3) + 6 * delta_n2 * M2 - 4 * delta_n * M3
        M3 = M3 + term1 * delta_n * (n - 2) - 3 * delta_n * M2
        M2 = M2 + term1
 
    kurtosis = (n*M4) / (M2*M2) - 3
    return kurtosis

Pébay [۱۲]سپس این نتایج را برای گشتاورهای مرکزیِ رده-اختیاری، برای موارد افزایشی و دوبه دو تعمیم داد. میتوان فرمولهای مشابهی را برای کواریانس یافت. Choi و Sweetman[۱۳] دو روش جایگزین برای محاسبه چولگی و کشیدگی پیشنهاد میدهند، که هر کدام میتواند مقدار قابل ملاحظه ای در حافظه مورد نیاز و زمان CPU در کاربردهای معینی، صرفه جویی کند. روش اول، برای محاسبه گشتاورهای آماری، با جداسازی داده ها، گشتاورها را از رابطه هندسیِ هیستوگرامِ به دست آمده محاسبه میکند، که به صورت بهینه ای تبدیل به یک الگوریتم یک مرحله ای برای گشتاورهای بالاتر میشود. یکی از مزایای کار این است که محاسبات گشتاورهای آماری میتواند با دقت اختیاری محاسبه شود، پس میتوان این دقت را با دقت فرمت انبار داده ویا سخت افزار اندازه گیری اصلی، منطبق کرد. هیستوگرام نسبی از یک متغیر تصادفی میتواند یک در این راه قراردادی ساخته شود: محدوده مقادیر بالقوه به بازه هایی تقسیم میشود و تعداد رخدادها در هر بازه شمارش و به گونه ای کشیده میشود که مساحت هر مستطیل برابر با سهم رخدادهای آن بازه شود:

 H(x_k)=\frac{h(x_k)}{A}

که h(x_k) و H(x_k) نشاندهنده ی میزان فرکانس و فرکانس نسبی در بازه ی x_k است و A= \sum_{k=1}^{K} h(x_k)
\,\Delta x_k مجموع تمام نواحی هیستوگرام است. بعد از نرمال سازی، n گشتاور خام و گشتاورهای مرکزیِ x(t) میتواند از هیستوگرام نسبی محاسبه شود:


 m_n^{(h)} = \sum_{k=1}^{K}  x_k^n \, H(x_k) \Delta x_k
            = \frac{1}{A} \sum_{k=1}^{K}  x_k^n \, h(x_k) \Delta x_k

 \theta_n^{(h)}= \sum_{k=1}^{K} \Big(x_k-m_1^{(h)}\Big)^n \, H(x_k)\Delta x_k
               = \frac{1}{A} \sum_{k=1}^{K} \Big(x_k-m_1^{(h)}\Big)^n \, h(x_k) \Delta x_k

که توان ^{(h)} نشان دهنده گشتاورهایی است که از هیستوگرام محاسبه شده است. برای بازه ی با ضخامت ثابت \Delta
x_k=\Delta x این دو عبارت میتواند با استفاده از I= A/\Delta x ساده شود:


 m_n^{(h)}= \frac{1}{I} {\sum_{k=1}^{K} x_k^n \, h(x_k)}

 \theta_n^{(h)}= \frac{1}{I}{\sum_{k=1}^{K} \Big(x_k-m_1^{(h)}\Big)^n \, h(x_k)}

دومین روش Choi و Sweetman[۱۳] یک متدولوژی تحلیلی برای ترکیب گشتاورهای آماری از بخشهای مجزا از یک بازه ی زمانی است که گشتاورهای کلی بدست آمده، همان گشتاورهای بازه ی زمانی کامل است. این متدولوژی میتواند برای محاسبه موازی گشتاورهای آماری با کمک ترکیب گشتاورهای پس آیند(بعدی) ویا ترکیب گشتاورهای آماری محاسبه شده در زمانهای پشت سرهم، استفاده شود. اگر تعداد Q مجموعه از گشتاورهای آماری به این صورت نمایش داده شوند: (\gamma_{0,q},\mu_{q},\sigma^2_{q},\alpha_{3,q},\alpha_{4,q})
\quad که در آن q=1,2,...,Q است، سپس هر \gamma_n میتواند براساس n گشتاور اولیه معادل نمایش داده شود:


 \gamma_{n,q}= m_{n,q} \gamma_{0,q} \qquad \quad \textrm{for} \quad n=1,2,3,4  \quad \text{ and } \quad q = 1,2, \dots ,Q

که \gamma_{0,q} به صورت کلی در طول q^{th} بازه ی زمانی بدست آمده استف و یا تعداد نقاطی که \Delta t ثابت است. مزیت نمایش گشتاورهای آماری براساس \gamma این است که Q مجموعه میتوانند با عمل جمع ترکیب شوند و محدودیتی در مقدار Q نیست.


 \gamma_{n,c}= \sum_{q=1}^{Q}\gamma_{n,q} \quad \quad \textrm{for} \quad n=0,1,2,3,4

در اینجا اندیس _c ترکیب بازه های زمانی یا \gammaی ترکیب شده را نشان میدهد. این مقادیر ترکیب شده از \gamma میتواند سپس به صورت معکوس به گشتاورهای خام که تمام بازه ی ترکیبی را نشان میدهند، تبدیل شوند:


 m_{n,c}=\frac{\gamma_{n,c}}{\gamma_{0,c}} \quad \textrm{for} \quad n=1,2,3,4

روابط مشاهده شده بین گشتاورهای خام (m_n) و گشتاورهای مرکزی ( \theta_n = E[(x-\mu)^n]))، سپس برای محاسبه گشتاورهای مرکزیِ ترکیبِ بازه زمانی مورد استفاده قرار می گیرد. در پایان، گشتاورهای آماری بازه ی ترکیبی از گشتاورهای مرکزی محاسبه میشوند:


 \mu_c=m_{1,c}
 \ \ \ \ \ \sigma^2_c=\theta_{2,c}
 \ \ \ \ \ \alpha_{3,c}=\frac{\theta_{3,c}}{\sigma_c^3}
 \ \ \ \ \ \alpha_{4,c}={\frac{\theta_{4,c}}{\sigma_c^4}}-3

کواریانس[ویرایش]

الگوریتمهای کاملاً مشابهی هم میتواند برای محاسبه کواریانس به کار برده شود. الگوریتم naive بدین صورت در می آید:

\operatorname{Cov}(X,Y) = \displaystyle\frac {\sum_{i=1}^n x_i y_i - (\sum_{i=1}^n x_i)(\sum_{i=1}^n y_i)/n}{n}. \!

برای الگوریتم فوق میتوان شبه کد زیر را استفاده کرد:

def naive_covariance(data1, data2):
    n = len(data1)
    sum12 = 0
    sum1 = sum(data1)
    sum2 = sum(data2)
 
    for i in range(n):
        sum12 += data1[i]*data2[i]
 
    covariance = (sum12 - sum1*sum2 / n) / n
    return covariance

یک الگوریتم دو مرحله ای مطمئن تر از لحاظ عددی، ابتدا میانه های میانگین را محاسبه میکند و سپس کواریانس را به دست می آورد:

\bar x = \displaystyle \sum_{i=1}^n x_i/n
\bar y = \displaystyle \sum_{i=1}^n y_i/n
\operatorname{Cov}(X,Y) = \displaystyle\frac {\sum_{i=1}^n (x_i - \bar x)(y_i - \bar y)}{n}. \!

الگوریتم دومرحله ای میتواند به صورت زیر نوشته شود:

def two_pass_covariance(data1, data2):
    n = len(data1)
 
    mean1 = sum(data1) / n
    mean2 = sum(data2) / n	
 
    covariance = 0
 
    for i in range(n):
        a = data1[i] - mean1		
        b = data2[i] - mean2
        covariance += a*b / n
 
    return covariance

یک نسخه ی جبرانی با دقت بالاتر، الگوریتم naive را برروی باقی‌مانده ها انجام میدهد. مجموع نهایی \textstyle\sum y_i و \textstyle\sum x_i باید صفر شود، اما مرحله دوم هرگونه خطای کوچکی را جبران کند. یک الگوریتم تک مرحله ای مطمئنی، مشابه با مورد بالا، وجود دارد، که گشتاورهای مشترک \textstyle C_n = \sum_{i=1}^n (x_i - \bar x_n)(y_i - \bar y_n) را محاسبه میکند:

\bar x_n = \bar x_{n-1} + \frac{x_n - \bar x_{n-1}}{n} \!
\bar y_n = \bar y_{n-1} + \frac{y_n - \bar y_{n-1}}{n} \!
C_n = C_{n-1} + (x_n - \bar x_n)(y_n - \bar y_{n-1}) = C_{n-1} + (y_n - \bar y_n)(x_n - \bar x_{n-1})

عدم تقارن آشکار در تساوی آخر بخاطر این است که \textstyle (x_n - \bar x_n) = \frac{n-1}{n}(x_n - \bar x_{n-1}) ، پس هر دو عبارت به روزرسانی برابرند با: \textstyle \frac{n-1}{n}(x_n - \bar x_{n-1})(y_n - \bar y_{n-1}) . با محاسبه ی میانگین ها میتوان به دقت بالاتری دست یافت، سپس بااستفاده از الگوریتم تک مرحله ای مطمئن برروی بقیه استفاده کرد. بنابراین ما میتوانیم کواریانس را به این صورت به دست آوریم:

\begin{align}
\operatorname{Cov}_N(X,Y) = \frac{C_N}{N} &= \frac{\operatorname{Cov}_{N-1}(X,Y)\cdot(N-1) + (x_n - \bar x_n)(y_n - \bar y_{n-1})}{N}\\
   &= \frac{\operatorname{Cov}_{N-1}(X,Y)\cdot(N-1) + (y_n - \bar y_n)(x_n - \bar x_{n-1})}{N}\\
   &= \frac{\operatorname{Cov}_{N-1}(X,Y)\cdot(N-1) + \frac{N-1}{N}(x_n - \bar x_{n-1})(y_n - \bar y_{n-1})}{N}.
\end{align}

به همین ترتیب، فرمولی برای ترکیب کواریانس دو مجموعه برای موازی سازی محاسبات وجود دارد:

C_X = C_A + C_B + (\bar x_A - \bar x_B)(\bar y_A - \bar y_B)\cdot\frac{n_A n_B}{n_X}.

منابع[ویرایش]

  1. ۱٫۰ ۱٫۱ Bo Einarsson (1 August 2005). Accuracy and Reliability in Scientific Computing. SIAM. p. 47. ISBN 978-0-89871-584-2. Retrieved 17 February 2013. 
  2. T.F.Chan, G.H. Golub and R.J. LeVeque (1983). "Algorithms for computing the sample variance: Analysis and recommendations", The American Statistician, 37. pp. 242–247. 
  3. Higham, Nicholas (2002). Accuracy and Stability of Numerical Algorithms (2 ed) (Problem 1.10). SIAM. 
  4. ۴٫۰ ۴٫۱ Chan, Tony F.; Golub, Gene H.; LeVeque, Randall J. (1979), "Updating Formulae and a Pairwise Algorithm for Computing Sample Variances.", Technical Report STAN-CS-79-773, Department of Computer Science, Stanford University .
  5. دانلد کنوت (1998). هنر برنامه‌نویسی رایانه, volume 2: Seminumerical Algorithms, 3rd edn., p. 232. Boston: Addison-Wesley.
  6. B. P. Welford (1962)."Note on a method for calculating corrected sums of squares and products". Technometrics 4(3):419–420.
  7. Chan, Tony F.; Golub, Gene H.; LeVeque, Randall J. (1983). Algorithms for Computing the Sample Variance: Analysis and Recommendations. The American Statistician 37, 242-247. http://www.jstor.org/stable/2683386
  8. Ling, Robert F. (1974). Comparison of Several Algorithms for Computing Sample Means and Variances. Journal of the American Statistical Association, Vol. 69, No. 348, 859-866. doi:10.2307/2286154
  9. http://www.johndcook.com/standard_deviation.html
  10. D. H. D. West (1979). Communications of the ACM, 22, 9, 532-535: Updating Mean and Variance Estimates: An Improved Method
  11. Terriberry, Timothy B. (2007), Computing Higher-Order Moments Online 
  12. Pébay, Philippe (2008), "Formulas for Robust, One-Pass Parallel Computation of Covariances and Arbitrary-Order Statistical Moments", Technical Report SAND2008-6212, Sandia National Laboratories 
  13. ۱۳٫۰ ۱۳٫۱ Choi, Muenkeun; Sweetman, Bert (2010), Efficient Calculation of Statistical Moments for Structural Health Monitoring 

پیوند به بیرون[ویرایش]