توابع آماری Excel در Delphi: توابع NORM ،CHISQ و BETA

عبارت =NORM.DIST(115,100,15,TRUE) را در یک سلول تایپ کنید و Excel بدون هیچ تشریفاتی مقدار 0.8413447 را برمی‌گرداند. این فراخوانی مانند یک جستجو به نظر می‌رسد. اما این‌طور نیست. در پشت آن یک عدد، توزیع نرمال تجمعی (cumulative normal distribution) قرار دارد، انتگرالی بدون فرم بسته، و در پشت CHISQ.INV.RT و BETA.DIST توابع ویژه‌ای قرار دارند که یک کتابخانه دقیق باید آن‌ها را ارزیابی کند، نه اینکه با دست تقریب بزند. یک کامپوننت صفحه گسترده که ادعای سازگاری با Excel را دارد، باید این مقادیر را تا آخرین رقمی که Excel نشان می‌دهد بازتولید کند، که این به معنای بازتولید روش‌های عددی است، نه فقط نام توابع.

مجموعه HotXLS بیش از پنجاه مورد از این توابع آماری را پیاده‌سازی می‌کند و کاری که آن‌ها را درست می‌کند تقریباً به طور کامل از نوار فرمول نامرئی است. این گردشی است در نحوه محاسبه آن‌ها توسط موتور: هسته تابع ویژه مشترک، تصمیمات شاخه‌ای که محاسبات را پایدار نگه می‌دارند و یک باگ معکوس-نرمال که برای مدت طولانی در دنباله پنهان بود زیرا حالت معمول هرگز خط خراب را لمس نمی‌کرد.

یک فراخوانی صفحه کاری، پنجاه توزیع در پشت آن

توابع، خانواده‌هایی را که یک کتاب کار آماری به سراغشان می‌رود در بر می‌گیرند. خانواده نرمال وجود دارد، NORM.DIST و NORM.S.DIST با معکوس‌هایشان؛ خانواده گما و کای‌اسکوئر (chi-square)، یعنی GAMMA.DIST ،CHISQ.DIST ،CHISQ.DIST.RT ،CHISQ.INV.RT؛ خانواده بتا، BETA.DIST و BETA.INV؛ توزیع‌های نمونه‌برداری T.DIST ،T.DIST.2T ،F.DIST و F.INV؛ جفت گسسته BINOM.DIST و POISSON.DIST؛ و هلپرهای استنتاج مانند CONFIDENCE.T و CONFIDENCE.NORM. از جایگاه فراخوان‌کننده، هر کدام یک فرمول واحد هستند. شما ورودی‌ها را در سلول‌ها تنظیم می‌کنید، از کتاب کار می‌خواهید ارزیابی کند و نتیجه را می‌خوانید.

var
  wb: IXLSWorkbook;
  sh: IXLSWorksheet;
begin
  wb := TXLSWorkbook.Create;
  sh := wb.Sheets.Add;
  sh.Range['A1', 'A1'].Value := 115;   // observation
  sh.Range['A2', 'A2'].Value := 100;   // mean
  sh.Range['A3', 'A3'].Value := 15;    // standard deviation

  Writeln(wb.Calculate('=NORM.DIST(A1;A2;A3;TRUE())'));   // 0.8413447
  Writeln(wb.Calculate('=CHISQ.INV.RT(0.05;10)'));        // 18.3070381
  Writeln(wb.Calculate('=BETA.DIST(0.5;2;3;TRUE())'));    // 0.6875
end;

متد Calculate در کتاب کار، یک فرمول موقت را در برابر صفحه زنده کامپایل و ارزیابی می‌کند و یک Variant را پس می‌دهد. یک جزئیات در اولین تلاش افراد را به دردسر می‌اندازد: پارسر فرمول پشت Calculate نقطه‌کاما را به عنوان جداکننده آرگومان خود می‌گیرد، بنابراین به صورت =SUM(A1;B1) است، نه =SUM(A1,B1). فرمول‌های سلول ذخیره‌شده کامای استاندارد Excel را حفظ می‌کنند. همان ارزیاب هر تابع آماری زیر را دیسپچ می‌کند، بنابراین وقتی یکی از این‌ها در Calculate کار کند، بقیه نیز همان مسیر را دنبال می‌کنند.

دو تابعی که بقیه چیزها بر اساس آن‌ها ساخته شده‌اند

بیشتر توابع توزیع تجمعی در این مجموعه با جمع یا انتگرال‌گیری از تعاریف خودشان محاسبه نمی‌شوند. آن‌ها از دو تابع ویژه محاسبه می‌گردند: گامای ناقص پایین نرمال‌شده (regularized lower incomplete gamma) که به صورت P(a, x) نوشته می‌شود، و بتای ناقص نرمال‌شده (regularized incomplete beta) که به صورت I_x(a, b) نوشته می‌شود. به طور داخلی، این‌ها هلپرهایی هستند که دیسپچرها به آن‌ها تکیه می‌کنند و زنجیره کوتاه است. تابع توزیع تجمعی (CDF) کای‌اسکوئر همان CDF گما با شکل df/2 و مقیاس ۲ است. CDF گما مستقیماً همان P(a, x) است. توابع تجمعی t، F و دوجمله‌ای همگی مقادیر بتای ناقص نرمال‌شده در آرگومان‌های صحیح هستند. CDF پوآسون همان گامای ناقص بالای Q است. توابع گما و بتا را به خوبی پیاده‌سازی کنید تا ده توزیع دقت خود را به طور رایگان به ارث ببرند.

کلمه «نرمال‌شده» (regularized) کل نکته است. گامای ناقص خام مانند یک فاکتوریل رشد می‌کند و انتگرال بتای خام می‌تواند خیلی قبل از پاسخ سرریز یا تحت‌ریز (underflow) کند. فرم‌های نرمال‌شده بر گما یا بتای کامل تقسیم می‌شوند، بنابراین کاملاً در بازه صفر تا یک زندگی می‌کنند که دقیقاً همان محدوده‌ای است که یک احتمال اشغال می‌کند. آن نرمال‌سازی همان چیزی است که به یک روتین اجازه می‌دهد به کای‌اسکوئر با دو درجه آزادی و یکی با دویست درجه بدون خارج شدن عبارات میانی از محدوده double خدمت کند. همچنین توضیح می‌دهد که چرا شما یک CDF را با جمع کردن دنباله طولانی از عبارات چگالی محاسبه نمی‌کنید: هر عبارت خطای گرد کردن خود را حمل می‌کند، خطاها با اجرای سری انباشته می‌شوند و تابع ویژه نرمال‌شده با ارزیابی یک سری سریعاً همگرا یا کسر مسلسل به جای آن، کاملاً از جمع کردن عبور می‌کند.

سری‌های زیر قطر، کسر مسلسل بالای آن

روتین گامای ناقص قبل از اینکه هر چیزی را محاسبه کند یک تصمیم می‌گیرد: x را با a + 1 مقایسه می‌کند. آن مرز تصادفی نیست. بسط سری توانی P(a, x) زمانی که x نسبت به a کوچک باشد سریع همگرا می‌شود و زمانی که x بزرگ باشد به آرامی و در نهایت بی‌فایده همگرا می‌گردد. کسر مسلسل ویژگی عکس دارد. بنابراین موتور از سری توانی برای x زیر a + 1 و از کسر مسلسل Lentz برای x در a + 1 یا بالاتر استفاده می‌کند و از هر شاخه خواسته می‌شود فقط کاری را انجام دهد که در آن تخصص دارد.

کسر مسلسل به یک گارد نیاز دارد. روش Lentz با حمل صورت و مخرج در حال اجرا و معکوس کردن مخرج در هر مرحله کار می‌کند و اگر هر کدام به صفر نزدیک شوند، معکوس کردن خراب می‌شود. راه حل یک کف بسیار کوچک است: هر زمان که یک عبارت میانی در اندازه زیر تقریباً 1e-30 قرار گیرد، به 1e-30 محدود (clamp) می‌شود که تکرار را بدون برهم زدن مقدار همگرا محدود نگه می‌دارد. همین محدودیت در کسر مسلسل بتای ناقص نیز به همین دلیل ظاهر می‌شود. این یک ثابت کوچک است که کار بزرگی را انجام می‌دهد، تفاوت بین یک ارزیابی پایدار و تقسیم بر چیزی که از صفر غیرقابل تشخیص است.

دنباله بالایی (upper tail)، یعنی Q(a, x)، به سادگی ۱ منهای P(a, x) است و این نحوه محاسبه شاخه تجمعی پوآسون است: احتمال حداکثر k رویداد با میانگین λ برابر است با Q(k + 1, λ). هدایت آن از طریق گامای ناقص بالایی به جای جمع کردن k + 1 عبارت پوآسون، مجدداً انتخابی برای ارزیابی یک عبارت همگرا به جای انباشتن بسیاری از عبارات کوچک است.

جرم‌های گسسته بدون سرریز فاکتوریل

توزیع‌های گسسته خطر متفاوتی را ایجاد می‌کنند. یک جرم احتمال دوجمله‌ای شامل یک ضریب دوجمله‌ای است و ضریب برای انتخاب ۲۶ از ۵۲ یک عدد صحیح بسیار بزرگ است. آن را مستقیماً شکل دهید و صورت کسر قبل از تقسیمی که آن را به یک احتمال منطقی بازگرداند، یک double را سرریز می‌کند. موتور هرگز آن را شکل نمی‌دهد. فاکتوریل‌ها را در فضای لگاریتمی از طریق تابع لگاریتم-گما (log-gamma) محاسبه می‌کند، لگاریتم‌ها را جمع و تفریق می‌کند، لگاریتم احتمال‌های موفقیت و شکست را اعمال می‌کند و در نهایت یک بار به توان می‌رساند.

// Binomial probability mass, evaluated entirely in log space.
result := Exp(LnGammaF(nt + 1) - LnGammaF(kk + 1) - LnGammaF(nt - kk + 1)
  + kk * Ln(pp) + (nt - kk) * Ln(1 - pp));

خود تابع لگاریتم-گما یک تقریب Lanczos است، که در کل محور مثبت دقیق بوده و ارزیابی آن ارزان است. از آنجا که هر مقدار بزرگ تا Exp نهایی به عنوان لگاریتم خود نگه داشته می‌شود، بزرگ‌ترین عددی که روتین تا کنون مادی می‌کند خود احتمال است که حداکثر یک است. تابع جرم پوآسون از همین دستورالعمل پیروی می‌کند، با این تفاوت که عبارت لگاریتم-گما واحد جایگزین فاکتوریل در مخرج می‌شود. فرم‌های بسته در لبه‌ها، جایی که p دقیقاً صفر یا یک است، به صورت حالت خاص در نظر گرفته می‌شوند تا کد هرگز Ln(0) را فراخوانی نکند. HotXLS مقدار 0.2460938 را برای BINOM.DIST(5,10,0.5,FALSE) و 0.6766764 را برای تجمعی POISSON.DIST(2,2,TRUE) برمی‌گرداند که با ارقام چاپی Excel مطابقت دارد.

معکوس‌ها با محدود کردن منحنی مستقیم

یک تابع توزیع معکوس سوال عکس را می‌پذیرد: با داشتن یک احتمال، مقداری را پیدا کنید که CDF آن با آن برابر باشد. تنها یک معکوس در این مجموعه دارای فرمول مستقیم سریع است. NORM.S.INV، یعنی معکوس نرمال استاندارد، از یک تقریب گویای Acklam استفاده می‌کند، جفتی از نسبت‌های چندجمله‌ای که تقریباً تا دقت یک double در کل محدوده دقیق هستند و به یک منطقه مرکزی و دو دنباله تقسیم شده‌اند. این یک ارزیاب فرم بسته بدون تکرار (iteration) است.

معکوس‌های دیگر چنین فرمول‌هایی ندارند، بنابراین موتور آن‌ها را به صورت عددی معکوس می‌کند. پاسخ را با یک باند پایین و بالا که از دامنه توزیع انتخاب شده محدود می‌کند، سپس دو نیم (bisect) می‌کند: CDF مستقیم را در نقطه میانی ارزیابی می‌کند، هر باندی را که احتمال هدف را محصور نگه می‌دارد جابه‌جا می‌کند و تا زمانی که بازه باریک شود تکرار می‌نماید. برای معکوس‌های گما و کای‌اسکوئر، این محدوده از صفر و یک تخمین بالایی سخاوتمندانه ساخته شده از شکل و مقیاس شروع می‌شود و در صورتی که احتمال هنوز محصور نشده باشد، باند بالایی را دو برابر می‌کند. معکوس t باندهای متقارنی را که به سمت بیرون گسترش می‌بندی محدود می‌کند؛ معکوس F در یک بازه غیرمنفی دو نیم می‌شود. هزینه کار چند ده ارزیابی CDF در هر فراخوانی است که در سرعت صفحه گسترده نامرئی است و مزیت آن این است که هر معکوس دقیقاً به اندازه تابع مستقیمی که معکوس می‌کند دقیق است. به همین دلیل است که یک سفر رفت و برگشت مانند CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) مقدار 0.7 را تا یک تار مو برمی‌گرداند.

لگاریتم مبنای 10 که در دنباله پنهان شده بود

در اینجا باگی وجود دارد که ارزش گفتن دارد، زیرا از آن‌هایی است که برای مدت طولانی زنده می‌ماند. روتین معکوس-نرمال Acklam دارای سه شاخه است. شاخه مرکزی گسترده که هر زمان احتمال بین تقریباً 0.025 و 0.975 باشد استفاده می‌شود، ورودی را از طریق یک نسبت چندجمله‌ای بدون هیچ لگاریتمی در هیچ کجا اجرا می‌کند. دو شاخه دنباله، برای احتمال‌های بسیار کوچک یا بسیار بزرگ، هر کدام ابتدا یک لگاریتم از ورودی می‌گیرند، زیرا دنباله مانند ریشه دوم منفی لگاریتم طبیعی p رفتار می‌کند.

نسخه اولیه شاخه دنباله، یک لگاریتم مبنای ۱۰ را در جایی که لگاریتم طبیعی به آن تعلق داشت گرفت. این دو با یک ضریب ثابت تقریباً ۲.۳۰ با هم تفاوت دارند، بنابراین نتایج دنباله با یک حاشیه ثابت و بزرگ اشتباه بودند. با این حال تابع در هر بررسی معمولی خوب به نظر می‌رسید، زیرا بررسی‌های معمولی در وسط زندگی می‌کنند. NORM.S.INV(0.5) صفر است، NORM.S.INV(0.975) مقدار کتاب درسی 1.959964 است و هر دوی آن‌ها از طریق چندجمله‌ای مرکزی اجرا می‌شوند که اصلاً لگاریتمی را فراخوانی نمی‌کند. خطا تنها زمانی ظاهر شد که یک احتمال وارد یک دنباله شد، مثلاً NORM.S.INV(0.001)، که باید -3.0902323 را برمی‌گرداند و در عوض با نسبت لگاریتم طبیعی در مقابل مبنای ۱۰ تفاوت داشت. هر تابعی که در دنباله خود به معکوس نرمال وابسته است، از جمله هلپرهای فاصله اطمینان، همین انحراف را به ارث برد. درس آموخته‌شده ساده و گران است: تابعی با ساختار شاخه‌ای نیاز به نقاط تست در داخل هر شاخه دارد، زیرا یک مسیر مشترک درست با خوشحالی یک مسیر خراب نادر را ماسک می‌کند. راه حل یک تغییر تک‌توکنی از لاگ مبنای ۱۰ به لاگ طبیعی بود و مقادیر دنباله با مقادیر Excel هم‌تراز شدند.

علامت x دنباله توزیع t را تعیین می‌کند

تابع تجمعی t استودنت ظرافتی را حمل می‌کند که به راحتی معکوس فهمیده می‌شود. مقدار آن از بتای ناقص نرمال‌شده ارزیابی‌شده در df / (df + x²) به دست می‌آید، اما آن مقدار بتا احتمال در دنباله فراتر از اندازه x است، نه احتمال تجمعی تا x. شکل متقارن توزیع t به این معنی است که تبدیل به این بستگی دارد که x در کدام سمت صفر قرار می‌گیرد.

// Student t CDF. ib is the regularized incomplete beta at df/(df+x*x),
ib := BetaIF(df / 2, 0.5, df / (df + x * x));
if x > 0 then
  result := 1 - 0.5 * ib        // above the mean: one minus half the tail
else if x < 0 then
  result := 0.5 * ib            // below the mean: half the tail
else
  result := 0.5;                // exactly at the mean

برای x بالای صفر، احتمال تجمعی یک منهای نصف دنباله متقارن است؛ برای x زیر صفر، نصف آن دنباله است؛ در صفر دقیقاً یک‌دوم است. مقدار بتا را مستقیماً برگردانید تا سمت اشتباه توزیع را گزارش کنید، که برای هر x غیرصفر به اندازه کل بدنه منحنی تفاوت دارد. واریانت‌های دنباله راست و دو دنباله بر روی همان شاخه ساخته می‌شوند، به همین دلیل است که T.DIST.2T(1,1) به عنوان 0.5 و T.DIST(1,1,TRUE) به صورت 0.75 بازمی‌گردد و معکوس T.INV در برابر این CDF اصلاح‌شده دو نیم می‌شود تا رفت و برگشت بسته شود.

هیچ‌کدام از این موارد از سلول قابل مشاهده نیست و این نتیجه مورد نظر است. شما فرمولی می‌نویسید و عددی را می‌خوانید که با Excel مطابقت دارد. اگر در حال گسترش موتور با منطق خودتان هستید، مکانیک ثبت یک تابع در راهنمای ما برای موتور فرمول و توابع سفارشی پوشش داده شده است، و نحوه دسترسی فرمول‌ها به صفحات و محدوده نام‌ها در مقاله نام‌های تعریف شده و فرمول‌های بین‌صفحه‌ای پوشش داده شده است. تمام این‌ها در داخل کامپوننت صفحه گسترده HotXLS برای Delphi و C++Builder در کنار APIهای خواندن، نوشتن، نمودار و قالب‌بندی که در بخش‌های دیگر این وبلاگ پوشش داده شده‌اند، ارائه می‌شوند.