الگوریتم محاسبه واریانس
الگوریتم Naïve
فرمول محاسبه واریانس جامعه آماری با اندازه N به این صورت است:
فرمول محاسبه تخمین unbiased از واریانس جامعه آماری از یک نمونه محدود از n مشاهده به این صورت میباشد:
بنابراین الگوریتم 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
میتوانند مقدار بسیار نزدیکی داشته باشند، تقریب زدن میتواند منجر شود که دقت نتیجه بسیار کمتر از دقت ذاتی موجود در محاسبات ممیز شناور شود. بنابراین این الگوریتم نباید در عمل مورد استفاده قرار گیرد. این موضوع بخصوص در مواردی که انحراف معیار نسبت به میانگین بسیار کوچک است، تشدید میشود. اگرچه، این الگوریتم میتواند با استفاده از روش میانگین مفروض بهبود یابد.
الگوریتم دومرحلهای
یک روش جایگزین، استفاده از فرمول دیگری برای محاسبه واریانس میباشد، به این صورت که ابتدا میانگین نمونه را محاسبه میکند:
- ,
و سپس مجموع مربعات فاصله از میانگین را بدست میآورد:
- ,
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
الگوریتم برخط
محاسبهٔ واریانس در یک مرحله، اغلب میتواند بسیار مفید باشد، به عنوان مثال، وقتی دادهها در حال جمعآوری هستند و فضای کافی برای ذخیره همه دادهها نداریم، یا موقعی که هزینه دسترسی به حافظه، بر هزینه محاسبات تأثیر زیادی دارد. برای چنین الگوریتم برخطی، رابطه رخداد مجدد بین کمیتهایی مورد استفاده در محاسبات، مورد نیاز میباشد تا محاسبات ما به صورت به روشی سازگار انجام شود.
فرمولهای زیر برای به روزرسانی میانگین و واریانس(تخمینی) دادهها، در زمان افزودن عنصر جدید
همچنین میتوان نشان داد که کمیت مناسبتر برای به روزرسانی مجموع مربعات فاصله از میانگین کنونی (
یک الگوریتم عددی سازگار به صورت زیر داده شدهاست، که البته میانگین را نیز محاسبه می نماید. این الگوریتم، منتسب به Knuth است که او نیز از الگوریتمی از Welford استفاده کرده و به صورت کامل تحلیل نمودهاست. همچنین معمول است که متغیرهای زیر به جای همدیگر استفاده شوند:
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. میگوید که سومین الگوریتم برخط فوق، نمونهای خاص از الگوریتم است که برای هر تقسیمبندی از نمونه
- .
این موضوع میتواند زمانی مفید باشد که برای مثال، چندین واحد پرداش برای قسمتهای مختلف از ورودی در نظر گرفته شده باشد.
شیوه Chan روشی برای تخمین میانگین، زمانی که
مثال
فرض کنید که عملیات ممیز شناور از استاندارد محاسباتی IEEE 754 double-precision استفاده میکند. نمونهٔ {4، 7، 13، 16} را از یک جامعه آماری نامحدود فرض کنید. براساس این نمونه، میانگین تخمینی جامعه برابر با 10 است و واریانس تخمینی بدون پیشقدرِ جامعه برابر با 30 است. هر دو الگوریتم اول و دوم این مقادیر را به درستی محاسبه میکند. سپس نمونه { 10 + 4, 10 + 7, 10 + 13, 10 + 16 } را در نظر بگیرید، که مقدار واریانسی برابر با نمونهٔ قبلی به ما میدهد. الگوریتم دوم این واریانس تخمینی را به درستی محاسبه میکند، اما الگوریتم اول مقدار 29.33333333 را به جای 30 برمیگرداند. اگرچه این کمبود دقت قابل تحمل است و به عنوان نقص کوچکی در الگوریتم اول مشاهده میشود، اما به سادگی دادههایی یافت که نقص اساسی در الگوریتم naive را آشکار میسازد: نمونهٔ { 10 + 4, 10 + 7, 10 + 13, 10 + 16 } را در نظر بگیرید. دوباره واریانس تخمینی برابر با 30 میباشد که توسط الگوریتم دوم به سادگی قابل محاسبه است، اما این دفعه الگوریتم naive مقدار −170.66666666666666 را به عنوان جواب برمیگرداند. این یک مشکل اساسی در الگوریتم اول میباشد و بخاطر تقریب زدن فاجعه آمیزی رخ داده است که در مرحله آخر الگوریتم، هنگام تفریق دو عدد مشابه صورت گرفتهاست.
آمارهای رده بالاتر
Terriberry فرمول Chan را برای محاسبه سومین و چهارمین گشتاور مرکزی تعمیم میدهد. برای مثال برای تخمین چولگی و کشیدگی مورد نیاز است:
اینجا
- چولگی:
- کشیدگی:
برای نسخه افزایشی (مثلاً
با نگهداشتن مقدار
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 در کاربردهای معینی، صرفه جویی کند. روش اول، برای محاسبه گشتاورهای آماری، با جداسازی دادهها، گشتاورها را از رابطه هندسیِ هیستوگرامِ به دست آمده محاسبه میکند، که به صورت بهینهای تبدیل به یک الگوریتم یک مرحلهای برای گشتاورهای بالاتر میشود. یکی از مزایای کار این است که محاسبات گشتاورهای آماری میتواند با دقت اختیاری محاسبه شود، پس میتوان این دقت را با دقت فرمت انبار داده ویا سخت افزار اندازهگیری اصلی، منطبق کرد. هیستوگرام نسبی از یک متغیر تصادفی میتواند یک در این راه قراردادی ساخته شود: محدوده مقادیر بالقوه به بازههایی تقسیم میشود و تعداد رخدادها در هر بازه شمارش و به گونهای کشیده میشود که مساحت هر مستطیل برابر با سهم رخدادهای آن بازه شود:
که
که توان
دومین روش Choi و Sweetman یک متدولوژی تحلیلی برای ترکیب گشتاورهای آماری از بخشهای مجزا از یک بازهٔ زمانی است که گشتاورهای کلی بدست آمده، همان گشتاورهای بازهٔ زمانی کامل است. این متدولوژی میتواند برای محاسبه موازی گشتاورهای آماری با کمک ترکیب گشتاورهای پس آیند(بعدی) ویا ترکیب گشتاورهای آماری محاسبه شده در زمانهای پشت سرهم، استفاده شود.
اگر تعداد
که
در اینجا اندیس
روابط مشاهده شده بین گشتاورهای خام (
کوواریانس
الگوریتمهای کاملاً مشابهی هم میتواند برای محاسبه کوواریانس به کار برده شود. الگوریتم naive بدین صورت در میآید:
برای الگوریتم فوق میتوان شبه کد زیر را استفاده کرد:
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
یک الگوریتم دو مرحلهای مطمئن تر از لحاظ عددی، ابتدا میانههای میانگین را محاسبه میکند و سپس کوواریانس را به دست میآورد:
الگوریتم دومرحلهای میتواند به صورت زیر نوشته شود:
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 را بر روی باقیماندهها انجام میدهد. مجموع نهایی
عدم تقارن آشکار در تساوی آخر بخاطر این است که
به همین ترتیب، فرمولی برای ترکیب کوواریانس دو مجموعه برای موازیسازی محاسبات وجود دارد:
- .
منابع
- ↑ Bo Einarsson (1 August 2005). Accuracy and Reliability in Scientific Computing. SIAM. p. 47. ISBN 978-0-89871-584-2. Retrieved 17 February 2013.
- ↑ T.F.Chan, G.H. Golub and R.J. LeVeque (1983). ""Algorithms for computing the sample variance: Analysis and recommendations", The American Statistician, 37" (PDF): 242–247.
- ↑ Higham, Nicholas (2002). Accuracy and Stability of Numerical Algorithms (2 ed) (Problem 1.10). SIAM.
- ↑ Chan, Tony F.; Golub, Gene H.; LeVeque, Randall J. (1979), "Updating Formulae and a Pairwise Algorithm for Computing Sample Variances." (PDF), Technical Report STAN-CS-79-773, Department of Computer Science, Stanford University.
- ↑ دانلد کنوت (1998). هنر برنامهنویسی رایانه, volume 2: Seminumerical Algorithms, 3rd edn., p. 232. Boston: Addison-Wesley.
- ↑ B. P. Welford (1962)."Note on a method for calculating corrected sums of squares and products". Technometrics 4(3):419–420.
- ↑ 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
- ↑ 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
- ↑ http://www.johndcook.com/standard_deviation.html
- ↑ D. H. D. West (1979). Communications of the ACM, 22, 9, 532-535: Updating Mean and Variance Estimates: An Improved Method
- ↑ Terriberry, Timothy B. (2007), Computing Higher-Order Moments Online, archived from the original on 23 April 2014, retrieved 19 December 2013
- ↑ Pébay, Philippe (2008), "Formulas for Robust, One-Pass Parallel Computation of Covariances and Arbitrary-Order Statistical Moments" (PDF), Technical Report SAND2008-6212, Sandia National Laboratories, archived from the original (PDF) on 9 اكتبر 2022, retrieved 19 December 2013
- ↑ Choi, Muenkeun; Sweetman, Bert (2010), Efficient Calculation of Statistical Moments for Structural Health Monitoring (PDF), archived from the original (PDF) on 3 March 2016, retrieved 19 December 2013