Technical Article

Fungsi Statistik Excel di Delphi: NORM, CHISQ, BETA

Ketik =NORM.DIST(115,100,15,TRUE) ke dalam sel dan Excel mengembalikan 0.8413447 tanpa basa-basi. Panggilan tersebut terbaca seperti pencarian biasa. Padahal tidak. Di balik satu angka itu adalah distribusi normal kumulatif, integral tanpa bentuk tertutup, dan di balik CHISQ.INV.RT serta BETA.DIST terdapat fungsi khusus yang harus dievaluasi oleh pustaka yang teliti, bukan didekati dengan tangan kosong. Komponen spreadsheet yang mengklaim kompatibilitas dengan Excel harus mereproduksi nilai-nilai ini hingga digit terakhir yang ditunjukkan Excel, yang berarti mereproduksi metode numerik, bukan hanya nama fungsinya saja.

HotXLS mengimplementasikan lebih dari lima puluh fungsi statistik ini, dan pekerjaan yang membuatnya benar hampir sepenuhnya tidak terlihat dari bilah rumus (formula bar). Ini adalah tur tentang bagaimana mesin menghitungnya: inti fungsi khusus bersama, keputusan percabangan yang menjaga aritmatika tetap stabil, dan satu bug normal-terbalik (inverse-normal) yang bersembunyi di ekor untuk waktu yang lama karena kasus umum tidak pernah menyentuh baris yang rusak.

Satu panggilan lembar kerja, lima puluh distribusi di belakangnya

Fungsi-fungsi tersebut mencakup keluarga-keluarga yang dicari oleh lembar kerja statistik. Ada keluarga normal, NORM.DIST dan NORM.S.DIST dengan kebalikannya (inverses); keluarga gamma dan chi-square, GAMMA.DIST, CHISQ.DIST, CHISQ.DIST.RT, CHISQ.INV.RT; keluarga beta, BETA.DIST dan BETA.INV; distribusi sampel T.DIST, T.DIST.2T, F.DIST, dan F.INV; pasangan diskrit BINOM.DIST dan POISSON.DIST; dan pembantu inferensi seperti CONFIDENCE.T dan CONFIDENCE.NORM. Dari kursi pemanggil, masing-masing adalah rumus tunggal. Anda menetapkan masukan dalam sel, meminta workbook untuk mengevaluasi, dan membaca hasilnya.

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;

Metode Calculate pada workbook mengompilasi dan mengevaluasi rumus ad-hoc terhadap lembar kerja yang aktif dan mengembalikan Variant. Satu detail menjebak orang pada percobaan pertama: parser rumus di balik Calculate menggunakan titik koma (semicolon) sebagai pemisah argumennya, sehingga rumusnya menjadi =SUM(A1;B1), bukan =SUM(A1,B1). Rumus sel yang disimpan mempertahankan koma standar Excel. Evaluator yang sama mengirimkan setiap fungsi statistik di bawah ini, jadi setelah salah satu dari fungsi ini berfungsi di Calculate, sisanya mengikuti jalur yang sama.

Dua fungsi dasar untuk semua hal lainnya

Sebagian besar distribusi kumulatif dalam set ini tidak dihitung dengan menjumlahkan atau mengintegrasikan definisi mereka sendiri. Mereka dihitung dari dua fungsi khusus: gamma tidak lengkap bawah yang diregulasi (regularized lower incomplete gamma), ditulis P(a, x), dan beta tidak lengkap yang diregulasi (regularized incomplete beta), ditulis Ix(a, b). Secara internal ini adalah pembantu yang diandalkan oleh dispatcher, dan rantainya pendek. CDF chi-square adalah CDF gamma dengan bentuk df/2 dan skala 2. CDF gamma adalah P(a, x) secara langsung. Fungsi kumulatif t, F, dan binomial semuanya merupakan nilai dari beta tidak lengkap yang diregulasi pada argumen yang tepat. CDF Poisson adalah gamma tidak lengkap atas Q (upper incomplete gamma Q). Terapkan fungsi gamma dan beta dengan baik, maka belasan distribusi akan mewarisi keakuratannya secara gratis.

Kata "diregulasi" (regularized) adalah poin utamanya. Gamma tidak lengkap yang mentah tumbuh seperti faktorial dan integral beta mentah dapat mengalami underflow atau overflow jauh sebelum jawabannya keluar. Bentuk yang diregulasi dibagi dengan gamma atau beta lengkap, sehingga mereka hidup sepenuhnya dalam interval dari nol hingga satu, yang merupakan rentang tepat yang ditempati oleh probabilitas. Normalisasi itulah yang memungkinkan rutinitas yang sama melayani chi-square dengan dua derajat kebebasan dan satu dengan dua ratus tanpa istilah perantara yang melampaui batas akhir double. Ini juga menjelaskan mengapa Anda tidak menghitung CDF dengan menjumlahkan ekor panjang dari istilah kepadatan (density terms): setiap istilah membawa kesalahan pembulatannya sendiri, kesalahan tersebut berakumulasi saat deret berjalan, dan fungsi khusus yang diregulasi menghindari jumlah tersebut sepenuhnya dengan mengevaluasi deret yang menyatu dengan cepat atau pecahan berkelanjutan (continued fraction) sebagai gantinya.

Deret di bawah diagonal, pecahan berkelanjutan di atasnya

Rutinitas gamma tidak lengkap membuat satu keputusan sebelum menghitung apa pun: ia membandingkan x terhadap a + 1. Batas tersebut tidak sewenang-wenang. Ekspansi deret-daya (power-series expansion) dari P(a, x) menyatu dengan cepat ketika x kecil relatif terhadap a, dan lambat, bahkan akhirnya tidak berguna, ketika x besar. Pecahan berkelanjutan (continued fraction) memiliki karakter sebaliknya. Jadi mesin menggunakan deret daya untuk x di bawah a + 1 dan pecahan berkelanjutan Lentz untuk x pada atau di atas a + 1, dan setiap cabang diminta untuk melakukan hanya pekerjaan yang dikuasainya.

Pecahan berkelanjutan membutuhkan satu pelindung (guard). Metode Lentz bekerja dengan membawa pembilang dan penyebut yang berjalan dan membalikkan penyebut pada setiap langkah, dan jika salah satu mendekati nol, pembalikan tersebut akan meledak. Perbaikannya adalah batas bawah yang sangat kecil: setiap kali istilah perantara jatuh di bawah sekitar 1e-30 dalam magnitudo, ia dijepit (clamped) ke 1e-30, yang menjaga rekurensi tetap terbatas tanpa mengganggu nilai yang dikonvergensikan. Penjepit yang sama muncul di pecahan berkelanjutan beta tidak lengkap karena alasan yang sama. Ini adalah konstanta kecil yang melakukan pekerjaan penahan beban, perbedaan antara evaluasi yang stabil dan pembagian oleh sesuatu yang tidak dapat dibedakan dari nol.

Ekor atas, Q(a, x), hanyalah 1 dikurangi P(a, x), dan begitulah cara cabang kumulatif Poisson dihitung: probabilitas paling banyak k peristiwa dengan rata-rata λ adalah Q(k + 1, λ). Merutekannya melalui gamma tidak lengkap atas alih-alih menjumlahkan k + 1 istilah Poisson adalah, sekali lagi, pilihan untuk mengevaluasi satu ekspresi konvergen alih-alih mengumpulkan banyak ekspresi kecil.

Massa diskrit tanpa luapan (overflow) faktorial

Distribusi diskrit menimbulkan bahaya yang berbeda. Massa probabilitas binomial melibatkan koefisien binomial, dan koefisien untuk lima-puluh-dua-pilih-dua-puluh-enam adalah integer yang sangat besar. Bentuk secara langsung dan pembilang akan meluap (overflow) dari tipe double sebelum pembagian yang akan membawanya kembali ke probabilitas yang masuk akal. Mesin tidak pernah membentuknya. Ia menghitung faktorial di ruang log (log space) melalui fungsi log-gamma, menambah dan mengurangi log, melipat log probabilitas keberhasilan dan kegagalan, dan mengeksponensialkan sekali di bagian paling akhir.

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

Fungsi log-gamma itu sendiri adalah perkiraan Lanczos, akurat di seluruh sumbu positif dan murah untuk dievaluasi. Karena setiap jumlah besar disimpan sebagai logaritmanya hingga Exp akhir, angka terbesar yang pernah diwujudkan oleh rutinitas adalah probabilitas itu sendiri, yang bernilai paling banyak satu. Fungsi massa Poisson mengikuti resep yang sama, dengan istilah log-gamma tunggal yang mewakili faktorial di penyebut. Bentuk tertutup dijadikan kasus khusus pada tepi-tepinya, di mana p adalah tepat nol atau satu, sehingga kode tidak pernah memanggil Ln(0). HotXLS mengembalikan 0.2460938 untuk BINOM.DIST(5,10,0.5,FALSE) and 0.6766764 untuk kumulatif POISSON.DIST(2,2,TRUE), cocok dengan Excel melalui digit yang dicetaknya.

Kebalikan (inverses) dengan mengurung kurva maju

Fungsi distribusi kebalikan (inverse) menanyakan pertanyaan sebaliknya: diberikan probabilitas, temukan nilai yang CDF-nya sama dengannya. Hanya satu kebalikan dalam set ini yang memiliki rumus langsung yang cepat. NORM.S.INV, kebalikan normal standar, menggunakan perkiraan rasional Acklam, sepasang rasio polinomial yang akurat hingga kira-kira presisi double di seluruh rentang, dibagi menjadi wilayah pusat dan dua ekor. Ini adalah evaluasi bentuk tertutup tanpa iterasi.

Kebalikan lainnya tidak memiliki rumus seperti itu, sehingga mesin membaliknya secara numerik. Ia mengurung (brackets) jawabannya dengan batas bawah dan atas yang dipilih dari dukungan distribusi (distribution's support), lalu membagi dua (bisects): mengevaluasi CDF maju pada titik tengah, memindahkan batas mana yang menjaga probabilitas target tetap tertutup, dan mengulangi sampai intervalnya sempit. Untuk kebalikan gamma dan chi-square, kurung dimulai dari nol dan perkiraan atas yang murah hati dibangun dari bentuk dan skala, menggandakan batas atas jika probabilitas belum tertutup. Kebalikan t mengurung batas simetris yang melebar ke luar; kebalikan F membagi dua pada interval non-negatif. Biayanya adalah beberapa lusin evaluasi CDF per panggilan, yang tidak terlihat pada kecepatan spreadsheet, dan manfaatnya adalah bahwa setiap kebalikan sama akuratnya dengan fungsi maju yang dibaliknya. Itulah mengapa perjalanan pulang pergi seperti CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) mengembalikan 0.7 dalam tingkat akurasi yang sangat tinggi.

Logaritma basis-10 yang tersembunyi di ekor

Berikut adalah bug yang layak diceritakan, karena merupakan jenis yang bertahan lama. Rutinitas normal-terbalik Acklam memiliki tiga cabang. Cabang tengah yang lebar, digunakan setiap kali probabilitas berada di antara sekitar 0.025 dan 0.975, menjalankan masukan melalui rasio polinomial tanpa logaritma di dalamnya. Dua cabang ekor, untuk probabilitas yang sangat kecil atau sangat besar, masing-masing mengambil logaritma dari masukan terlebih dahulu, karena perilaku ekor bertindak seperti akar kuadrat dari minus log alami p.

Versi awal dari cabang ekor mengambil logaritma basis-10 di tempat yang seharusnya merupakan logaritma alami. Keduanya berbeda dengan faktor konstan sekitar 2.30, sehingga hasil ekor salah dengan selisih yang konsisten dan cukup besar. Namun, fungsinya terlihat baik-baik saja dalam setiap pemeriksaan kasual, karena pemeriksaan kasual berada di tengah. NORM.S.INV(0.5) adalah nol, NORM.S.INV(0.975) adalah teks buku 1.959964, dan keduanya dijalankan melalui polinomial pusat yang tidak pernah memanggil logaritma sekali pun. Kesalahan baru muncul setelah probabilitas menyeberang ke ekor, katakanlah NORM.S.INV(0.001), yang harus mengembalikan -3.0902323 dan sebaliknya kembali dengan kesalahan akibat rasio log-alami-versus-basis-10. Fungsi apa pun yang bergantung pada normal terbalik di ekornya, termasuk pembantu interval kepercayaan (confidence-interval helpers), mewarisi kemiringan yang sama. Pelajarannya biasa saja dan mahal: fungsi dengan struktur cabang membutuhkan titik uji di dalam setiap cabang, karena jalur umum yang benar akan dengan senang hati menyembunyikan jalur langka yang rusak. Perbaikannya adalah perubahan satu token dari log basis-10 ke log alami, dan nilai ekor langsung terkunci sesuai dengan Excel.

Tanda x menentukan ekor distribusi-t

Fungsi kumulatif t Student membawa kehalusan yang mudah terbalik. Nilainya berasal dari beta tidak lengkap yang diregulasi yang dievaluasi pada df / (df + x²), tetapi nilai beta tersebut adalah probabilitas di ekor di luar besaran x, bukan probabilitas kumulatif hingga x. Bentuk simetris dari distribusi-t berarti konversi bergantung pada sisi mana dari nol x jatuh.

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

Untuk x di atas nol, probabilitas kumulatif adalah satu dikurangi setengah ekor simetris; untuk x di bawah nol, itu adalah setengah dari ekor tersebut; pada nol, itu adalah tepat setengah. Kembalikan nilai beta secara langsung dan Anda melaporkan sisi distribusi yang salah, meleset oleh seluruh badan kurva untuk x non-nol. Varian ekor kanan dan dua ekor dibangun di atas cabang yang sama, itulah sebabnya T.DIST.2T(1,1) kembali sebagai 0.5 dan T.DIST(1,1,TRUE) sebagai 0.75, dan kebalikan T.INV membagi dua terhadap CDF yang dikoreksi ini sehingga perjalanan pulang pergi ditutup.

Semua ini tidak terlihat dari sel, dan itulah hasil yang dimaksudkan. Anda menulis rumus dan membaca angka yang sesuai dengan Excel. Jika Anda memperluas mesin dengan logika Anda sendiri, mekanisme pendaftaran fungsi dibahas dalam panduan kami tentang mesin rumus dan fungsi kustom, dan cara rumus menjangkau seluruh lembar kerja dan rentang bernama dibahas dalam artikel tentang nama yang ditentukan dan rumus lintas lembar kerja. Semua ini dikirimkan di dalam komponen spreadsheet HotXLS untuk Delphi dan C++Builder, bersama dengan API pembacaan, penulisan, pembuatan bagan, dan pemformatan yang dibahas di bagian lain di blog ini.