פונקציות סטטיסטיות של Excel ב-Delphi: ‏NORM, CHISQ, BETA

הקלד =NORM.DIST(115,100,15,TRUE) לתוך תא ו-Excel יחזיר 0.8413447 ללא גינונים. הקריאה נראית כמו חיפוש פשוט. היא לא. מאחורי המספר הבודד הזה עומדת התפלגות נורמלית מצטברת, אינטגרל ללא צורה סגורה, ומאחורי CHISQ.INV.RT ו-BETA.DIST יושבות פונקציות מיוחדות שספרייה זהירה חייבת להעריך, ולא לקרב באופן ידני. רכיב גיליונות אלקטרוניים הטוען לתאימות ל-Excel חייב לשחזר את הערכים הללו עד לספרה האחרונה ש-Excel מציג, מה שאומר שחזור השיטות המספריות, ולא רק שמות הפונקציות

התוכנה HotXLS מממשת יותר מחמישים פונקציות סטטיסטיות אלו, והעבודה שהופכת אותן לנכונות כמעט ואינה נראית משורת הנוסחאות. זהו סיור באופן שבו המנוע מחשב אותן: ליבת הפונקציה המיוחדת המשותפת, החלטות ההסתעפות השומרות על יציבות החשבון, ובאג אחד בהתפלגות נורמלית הפוכה שהסתתר בזנב במשך זמן רב מכיוון שהמקרה הנפוץ מעולם לא נגע בשורה השבורה

קריאת גיליון עבודה אחת, חמישים התפלגויות מאחוריה

הפונקציות מקיפות את המשפחות שחוברת עבודה סטטיסטית משתמשת בהן. ישנה המשפחה הנורמלית, NORM.DIST ו-NORM.S.DIST עם ההפוכים שלהן; משפחת גמא וכי בריבוע, 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

  // The XLS formula parser uses ';' as the argument separator.
  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), ופונקציית בטא הלא שלמה המנורמלת, הנכתבת כ-I_x(a, b). באופן פנימי אלו הם העוזרים שהמנתבים נשענים עליהם, והשרשרת קצרה. ה-CDF של כי בריבוע הוא ה-CDF של גמא עם צורה df/2 וקנה מידה 2. ה-CDF של גמא הוא P(a, x) ישירות. פונקציות ההצטברות של t, F ובינומיות הן כולן ערכים של בטא הלא שלמה המנורמלת בארגומנטים הנכונים. ה-CDF של פואסון הוא גמא הלא שלמה העליונה Q. ממש את פונקציות הגמא והבטא היטב ומספר רב של התפלגויות יירשו את הדיוק שלהן בחינם

המילה "מנורמלת" (regularized) היא כל העניין. פונקציית גמא הלא שלמה הגולמית גדלה כמו עצרת ואינטגרל בטא הגולמי יכול לעבור תת-גלישה (underflow) או גלישה זמן רב לפני שהתשובה עושה זאת. הצורות המנורמלות מחולקות דרך הגמא או הבטא השלמות, כך שהן חיות לחלוטין במרווח שבין אפס לאחד, שהוא בדיוק הטווח שהסתברות תופסת. הנורמליזציה הזו היא שמאפשרת לאותה שגרה לשרת כי בריבוע עם שתי דרגות חופש ואחת עם מאתיים מבלי שהמונחים האמצעיים יגלשו מקצה ה-double. היא גם מסבירה מדוע אינך מחשב CDF על ידי חיבור זנב ארוך של מונחי צפיפות: כל מונח נושא שגיאת עיגול משלו, השגיאות מצטברות ככל שהסדרה רצה, והפונקציה המיוחדת המנורמלת עוקפת את הסכום לחלוטין על ידי הערכת סדרה מתכנסת מהר או שבר משולב (continued fraction) במקום זאת

סדרה מתחת לאלכסון, שבר משולב מעליו

שגרת גמא הלא שלמה מקבלת החלטה אחת לפני שהיא מחשבת דבר: היא משווה את x מול a + 1. הגבול הזה אינו שרירותי. פיתוח סדרת החזקות של P(a, x) מתכנס במהירות כאשר x קטן יחסית ל-a, ובאיטיות, עד כדי חוסר תועלת בסופו של דבר, כאשר x גדול. לשבר המשולב יש תכונה הפוכה. לכן המנוע משתמש בסדרת החזקות עבור x מתחת ל-a + 1 ובשבר משולב של Lentz עבור x ב-a + 1 או מעליו, וכל ענף מתבקש לעשות רק את העבודה שהוא טוב בה

השבר המשולב זקוק לשומר אחד. השיטה של Lentz עובדת על ידי נשיאת מונה ומכנה רצים והיפוך המכנה בכל שלב, ואם אחד מהם מתקרב לאפס ההיפוך מתפוצץ. התיקון הוא רצפה זעירה: בכל פעם שמונח ביניים יורד מתחת ל-1e-30 בגודלו הוא ננעל ל-1e-30, מה ששומר על הנסיגה סופית מבלי להפריע לערך שהתכנס. אותו מנעול מופיע בשבר המשולב של בטא הלא שלמה מאותה סיבה. זהו קבוע קטן שעושה עבודה כבדה, ההבדל בין הערכה יציבה לחלוקה במשהו שאינו ניתן להבחנה מאפס

הזנב העליון, Q(a, x), הוא פשוט 1 פחות P(a, x), וכך מחושב הענף המצטבר של פואסון: ההסתברות לכל היותר k אירועים עם ממוצע λ היא Q(k + 1, λ). ניתוב זה דרך הגמא העליונה הלא שלמה במקום סיכום של k + 1 מונחי פואסון הוא, שוב, בחירה להעריך ביטוי מתכנס אחד במקום לצבור מונחים קטנים רבים

מסות בדידות ללא גלישת עצרת

ההתפלגויות הבדידות מעלות סכנה שונה. מסת הסתברות בינומית (binomial probability mass) כוללת מקדם בינומי, והמקדם עבור חמישים ושניים-בחר-עשרים ושש הוא מספר שלם עצום. אם תיצור אותו ישירות, המונה יגלוש מעבר ל-double לפני החלוקה שתחזיר אותו להסתברות הגיונית. המנוע לעולם אינו יוצר אותו. הוא מחשב את העצרות במרחב לוגריתמי דרך פונקציית לוג-גמא, מחבר ומחסיר את הלוגריתמים, משלב את הלוג של הסתברויות ההצלחה והכישלון, ומבצע אקספוננציאציה פעם אחת ממש בסוף

// Binomial probability mass, evaluated entirely in log space.
// LnGammaF(n+1) is ln(n!); the three log-factorials form ln(C(n,k)),
// and the whole exponent is built before a single Exp call.
//   ln P(X=k) = ln(n!) - ln(k!) - ln((n-k)!) + k*ln(p) + (n-k)*ln(1-p)
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 לאורך כל הטווח, המחולק לאזור מרכזי ושני זנבות. זו הערכה בצורה סגורה ללא איטרציות

להפוכים האחרים אין נוסחה כזו, ולכן המנוע הופך אותם מספרית. הוא כובל את התשובה בגבול נמוך וגבוה שנבחרו מתחום ההתפלגות, ואז חוצה (bisects): מעריך את ה-CDF הישיר בנקודת האמצע, מזיז את הגבול ששומר על הסתברות היעד כלואה, וחוזר חלילה עד שהמרווח צר. עבור הפוכי גמא וכי בריבוע הכליאה מתחילה באפס והערכה עליונה נדיבה שנבנית מהצורה וקנה המידה, ומכפילה את הגבול העליון אם ההסתברות עדיין אינה כלואה. הפוך t כובל גבולות סימטריים המתרחבים החוצה; הפוך F חוצה במרווח אי-שלילי. העלות היא כמה עשרות הערכות CDF לכל קריאה, מה שאינו מורגש במהירות של גיליון אלקטרוני, והיתרון הוא שכל הפוך מדויק בדיוק כמו הפונקציה הישירה שהוא הופך. זו הסיבה שסבב מלא כגון CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) מחזיר 0.7 בדיוק רב

הלוגריתם בבסיס 10 שהסתתר בזנב

גרסה מוקדמת של ענף הזנב לקחה לוגריתם בבסיס 10 במקום שהלוגריתם הטבעי היה שייך לו. השניים נבדלים בגורם קבוע של כ-2.30, כך שתוצאות הזנב היו שגויות בהפרש עקבי ומשמעותי. ועדיין הפונקציה נראתה תקינה בכל בדיקה אקראית, מכיוון שבדיקות אקראיות חיות באמצע. NORM.S.INV(0.5) הוא אפס, NORM.S.INV(0.975) הוא ה-1.959964 מהספרים, ושניהם רצים דרך הפולינום המרכזי שלעולם אינו קורא ללוגריתם כלל. השגיאה הופיעה רק ברגע שההסתברות חצתה לתוך זנב, נניח NORM.S.INV(0.001), שחייב להחזיר -3.0902323 ובמקום זאת חזר ליד יחס הלוג הטבעי מול בסיס 10 משובש. כל פונקציה התלויה בנורמלי ההפוך בזנב שלו, כולל עוזרי מרווח בר-סמך, ירשה את אותה הסטייה. הלקח הוא יומיומי ויקר: פונקציה עם מבנה ענפים זקוקה לנקודות בדיקה בתוך כל ענף, מכיוון שנתיב נפוץ נכון ימסך בשמחה נתיב נדיר שבור. התיקון היה שינוי של אסימון (token) בודד מהלוג של בסיס 10 ללוג הטבעי, וערכי הזנב התיישרו בדיוק לאלו של 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),
// which measures the symmetric tail. The cumulative value depends on
// the sign of x; returning ib unconverted gives the wrong tail.
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 spreadsheet component עבור Delphi ו-C++Builder, לצד ממשקי לקריאה, כתיבה, שרטוט תרשימים ועיצוב המכוסים במקומות אחרים בבלוג זה