اكتب =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، تتبع البقية نفس المسار.
الدالتان اللتان يبنى عليهما كل شيء آخر
معظم التوزيعات التراكمية في هذه المجموعة لا تُحسب بجمع أو تكامل تعريفاتها الخاصة. بل تُحسب من دالتين خاصتين: دالة جاما غير المكتملة السفلية المنظمة، وتُكتب P(a, x), ودالة بيتا غير المكتملة المنظمة، وتُكتب Ix(a, b). داخلياً، هذه هي المساعدات التي تعتمد عليها الموزعات، والسلسلة قصيرة. دالة التوزيع التراكمي (CDF) لكاي تربيع هي CDF لجاما مع معلمة الشكل df/2 والمقياس 2. وتكون CDF لجاما هي P(a, x) مباشرة. ودوال التراكم لـ t وF وثنائي الحدين هي جميعها قيم لدالة بيتا غير المكتملة المنظمة عند الوسائط الصحيحة. وتكون CDF لبواسون هي دالة جاما غير المكتملة العلوية Q. نفذ دالتي جاما وبيتا جيداً وستحصل عشرات التوزيعات على دقتها مجاناً.
كلمة "المنظمة" (regularized) هي النقطة الأساسية. تنمو دالة جاما غير المكتملة الخام مثل المضروب (factorial) ويمكن أن يقل تدفق تكامل بيتا الخام أو يتجاوز الحد قبل وقت طويل من وصول الإجابة. تُقسم الأشكال المنظمة على جاما أو بيتا الكاملة، لذا تعيش بالكامل في الفاصل من صفر إلى واحد، وهو بالضبط النطاق الذي يشغله الاحتمال. هذا التنظيم هو ما يسمح للروتين نفسه بخدمة كاي تربيع مع درجتين من الحرية وتوزيع آخر بمائتين دون أن تخرج الحدود الوسيطة عن نهاية عدد مزدوج الدقة. وهو يوضح أيضاً لماذا لا تحسب CDF بجمع ذيل طويل من حدود الكثافة: يحمل كل حد خطأ التقريب الخاص به، وتتراكم الأخطاء مع تشغيل السلسلة، وتتجنب الدالة الخاصة المنظمة المجموع بالكامل عن طريق تقييم سلسلة سريعة التقارب أو كسر مستمر بدلاً من ذلك.
السلسلة تحت القطر، والكسر المستمر فوقه
يتخذ روتين جاما غير المكتمل قراراً واحداً قبل أن يحسب أي شيء: يقارن 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 من حدود بواسون هو، مرة أخرى، اختيار لتقييم تعبير متقارب واحد بدلاً من تراكم العديد من التعبيرات الصغيرة.
الكتل المتقطعة بدون تجاوز المضروب
تثير التوزيعات المتقطعة خطراً مختلفاً. تشتمل كتلة الاحتمال ثنائي الحدين على معامل ثنائي الحدين، والمعامل لخمسين اختياراً لستة وعشرين هو عدد صحيح هائل. شكله مباشرة وسيتجاوز البسط قيمة 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 عبر النطاق بأكمله، ومقسمة إلى منطقة مركزية وذيلين. إنه تقييم ذو شكل مغلق بدون تكرار.
المقلوبات الأخرى ليس لها مثل هذه الصيغة، لذا يقوم المحرك بعكسها عددياً. فهو يحصر الإجابة بحد أدنى وأقصى يتم اختيارهما من دعم التوزيع، ثم ينصف: يقيم CDF المباشرة عند نقطة المنتصف، وينقل أي حد يحافظ على الاحتمال المستهدف محصوراً، ويكرر ذلك حتى يضيق الفاصل الإضافي. بالنسبة لمقلوبات جاما وكاي تربيع، يبدأ الحصر عند الصفر وتقدير علوي سخي مبني على الشكل والمقياس، مع مضاعفة الحد العلوي إذا لم يكن الاحتمال محصوراً بعد. وييحصر مقلوب t حدوداً متناظرة تتسع للخارج، وينصف مقلوب F على فاصل غير سالب. التكلفة هي بضع عشرات من تقييمات CDF لكل استدعاء، وهو أمر غير مرئي بسرعة جدول البيانات، والفائدة هي أن كل مقلوب دقيق تماماً مثل الدالة المباشرة التي يعكسها. وهذا هو السبب في أن رحلة ذهاب وإياب مثل CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) تعيد 0.7 بدقة متناهية.
اللوغاريتم ذو الأساس 10 الذي اختبأ في الذيل
إليك الخطأ البرمجي الذي يستحق الحديث عنه، لأنه من النوع الذي يعيش طويلاً. يحتوي روتين Acklam للتوزيع الطبيعي العكسي على ثلاثة فروع. الفرع المركزي الواسع، المستخدم عندما يقع الاحتمال بين 0.025 و0.975 تقريباً، يمرر الإدخال عبر نسبة متعددة الحدود دون وجود لوغاريتم في أي مكان. ويأخذ فرعا الذيلين، للاحتمالات الصغيرة جداً أو الكبيرة جداً، لوغاريتم الإدخال أولاً، لأن الذيل يتصرف مثل الجذر التربيعي لسالب اللوغاريتم الطبيعي لـ p.
أخذت نسخة مبكرة من فرع الذيل لوغاريتم الأساس 10 حيث يجب أن يكون اللوغاريتم الطبيعي. يختلف الاثنان بمعامل ثابت يبلغ حوالي 2.30, لذا كانت نتائج الذيل خاطئة بهامش ثابت وكبير. ومع ذلك، بدت الدالة جيدة في كل فحص عشوائي، لأن الفحوصات العشوائية تعيش في المنتصف. NORM.S.INV(0.5) هي صفر، وNORM.S.INV(0.975) هي القيمة القياسية 1.959964، وكلاهما يمر عبر متعددة الحدود المركزية التي لا تستدعي لوغاريتم على الإطلاق. ظهر الخطأ فقط بمجرد عبور الاحتمال إلى الذيل، على سبيل المثال NORM.S.INV(0.001)، والتي يجب أن تعيد -3.0902323 وجاءت بدلاً من ذلك قريبة من نسبة اللوغاريتم الطبيعي مقابل الأساس 10 الخاطئة. وورثت أي دالة تعتمد على التوزيع الطبيعي العكسي في ذيلها، بما في ذلك مساعدات فترات الثقة، نفس الانحراف. الدرس بسيط ومكلف: الدالة ذات بنية الفروع تحتاج إلى نقاط اختبار داخل كل فرع، لأن المسار المشترك الصحيح سيقنعك بصحته ويخفي فرعاً معيباً نادراً. وكان الإصلاح عبارة عن تغيير رمز واحد من لوغاريتم الأساس 10 إلى اللوغاريتم الطبيعي، واستقرت قيم الذيل لتطابق قيم Excel.
إشارة x تحدد ذيل توزيع t
تحمل الدالة التراكمية لتوزيع Student 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 لـ Delphi وC++Builder، إلى جانب واجهات برمجة تطبيقات القراءة والكتابة والرسم البياني والتنسيق المغطاة في مكان آخر من هذه المدونة.