Vìøit statistickému software? Josef Tvrdík
Ostravská universita 1
Abstrakt. Pøíspìvek se zabývá selháními softwarových statistických procedur, zji¹tìnými v prùbìhu jejich dlouhodobého u¾ívání. Jsou diskutovány nìkteré chyby v jednoduchých popisných metodách nalezené v Excelu, nesprávné odhady parametrù nelineárních regresních modelù vyskytující se a¾ pøíli¹ èasto v bì¾nì prodávaných statistických paketech (NCSS, SYSTAT, SPSS, S-PLUS) a také numerické nesrovnalosti i triviální implementaèní chyby ve statistikách pro test shody funkcí pøe¾ití v NCSS. Klíèová slova: statistický software, Excel 97, algoritmy, chyby, nelineární regrese.
Úvod
Text je roz¹íøenou versí sdìlení pøedneseného na Statistických dnech Èeské statistické spoleènosti v Ostravì v èervnu 2000. Pod názvem Opravdu jen drobné vady na kráse statistického software?\ byl pøed"nesen také na letní ¹kole ROBUST 2000 v Neètinách v záøí tého¾ roku. Obì tato upozornìní na nìkteré zjevné i ménì zjevné numerické nesrovnalosti èi chyby zji¹tìné pøi u¾ívání známých statistických paketù a tabulkového procesoru Excel 97 vyvolala dosti ¾ivý ohlas u úèastníkù zmínìných akcí. Po dohodì s editory sborníku ROBUST je tento text zveøejnìn v Bulletinu Èeské statistické spoleènosti, kde snad má pøíle¾itost se dostat k vìt¹ímu okruhu zainteresovaných ètenáøù. Excel je pøidán ke kritizovaným statistickým programùm, nebo» je pro statistické výpoèty velmi èasto u¾íván zejména lidmi potøebujícími statistiku pouze obèas a u nich je rozpoznání chybného výsledku je¹tì ménì pravdìpodobné ne¾ u zku¹eného statistika. Chyby v Excelu mohou tedy pùsobit ¹kodu velice èasto. Navíc ¹kody pùsobí také èeská lokalizace Excelu, viz Tvrdík (1998). Po tìchto zku¹enostech s nespolehlivostí statistického software se stává naléhavou otázka, zda úsilí, které statistici vìnují hledání rigorosních øe¹ení rùzných statistických problémù (obèas i dosti vyumìlkova1 Tato práce byla podporována z grantu 402/00/1165 GA ÈR a z projektu institucionálního výzkumu CEZ: J09/98:179000002.
ných) není z hlediska aplikací statistiky zbyteèné a zda by podobné úsilí nemìlo být orientováno na výbìr vhodnìj¹ích, numericky spolehlivìj¹ích algoritmù a dùkladnìj¹ímu testování jejich implementace. Naprostá vìt¹ina aplikací statistiky je opøena o výpoèty provedené s vyu¾itím statistického software a pokud jsou jejich výsledky numericky chybné, jsou so stikované statistické metody na nic.
Testy numerické spolehlivosti
Testováním spolehlivosti statistických programù se zabývají v¹ichni výrobci software, ale patrnì nìkteré výsledky si nechávají pro sebe. Objektivní pohled na spolehlivost statistických procedur je asi dost problematický. O jednu z mo¾ných cest se pokou¹í americký Státní institut standardù a technologie (National Institute of Standards and Technology, NIST, viz citaci na jeho web-stránku). Tam je shromá¾dìna sada testovacích úloh z nìkolika oblastí statistiky, u kterých jsou známé jejich výsledky na jistý (tzv. certi kovaný) poèet platných cifer. Pøehled je uveden v tab. 1. Tabulka 1: Standardní referenèní úlohy NIST { pøehled Druh poèet certi kovaný poèet úloh úloh platných míst jednorozmìrné statistiky 9 15 lineární regrese 11 15 analýza rozptylu 11 15 nelineární regrese 27 11 Numerickou správností výsledkù nìkterých statistických procedur v Excelu 97 se nedávno zabývali McCullough a Wilson (1999). Zji¹»ovali, jak se shodují výsledky získané Excelem s certi kovanými výsledky NIST. Ponìkud paradoxní je, ¾e v jejich èlánku je chyba v de nici velièiny, kterou sledovali. Podle slovního popisu má velièina vyjadøovat míru shody výsledné hodnoty x spoèítané Excelem s certi kovanou hodnotou c a znamená vlastnì poèet platných èíslic shodných s certi kovaným výsledkem. V recenzovaném èlánku v pomìrnì renomovaném èasopisu jim pro¹la následující de nièní rovnice = log10 (jx
cj) =jcj
(1)
Správná de nice má mít zøejmì tvar 8 > > <
kdy¾ jxjcjcj 1 jx cj 15 = (2) kdy¾ jcj < 1 10 > > : log10 jx cj jinak jcj V úlohách, kde se poèítá více ne¾ jeden parametr, je výsledná míra shody pro úlohu chápána jako 0 15
= min(1 ; 2 ; : : : ; k );
(3)
kde k je poèet vypoèítávaných parametrù.
Jednorozmìrné statistiky
Pøes vý¹e uvedené výhrady v¹ak McCullough a Wilson (1999) dùvìryhodnì zjistili, ¾e Excel selhává i v jednorozmìrných statistikách, kdy dokonce mezi úlohami v sadì NIST byly nalezeny takové. pro které = 0. Bylo to zpùsobeno vìt¹inou u¾itím nevhodného algoritmu pro výpoèet výbìrového rozptylu. V mnoha uèebnicích základních statistických metod se tradiènì uvádí, ¾e výbìrový rozptyl je sx 2 =
n
2 n 1 4X xi 2
n 1 X
2 1 i=1 (xi x) = n 1
i=1
n 1 X x n i=1 i
!2 3 5;
(4)
pøièem¾ výraz za druhým rovnítkem se doporuèuje jako výpoèetnì vhodnìj¹í. Na o¹idnost toho doporuèení upozoròuje Ekblom (1994). Jak je rovnì¾ v uèebnicích elementárních statistických metod uvádìno, rozptyl je invariantní vùèi posunu, tj. pro yi = a + xi , a je konstantní, je pak sy 2 = sx 2 . Tento vztah má v¹ak pøi numerických výpoètech omezenou platnost, nebo» musíme uva¾ovat chyby ze zaokrouhlování. Pokud je prùmìr x velký a rozptyl malý, pak druhá rovnost v rov. 4 platí jen pøibli¾nì, za jistých okolností mù¾e být poèítaèová hodnota výrazu S=
n X i=1
xi 2
n 1 X xi n i=1
!2
(5)
dokonce záporná. Podle experimentálních výsledkù nìkolika testovacích pøíkladù lze usoudit, ¾e v Excelu je tento problém vyøe¹en vskutku velice úspornì. Byl pøímo ukázkovì uplatnìn i odjinud známý racionální pøístup Microsoftu: Pokud nastane situace, ¾e hodnota výrazu S z rov. 5 je men¹í ne¾ 0, pak výsledný výbìrový rozptyl v Excelu je roven nule. Ve standardním statistickém software podobná hrubá implementaèní chyba zji¹tìna nebyla, napø. NCSS dává správné výsledky pro v¹echny úlohy z testovací sady NIST.
Nelineární regrese
Na sadì pøíkladù NIST McCullough a Wilson (1999) ovìøovali rovnì¾ numerickou správnost odhadù parametrù nelineárních regresních modelù získaných Excelem. Za selhání programu se pova¾uje, kdy¾ program skonèí výpoèet v lokálním minimu souètu ètvercù residuálních odchylek, tzn. = 0. Na stejných úlohách testovala Valchaøová (2000) statistický paket NCSS 6.0, pro ka¾dou úlohu 10 opakování s náhodnì volenými poèáteèními hodnotami odhadù z jejich pøijatelného oboru. Za selhání se pova¾uje, kdy¾ ve více ne¾ tøetinì opakování výpoètu pro danou úlohu skonèí odhad parametrù v lokálním minimu. Výsledky uvedeny v tab. 2. Pohled na tabulku dùvìru ve spolehlivost statistických programù nepovzbudí. Tabulka 2: Nelineární regrese { poèty selhání obtí¾nost poèet Excel NCSS úloh úloh McCullough, Wilson Valchaøová nízká 8 4 2 støední 11 10 6 vysoká 8 7 3 celkem 27 21 11 Na 8 z 11 v NCSS selhávajících úloh zkusila Valchaøová u¾ít pro odhad parametrù stochastické algoritmy, popis algoritmù viz Tvrdík a Køivý (1999). Aè byly výpoèty provádìny v èasové tísni pøed termínem odevzdání diplomové práce bez hlub¹ího rozmyslu (ale to je dosti èastý pøístup u¾ivatelù statistického software), v polovinì z tìchto úloh stochastické algoritmy neselhaly ani jednou, tak¾e z testu vy¹ly pøece jen o trochu lépe ne¾ softwarová klasika.
Zajímavé srovnání úspì¹nosti statistického software (Excel, S-Plus, SPSS, SAS, STATA, Mathematica) na sadì úloh NIST publikoval nedávno McCoullogh (2000). Jediným plnì úpì¹ným softwarovým produktem na v¹ech 59 úlohách NIST byla Mathematica. Tabulka 3: Procento selhání standardních statistických paketù Model NCSS SYSTAT SYSTAT S-PLUS SPSS Prùmìr G-N Simplex 1 0 3 0 6 97 21.2 90 35 37 73 77 62.2 2 3 89 69 67 100 100 85.0 4 4 0 76 0 0 16.0 0 16 57 3 35.2 5 100 6 45 8 0 100 0 30.6 7 100 100 100 81 69 90.0 78 11 30 18 0 27.4 8 9 0 0 3 0 0 0.6 81 2 75 79 76 62.6 10 11 59 7 49 34 20 33.8 12 33 100 1 51 9 38.8 13 68 100 100 64 62 78.8 14 0 2 36 37 8 16.8 Prùmìr 53.4 31.2 42.1 49.9 37.2 42.8 Pro testování stochastických algoritmù globální optimalizace byla v 90. létech sestavena sada 14 obtí¾ných úloh nelineární regrese. Nìkteré z nich jsou léta v literatuøe zmiòované - Jennrich a Sampson (1968), Meyer a Roth (1972), Militký a Meloun (1994), èást z tìchto úloh pochází z nepublikovaných pøíkladù Militkého. Pøehled modelù je uveden v Køivý et al (2000), úplná data byla uvedena v Tvrdík a Køivý (1995) a v Tvrdík a Køivý (1998), v elektronické formì je mù¾ete obdr¾et na adrese
[email protected]. Na této sadì úloh byly testovány komerènì dostupné statistické pakety NCSS 2000, SYSTAT 8.0, S-PLUS 4.5 a SPSS 8.0. Pro odhad parametrù nelineárních modelù je v S-PLUS a SYSTATu u¾íván Gauss-Newtonùv algoritmus, v SYSTATu je i mo¾nost u¾ití simplexové metody, NCSS u¾ívá Levenberg-Marquartùv algoritmus a SPSS modi kovaný Levenberg-Marquartùv algoritmus. Pro ka¾dou úlohu byla
vygenerována náhodnì stovka k-tic poèateèních odhadù parametrù (k je poèet parametrù modelu, 2 7 pro úlohy z této sady). Za selhání se pova¾uje, kdy¾ program skonèí v lokálním minimu (hodnota kriteriální funkce, tj. residuální souèet ètvercù je o více jak 5 % vìt¹í ne¾ hodnota v globálním minimu). Výsledky pro statistický software dosti smutné jsou uvedeny v tabulce 3 (Krpec, 1999, Køivý et al, 2000). Globální optimalizace multimodálních funkcí je tì¾ký problém, není znám algoritmus, který by tento problém obecnì øe¹il v polynomiálním èase. Jak v¹ak ukazuje tabulka 4, lze u¾ít spolehlivìj¹í algoritmy ne¾ ty, které jsou bì¾nì implementovány ve statistickém software. V tabulce 4 jsou uvedena procenta selhání dvou stochastických algoritmù (MCRS, ES2, Køivý et al, 2000). Kromì metody nejmen¹ích ètvercù (sloupec RSS) byly jako kriteriální funkce u¾ity i nejmen¹í uøezáváné ètverce (sloupec LTS) a souèet absolutních odchylek (SAD). V tabulce jsou uvedeny jen modely, u kterých procento selhání pøi testování bylo nenulové. Tabulka 4: Procento selhání stochastických algoritmù MCRS ES2 Model RSS LTS SAD RSS LTS SAD 2 1 88 0 0 97 0 5 0 3 0 0 0 0 8 0 100 0 0 85 0 11 24 19 20 0 0 5 13 0 20 0 0 0 0
Testy shody funkcí pøe¾ití v NCSS
Pozoruhodnì podivné výsledky NCSS 2000 pøi testech shody funkcí pøe¾ití byly nedávno objeveny shodou náhod. Pan primáø Vodváøka z radioterapie FNsP v Ostravì potøeboval jen "takovou drobnost na poèkání\, tak¾e jsme nìkolik hodin u poèítaèe pøeskupovali data a chrlili funkce pøe¾ití a výsledky testù jejich shody. On s neutuchající pozorností nahlí¾el na výsledky na obrazovce a neunikla mu následující nesrovnalost ve výstupu z programu: :::
Gehans-Wilcoxon Section: : : : Chi Square = 0.63 DF = 2
Prob>CS = 0.730450
Peto-Wilcoxon Section: : : : Chi Square = 0.63 DF = 2 Log-Rank Section: : : : Chi Square = 1.76 DF = 2
Prob>CS = 0.000000 Prob>CS = 0.415603
Na první pohled je zøejmé, ¾e hodnota P = 0:000000 u Peto-Wilcoxonovy statistiky 22 = 0:63 je nesprávná. Otázkou je, zda je dobøe spoèítaná hodnota statistiky. Porovnání s výsledky získanými jinými pakety v¹ak pøineslo dal¹í pochybnosti, viz tab. 5. Hodnoty statistik shodné u S-Plus a STATA se li¹í od NCSS. Má se snad pøi statistické analýze dat u¾ívat v¾dy více programù a o správném výsledku má rozhodovat vìt¹inová shoda? Tabulka 5: Hodnoty statistik 2 soubor1 soubor2 soubor3 NCSS 2000 6.53 8.26 0.63 8.02 0.62 Wilcoxon STATA 6.0 7.80 S-Plus 4.5 7.8 8.0 0.6 NCSS 2000 6.08 11.32 1.76 Log-Rank STATA 6.0 10.15 8.47 1.35 S-Plus 4.5 10.1 8.5 1.3 Reklamoval jsem zji¹tìné nesrovnalosti v NCSS pøes dodavatele tohoto software (Statistical Solutions, Cork) u výrobce. Netu¹il jsem, ¾e se tím stávám podezøelým a budu se muset po tøi mìsíce obhajovat. První reakce J. Hintze byla, ¾e rozdílné výsledky byly získány na rùzných datech. Reklamované výsledky byly toti¾ spoèítány jako podskupiny jednoho souboru pomocí funkce FILTER a zøejmì ani autor NCSS nevìøí ve spolehlivost její implementace. Dal¹í reklamace tentokrát u¾ s daty rozdìlenými do více souborù (výsledky byly shodné s pøedchozími) pøinesla jediný pozitivní výsledek celého dlouhého reklamaèního procesu. J. Hintze pøipustil, ¾e hodnota P = 0:000000 u Peto-Wilcoxonovy statistiky 22 = 0:63 je chyba NCSS a py¹nì sdìlil, ¾e ji xoval\: Je zpùsobena nulovým poètem cenzorovaných pozorování. K "neshodám v hodnotách statistik oznámil, ¾e rozdíly ve statistikách Wilcoxonova typu jsou zpùsobeny odli¹nými variantami tìchto statistik v rùzných programech (co¾ jsem akceptoval, i kdy¾ manuály zmínìných statistických paketù jednoznaènou odpovìï nedávají), kromì toho oznámil, ¾e NCSS byl znovu
provìøen na pøíkladech z knihy Lee (1992) a bylo shledáno v¹e v poøádku. O rozdílech v log-rank testu pomlèel a tak to zùstalo. Asi jsme se dostali do nekoneèného cyklu, na otázku, proè se li¹í NCSS v log-rank testu, v¾dy pøi¹la pøedchozí odpovìï. Po dvou mìsících jsem rezignoval. Usoudil jsem, ¾e softwarové rmy více zajímá, zda zákazník platí, ne¾ to, zda mají chybu v programu. Èasy se mìní, pamatuji se, ¾e rychlé opravení reklamované chyby bylo pova¾ováno za nutnou ohajobu øemeslnické cti programátora a také tomu bylo vìnováno patøièné osobní úsilí.
Závìr
U¾ívání statistického software nepøiná¹í jen pohodlí, ale obèas také jistou frustraci ze zboøených pocitù jistoty a dùvìry. Upozornìní na problémy a chyby v implementacích generování náhody a jejich mo¾né dùsledky (Antoch, 1998) byly zaruèenì silnou ranou. Podobné rány nás v¹ak mohou potkávat i v situacích deterministických, kde bychom je oèekávali je¹tì ménì. Pokud laskavý ètenáø doèetl text a¾ sem, je patrnì zvìdav, zda mu bude nabídnuta nìjaká odpovìï na otázku vyslovenou v názvu èlánku. Vyèerpávající odpovìï asi neèeká, ale snad dvì mo¾nosti se nabízejí. Pro optimisty: Doveraj, no proveraj! A doufejme, ¾e dùsledná aplikace tohoto pøístupu pomù¾e zvý¹it spolehlivost statistického software, kdy¾ v 80. létech dokázala zmìnit svìt. Pro zdravì skeptické realisty: Nevìøte nièemu!
Literatura:
Antoch, J., Jak pomocí simulací dokázat nemo¾né, Informaèní Bulletin Èeské statistické spoleènosti, 9(1), 1{14, 1998 Ekblom, H., What can numerical analysis do for statistics, COMPSTAT 1994, Proceedings in Computational Statistics (eds R.Dutter and W. Grossmann), 31{45, Physica Verlag, 1994 Jennrich, R. I. and Sampson, P. F., Application of stepwise regression to non-linear estimation, Technometrics, 10(1), 63{72, 1968 Krpec, R., Optimalizace nelineárních regresních modelù, In: Sborník semináøe Moderní matematické metody v in¾enýrství, V©B-TUO, 66{ 69, 1999 Køivý, I., Tvrdík, J., Krpec, R., Stochastic algorithms in nonlinear regression, Comput. Statist. Data Anal. 33, 278{290, 2000
Lee, Elisa T., Statistical Methods for Survival Data Analysis, Second Edition, Wiley-Interscience, 1992 McCullough, B.D., The Accuracy of Mathematica 4 as a Statistical Package, Computational Statistics, 2000 (September), viz http://www.wolfram.com/news/statistics.html McCullough, B.D., Wilson, B., On the accuracy of statistical procedures in Microsoft Excel 97, Comput. Statist. Data Anal. 31, 27{37, 1999 Meyer, R. R. and Roth, P. M., Modi ed damped least squares: An algorithm for non-linear estimation, J. Inst. Math. Applics., 9, 218{233, 1972 Militký, J., Meloun, M.: Modus operandi of the least squares algorithm MINOPT. Talanta, 40(2), 269{277, 1994 NCSS 97, Statistical System for Windows, Number Cruncher Statistical Systems, Dr. Jerry Hintze, Kaysville, Utah, 1997 NIST, Statistical Reference Datasets, http://www.itl.nist.gov/div898/strd S-PLUS 4.5, Data Analysis Products Division, MathSoft, Seattle, 1998 SPSS ver. 8.0, SPSS Inc., Michigan, 1998 STATA 6.0, StataCorp. College Station, TX, 1999 SYSTAT 8.0, SYSTAT, Chicago, 1997 Tvrdík, J., Excel, statistika, lokalizace a zmatek, Informaèní Bulletin Èeské statistické spoleènosti, 9(2), 13{20, 1998 Tvrdík, J., Køivý, I., Stochastic algorithms in estimating regression parameters, in: J. Hanèlová (Ed.), Proceedings of the MME'95 Symposium, AIMES Press, Ostrava, 217{228, 1995 Tvrdík, J., Køivý, I., Evoluèní algoritmy a odhad parametrù nelineárních regresních modelù, In: Sborník konference Analýza dat'98, 56{69, Trilobyte, Pardubice, 1998 Tvrdík, J. and Køivý, I., Simple Evolutionary Heuristics for Global Optimization, Comput. Statist. Data Anal. (in SSN), 30, 345-352, 1999 Valchaøová, A., Aplikace evoluèních algoritmù v odhadech parametrù nelineárních regresních modelù, diplomová práce, Ostravská universita, Pøírodovìdecká fakulta, 2000
