Στατιστικές συναρτήσεις Excel στο Delphi: NORM, CHISQ, BETA

Πληκτρολογήστε =NORM.DIST(115,100,15,TRUE) σε ένα κελί και το Excel επιστρέφει 0,8413447 χωρίς περαιτέρω διαδικασίες. Η κλήση διαβάζεται σαν αναζήτηση. Δεν είναι. Πίσω από αυτόν τον αριθμό κρύβεται η αθροιστική κανονική κατανομή (cumulative normal distribution), ένα ολοκλήρωμα χωρίς κλειστή μορφή, και πίσω από τις CHISQ.INV.RT και BETA.DIST βρίσκονται ειδικές συναρτήσεις τις οποίες μια προσεγμένη βιβλιοθήκη πρέπει να αξιολογήσει, όχι να προσεγγίσει με το χέρι. Ένα υπολογιστικό στοιχείο (spreadsheet component) που ισχυρίζεται συμβατότητα με το Excel πρέπει να αναπαράγει αυτές τις τιμές μέχρι το τελευταίο ψηφίο που εμφανίζει το Excel, πράγμα που σημαίνει αναπαραγωγή των αριθμητικών μεθόδων, όχι μόνο των ονομάτων των συναρτήσεων.

Το HotXLS υλοποιεί περισσότερες από πενήντα από αυτές τις στατιστικές συναρτήσεις, και η εργασία που τις καθιστά σωστές είναι σχεδόν εντελώς αόρατη από τη γραμμή τύπων. Αυτή είναι μια περιήγηση στον τρόπο με τον οποίο τις υπολογίζει η μηχανή: ο κοινός πυρήνας ειδικών συναρτήσεων, οι αποφάσεις διακλάδωσης που διατηρούν την αριθμητική σταθερή, και ένα σφάλμα αντίστροφης κανονικής κατανομής που κρυβόταν στην ουρά (tail) για μεγάλο χρονικό διάστημα επειδή η κοινή περίπτωση δεν άγγιζε ποτέ τη χαλασμένη γραμμή.

Μία κλήση υπολογιστικού φύλλου, πενήντα κατανομές πίσω της

Οι συναρτήσεις καλύπτουν τις οικογένειες που αναζητά ένα βιβλίο εργασίας στατιστικής. Υπάρχει η κανονική οικογένεια, η 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 στο βιβλίο εργασίας μεταγλωττίζει και αξιολογεί έναν ad-hoc τύπο έναντι του ζωντανού φύλλου και επιστρέφει ένα Variant. Μια λεπτομέρεια μπερδεύει τους χρήστες στην πρώτη προσπάθεια: ο parser τύπων πίσω από τη Calculate δέχεται το ελληνικό ερωτηματικό (semicolon) ως διαχωριστικό ορισμάτων, επομένως είναι =SUM(A1;B1), όχι =SUM(A1,B1). Οι αποθηκευμένοι τύποι κελιών διατηρούν το τυπικό κόμμα του Excel. Ο ίδιος αξιολογητής αποστέλλει κάθε στατιστική συνάρτηση παρακάτω, οπότε μόλις μία από αυτές λειτουργήσει στη Calculate, οι υπόλοιπες ακολουθούν την ίδια διαδρομή.

Οι δύο συναρτήσεις στις οποίες βασίζονται όλα τα άλλα

Οι περισσότερες αθροιστικές κατανομές σε αυτό το σύνολο δεν υπολογίζονται με άθροιση ή ολοκλήρωση των δικών τους ορισμών. Υπολογίζονται από δύο ειδικές συναρτήσεις: τη regularized lower incomplete gamma, γραμμένη ως P(a, x), και τη regularized incomplete beta, γραμμένη ως I_x(a, b). Εσωτερικά αυτοί είναι οι βοηθοί στους οποίους βασίζονται οι dispatchers, και η αλυσίδα είναι σύντομη. Η chi-square CDF είναι η gamma CDF με σχήμα df/2 και κλίμακα 2. Η gamma CDF είναι η P(a, x) απευθείας. Οι αθροιστικές συναρτήσεις t, F και διωνυμική είναι όλες τιμές της regularized incomplete beta στα σωστά ορίσματα. Η CDF της Poisson είναι η upper incomplete gamma Q. Υλοποιήστε τις συναρτήσεις γάμμα και βήτα σωστά και δώδεκα κατανομές θα κληρονομήσουν την ακρίβειά τους δωρεάν.

Η λέξη "regularized" είναι όλη η ουσία. Η ακατέργαστη incomplete gamma αυξάνεται σαν παραγοντικό και το ακατέργαστο ολοκλήρωμα βήτα μπορεί να παρουσιάσει υποχείλιση (underflow) ή υπερχείλιση (overflow) πολύ πριν το κάνει η απάντηση. Οι regularized μορφές διαιρούνται με την πλήρη γάμμα ή βήτα, έτσι ώστε να ζουν εξ ολοκλήρου στο διάστημα από μηδέν έως ένα, που είναι ακριβώς το εύρος που καταλαμβάνει μια πιθανότητα. Αυτή η κανονικοποίηση είναι που επιτρέπει στην ίδια ρουτίνα να εξυπηρετεί ένα χι-τετράγωνο με δύο βαθμούς ελευθερίας και ένα με διακόσιους χωρίς οι ενδιάμεσοι όροι να ξεπερνούν το όριο ενός double. Εξηγεί επίσης γιατί δεν υπολογίζετε μια CDF προσθέτοντας μια μακρά ουρά όρων πυκνότητας: κάθε όρος φέρει το δικό του σφάλμα στρογγυλοποίησης, τα σφάλματα συσσωρεύονται καθώς εκτελείται η σειρά, και η regularized ειδική συνάρτηση παρακάμπτει εντελώς το άθροισμα αξιολογώντας αντ' αυτού μια γρήγορα συγκλίνουσα σειρά ή ένα συνεχές κλάσμα.

Σειρά κάτω από τη διαγώνιο, συνεχές κλάσμα πάνω από αυτήν

Η incomplete gamma ρουτίνα λαμβάνει μια απόφαση πριν υπολογίσει οτιδήποτε: συγκρίνει το x με το a + 1. Αυτό το όριο δεν είναι αυθαίρετο. Η ανάπτυξη σε δυναμοσειρά της P(a, x) συγκλίνει γρήγορα όταν το x είναι μικρό σε σχέση με το a, και αργά, τελικά άχρηστα, όταν το x είναι μεγάλο. Το συνεχές κλάσμα έχει τον αντίθετο χαρακτήρα. Έτσι, η μηχανή χρησιμοποιεί τη δυναμοσειρά για x κάτω από a + 1 και ένα συνεχές κλάσμα Lentz για x ίσο ή πάνω από a + 1, και από κάθε κλάδο ζητείται να κάνει μόνο τη δουλειά στην οποία είναι καλός.

Το συνεχές κλάσμα χρειάζεται μια προστασία. Η μέθοδος του Lentz λειτουργεί μεταφέροντας έναν τρέχοντα αριθμητή και παρονομαστή και αντιστρέφοντας τον παρονομαστή σε κάθε βήμα, και αν κάποιο από τα δύο πλησιάσει το μηδέν, η αντιστροφή ανατινάζεται. Η διόρθωση είναι ένα μικροσκοπικό κατώτατο όριο: όποτε ένας ενδιάμεσος όρος πέφτει κάτω από περίπου 1e-30 σε μέγεθος, περιορίζεται (clamped) στο 1e-30, γεγονός που διατηρεί την επανάληψη πεπερασμένη χωρίς να διαταράσσει τη συγκλίνουσα τιμή. Ο ίδιος περιορισμός εμφανίζεται στο συνεχές κλάσμα της incomplete beta για τον ίδιο λόγο. Είναι μια μικρή σταθερά που κάνει σημαντική δουλειά, η διαφορά μεταξύ μιας σταθερής αξιολόγησης και μιας διαίρεσης με κάτι που δεν διακρίνεται από το μηδέν.

Η άνω ουρά, Q(a, x), είναι απλώς 1 μείον P(a, x), και έτσι υπολογίζεται ο αθροιστικός κλάδος της Poisson: η πιθανότητα το πολύ k συμβάντων με μέσο όρο λ είναι Q(k + 1, λ). Η δρομολόγησή της μέσω της άνω incomplete gamma αντί για το άθροισμα k + 1 όρων Poisson είναι, και πάλι, μια επιλογή για την αξιολόγηση μιας συγκλίνουσας έκφρασης αντί για τη συσσώρευση πολλών μικρών.

Διακριτές μάζες χωρίς υπερχείλιση παραγοντικού

Οι διακριτές κατανομές εγείρουν έναν διαφορετικό κίνδυνο. Μια διωνυμική συνάρτηση μάζας πιθανότητας περιλαμβάνει έναν διωνυμικό συντελεστή, και ο συντελεστής για το πενήντα δύο ανά είκοσι έξι είναι ένας τεράστιος ακέραιος. Δημιουργήστε τον απευθείας και ο αριθμητής θα προκαλέσει υπερχείλιση σε double πριν από τη διαίρεση που θα τον επανέφερε σε μια λογική πιθανότητα. Η μηχανή δεν τον δημιουργεί ποτέ. Υπολογίζει τα παραγοντικά σε λογαριθμικό χώρο (log space) μέσω της συνάρτησης log-gamma, προσθέτει και αφαιρεί τους λογαρίθμους, ενσωματώνει το λογάριθμο των πιθανοτήτων επιτυχίας και αποτυχίας, και εκτελεί την ύψωση σε εκθέτη μία φορά στο τέλος.

// 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));

Η ίδια η συνάρτηση log-gamma είναι μια προσέγγιση Lanczos, ακριβής σε ολόκληρο τον θετικό άξονα και οικονομική στην αξιολόγηση. Επειδή κάθε μεγάλη ποσότητα διατηρείται ως λογάριθμός της μέχρι το τελικό Exp, ο μεγαλύτερος αριθμός που υλοποιεί ποτέ η ρουτίνα είναι η ίδια η πιθανότητα, η οποία είναι το πολύ ένα. Η συνάρτηση μάζας της Poisson ακολουθεί την ίδια συνταγή, με τον ενιαίο όρο log-gamma να αντικαθιστά το παραγοντικό στον παρονομαστή. Οι κλειστές μορφές αντιμετωπίζονται ως ειδικές περιπτώσεις στα άκρα, όπου το 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 σε ολόκληρο το εύρος, χωρισμένο σε μια κεντρική περιοχή και δύο ουρές (tails). Πρόκειται για μια αξιολόγηση κλειστής μορφής χωρίς επαναλήψεις.

Οι άλλες αντίστροφες δεν έχουν τέτοιο τύπο, επομένως η μηχανή τις αντιστρέφει αριθμητικά. Οριοθετεί (brackets) την απάντηση με ένα χαμηλό και υψηλό όριο επιλεγμένο από την υποστήριξη της κατανομής, και στη συνέχεια διχοτομεί (bisects): αξιολογεί την CDF προς τα εμπρός στο μέσο, μετακινεί όποιο όριο διατηρεί την πιθανότητα-στόχο εσώκλειστη, και επαναλαμβάνει μέχρι το διάστημα να γίνει στενό. Για τις αντίστροφες γάμμα και χι-τετράγωνο, η οριοθέτηση ξεκινά από το μηδέν και μια γενναιόδωρη άνω εκτίμηση βασισμένη στο σχήμα και την κλίμακα, διπλασιάζοντας το άνω όριο εάν η πιθανότητα δεν είναι ακόμη εσώκλειστη. Η αντίστροφη t οριοθετεί συμμετρικά όρια που διευρύνονται προς τα έξω. Η αντίστροφη F διχοτομεί σε ένα μη αρνητικό διάστημα. Το κόστος είναι μερικές δεκάδες αξιολογήσεις CDF ανά κλήση, κάτι που είναι αόρατο στην ταχύτητα του υπολογιστικού φύλλου, και το όφελος είναι ότι κάθε αντίστροφη είναι ακριβώς τόσο ακριβής όσο η συνάρτηση προς τα εμπρός που αντιστρέφει. Γι' αυτό ένα round trip όπως το CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) επιστρέφει 0.7 με εξαιρετική ακρίβεια.

Ο δεκαδικός λογάριθμος που κρυβόταν στην ουρά

Εδώ είναι το σφάλμα που αξίζει να αναφερθεί, επειδή είναι το είδος που επιβιώνει για μεγάλο χρονικό διάστημα. Η ρουτίνα αντίστροφης κανονικής κατανομής Acklam έχει τρεις κλάδους. Ο ευρύς κεντρικός κλάδος, που χρησιμοποιείται όποτε η πιθανότητα βρίσκεται μεταξύ περίπου 0,025 και 0,975, εκτελεί την εισαγωγή μέσω μιας αναλογίας πολυωνύμων χωρίς λογάριθμο πουθενά σε αυτήν. Οι δύο κλάδοι ουράς, για πολύ μικρές ή πολύ μεγάλες πιθανότητες, λαμβάνουν πρώτα έναν λογάριθμο της εισαγωγής, επειδή η ουρά συμπεριφέρεται όπως η τετραγωνική ρίζα του μείον του φυσικού λογαρίθμου του p.

Μια πρώιμη έκδοση του κλάδου ουράς λάμβανε έναν δεκαδικό λογάριθμο (base-10 log) εκεί όπου ανήκε ο φυσικός λογάριθμος. Οι δύο διαφέρουν κατά έναν σταθερό συντελεστή περίπου 2,30, επομένως τα αποτελέσματα της ουράς ήταν λάθος κατά ένα σταθερό, σημαντικό περιθώριο. Κι όμως, η συνάρτηση φαινόταν μια χαρά σε κάθε τυχαίο έλεγχο, επειδή οι τυχαίοι έλεγχοι ζουν στη μέση. NORM.S.INV(0.5) είναι μηδέν, NORM.S.INV(0.975) είναι το τυπικό 1.959964, και τα δύο αυτά εκτελούνται μέσω του κεντρικού πολυωνύμου που δεν καλεί καθόλου λογάριθμο. Το σφάλμα εμφανιζόταν μόνο όταν μια πιθανότητα περνούσε σε μια ουρά, για παράδειγμα η NORM.S.INV(0.001), η οποία πρέπει να επιστρέφει -3,0902323 και αντίθετα επέστρεφε λάθος τιμή λόγω της αναλογίας φυσικού-έναντι-δεκαδικού λογαρίθμου. Οποιαδήποτε συνάρτηση εξαρτάται από την αντίστροφη κανονική κατανομή στην ουρά της, συμπεριλαμβανομένων των βοηθών διαστήματος εμπιστοσύνης, κληρονομούσε την ίδια απόκλιση. Το μάθημα είναι κοινότοπο και ακριβό: μια συνάρτηση με δομή διακλάδωσης χρειάζεται σημεία δοκιμής μέσα σε κάθε κλάδο, επειδή μια σωστή κοινή διαδρομή θα αποκρύψει ευχάριστα μια χαλασμένη σπάνια. Η διόρθωση ήταν μια αλλαγή ενός token από τον δεκαδικό λογάριθμο στον φυσικό λογάριθμο, και οι τιμές της ουράς ευθυγραμμίστηκαν με αυτές του Excel.

Το πρόσημο του x αποφασίζει την ουρά της t-κατανομής

Η αθροιστική συνάρτηση Student t φέρει μια λεπτομέρεια που είναι εύκολο να κατανοηθεί ανάποδα. Η τιμή της προέρχεται από τη regularized incomplete beta αξιολογημένη στο 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 ώστε το round trip να κλείνει.

Τίποτα από αυτά δεν είναι ορατό από το κελί, και αυτό είναι το επιδιωκόμενο αποτέλεσμα. Γράφετε έναν τύπο και διαβάζετε έναν αριθμό που συμφωνεί με το Excel. Εάν επεκτείνετε τη μηχανή με τη δική σας λογική, οι μηχανισμοί καταχώρισης μιας συνάρτησης καλύπτονται στον οδηγό μας για τη μηχανή τύπων και τις προσαρμοσμένες συναρτήσεις, και ο τρόπος με τον οποίο οι τύποι εκτείνονται σε φύλλα και ονομαστικά εύρη καλύπτεται στο άρθρο για τα ορισμένα ονόματα και τους τύπους μεταξύ φύλλων. Όλα αυτά αποστέλλονται μέσα στο HotXLS spreadsheet component για Delphi και C++Builder, μαζί με τα API ανάγνωσης, εγγραφής, δημιουργίας διαγραμμάτων και μορφοποίησης που καλύπτονται αλλού σε αυτό το ιστολόγιο.