# Copyright (c) 2011, Roger Lew [see LICENSE.txt] # This software is funded in part by NIH Grant P20 RR016454. """ Implementation of Gleason's (1999) non-iterative upper quantile studentized range approximation. According to Gleason this method should be more accurate than the AS190 FORTRAN algorithm of Lund and Lund (1983) and works from .5 <= p <= .999 (The AS190 only works from .9 <= p <= .99). It is more efficient then the Copenhaver & Holland (1988) algorithm (used by the _qtukey_ R function) although it requires storing the A table in memory. (q distribution) approximations in Python. see: Gleason, J. R. (1999). An accurate, non-iterative approximation for studentized range quantiles. Computational Statistics & Data Analysis, (31), 147-158. Gleason, J. R. (1998). A table of quantile points of the Studentized range distribution. http://www.stata.com/stb/stb46/dm64/sturng.pdf """ from statsmodels.compat.python import lrange import math import scipy.stats import numpy as np from scipy.optimize import fminbound inf = np.inf __version__ = '0.2.3' # changelog # 0.1 - initial release # 0.1.1 - vectorized # 0.2 - psturng added # 0.2.1 - T, R generation script relegated to make_tbls.py # 0.2.2 # - select_points refactored for performance to select_ps and # select_vs # - pysturng tester added. # 0.2.3 - uses np.inf and np.isinf # Gleason's table was derived using least square estimation on the tabled # r values for combinations of p and v. In total there are 206 # estimates over p-values of .5, .75, .9, .95, .975, .99, .995, # and .999, and over v (degrees of freedom) of (1) - 20, 24, 30, 40, # 60, 120, and inf. combinations with p < .95 do not have coefficients # for v = 1. Hence the parentheses. These coefficients allow us to # form f-hat. f-hat with the inverse t transform of tinv(p,v) yields # a fairly accurate estimate of the studentized range distribution # across a wide range of values. According to Gleason this method # should be more accurate than algorithm AS190 of Lund and Lund (1983) # and work across a wider range of values (The AS190 only works # from .9 <= p <= .99). R's qtukey algorithm was used to add tables # at .675, .8, and .85. These aid approximations when p < .9. # # The code that generated this table is called make_tbls.py and is # located in version control. A = {(0.1, 2.0): [-2.2485085243379075, -1.5641014278923464, 0.55942294426816752, -0.060006608853883377], (0.1, 3.0): [-2.2061105943901564, -1.8415406600571855, 0.61880788039834955, -0.062217093661209831], (0.1, 4.0): [-2.1686691786678178, -2.008196172372553, 0.65010084431947401, -0.06289005500114471], (0.1, 5.0): [-2.145077200277393, -2.112454843879346, 0.66701240582821342, -0.062993502233654797], (0.1, 6.0): [-2.0896098049743155, -2.2400004934286497, 0.70088523391700142, -0.065907568563272748], (0.1, 7.0): [-2.0689296655661584, -2.3078445479584873, 0.71577374609418909, -0.067081034249350552], (0.1, 8.0): [-2.0064956480711262, -2.437400413087452, 0.76297532367415266, -0.072805518121505458], (0.1, 9.0): [-2.3269477513436061, -2.0469494712773089, 0.60662518717720593, -0.054887108437009016], (0.1, 10.0): [-2.514024350177229, -1.8261187841127482, 0.51674358077906746, -0.044590425150963633], (0.1, 11.0): [-2.5130181309130828, -1.8371718595995694, 0.51336701694862252, -0.043761825829092445], (0.1, 12.0): [-2.5203508109278823, -1.8355687130611862, 0.5063486549107169, -0.042646205063108261], (0.1, 13.0): [-2.5142536438310477, -1.8496969402776282, 0.50616991367764153, -0.042378379905665363], (0.1, 14.0): [-2.3924634153781352, -2.013859173066078, 0.56421893251638688, -0.048716888109540266], (0.1, 15.0): [-2.3573552940582574, -2.0576676976224362, 0.57424068771143233, -0.049367487649225841], (0.1, 16.0): [-2.3046427483044871, -2.1295959138627993, 0.59778272657680553, -0.051864829216301617], (0.1, 17.0): [-2.2230551072316125, -2.2472837435427127, 0.64255758243215211, -0.057186665209197643], (0.1, 18.0): [-2.3912859179716897, -2.0350604070641269, 0.55924788749333332, -0.047729331835226464], (0.1, 19.0): [-2.4169773092220623, -2.0048217969339146, 0.54493039319748915, -0.045991241346224065], (0.1, 20.0): [-2.4264087194660751, -1.9916614057049267, 0.53583555139648154, -0.04463049934517662], (0.1, 24.0): [-2.3969903132061869, -2.0252941869225345, 0.53428382141200137, -0.043116495567779786], (0.1, 30.0): [-2.2509922780354623, -2.2309248956124894, 0.60748041324937263, -0.051427415888817322], (0.1, 40.0): [-2.1310090183854946, -2.3908466074610564, 0.65844375382323217, -0.05676653804036895], (0.1, 60.0): [-1.9240060179027036, -2.6685751031012233, 0.75678826647453024, -0.067938584352398995], (0.1, 120.0): [-1.9814895487030182, -2.5962051736978373, 0.71793969041292693, -0.063126863201511618], (0.1, inf): [-1.913410267066703, -2.6947367328724732, 0.74742335122750592, -0.06660897234304515], (0.5, 2.0): [-0.88295935738770648, -0.1083576698911433, 0.035214966839394388, -0.0028576288978276461], (0.5, 3.0): [-0.89085829205846834, -0.10255696422201063, 0.033613638666631696, -0.0027101699918520737], (0.5, 4.0): [-0.89627345339338116, -0.099072524607668286, 0.032657774808907684, -0.0026219007698204916], (0.5, 5.0): [-0.89959145511941052, -0.097272836582026817, 0.032236187675182958, -0.0025911555217019663], (0.5, 6.0): [-0.89959428735702474, -0.098176292411106647, 0.032590766960226995, -0.0026319890073613164], (0.5, 7.0): [-0.90131491102863937, -0.097135907620296544, 0.032304124993269533, -0.0026057965808244125], (0.5, 8.0): [-0.90292500599432901, -0.096047500971337962, 0.032030946615574568, -0.0025848748659053891], (0.5, 9.0): [-0.90385598607803697, -0.095390771554571888, 0.031832651111105899, -0.0025656060219315991], (0.5, 10.0): [-0.90562524936125388, -0.093954488089771915, 0.031414451048323286, -0.0025257834705432031], (0.5, 11.0): [-0.90420347371173826, -0.095851656370277288, 0.0321150356209743, -0.0026055056400093451], (0.5, 12.0): [-0.90585973471757664, -0.094449306296728028, 0.031705945923210958, -0.0025673330195780191], (0.5, 13.0): [-0.90555437067293054, -0.094792991050780248, 0.031826594964571089, -0.0025807109129488545], (0.5, 14.0): [-0.90652756604388762, -0.093792156994564738, 0.031468966328889042, -0.0025395175361083741], (0.5, 15.0): [-0.90642323700400085, -0.094173017520487984, 0.031657517378893905, -0.0025659271829033877], (0.5, 16.0): [-0.90716338636685234, -0.093785178083820434, 0.031630091949657997, -0.0025701459247416637], (0.5, 17.0): [-0.90790133816769714, -0.093001147638638884, 0.031376863944487084, -0.002545143621663892], (0.5, 18.0): [-0.9077432927051563, -0.093343516378180599, 0.031518139662395313, -0.0025613906133277178], (0.5, 19.0): [-0.90789499456490286, -0.09316964789456067, 0.031440782366342901, -0.0025498353345867453], (0.5, 20.0): [-0.90842707861030725, -0.092696016476608592, 0.031296040311388329, -0.0025346963982742186], (0.5, 24.0): [-0.9083281347135469, -0.092959308144970776, 0.031464063190077093, -0.0025611384271086285], (0.5, 30.0): [-0.90857624050016828, -0.093043139391980514, 0.031578791729341332, -0.0025766595412777147], (0.5, 40.0): [-0.91034085045438684, -0.091978035738914568, 0.031451631000052639, -0.0025791418103733297], (0.5, 60.0): [-0.91084356681030032, -0.091452675572423425, 0.031333147984820044, -0.0025669786958144843], (0.5, 120.0): [-0.90963649561463833, -0.093414563261352349, 0.032215602703677425, -0.0026704024780441257], (0.5, inf): [-0.91077157500981665, -0.092899220350334571, 0.032230422399363315, -0.0026696941964372916], (0.675, 2.0): [-0.67231521026565144, -0.097083624030663451, 0.027991378901661649, -0.0021425184069845558], (0.675, 3.0): [-0.65661724764645824, -0.08147195494632696, 0.02345732427073333, -0.0017448570400999351], (0.675, 4.0): [-0.65045677697461124, -0.071419073399450431, 0.020741962576852499, -0.0015171262565892491], (0.675, 5.0): [-0.64718875357808325, -0.064720611425218344, 0.019053450246546449, -0.0013836232986228711], (0.675, 6.0): [-0.64523003702018655, -0.059926313672731824, 0.017918997181483924, -0.0012992250285556828], (0.675, 7.0): [-0.64403313148478836, -0.056248191513784476, 0.017091446791293721, -0.0012406558789511822], (0.675, 8.0): [-0.64325095865764359, -0.053352543126426684, 0.016471879286491072, -0.0011991839050964099], (0.675, 9.0): [-0.64271152754911653, -0.051023769620449078, 0.01599799600547195, -0.0011693637984597086], (0.675, 10.0): [-0.64232244408502626, -0.049118327462884373, 0.015629704966568955, -0.0011477775513952285], (0.675, 11.0): [-0.64203897854353564, -0.047524627960277892, 0.015334801262767227, -0.0011315057284007177], (0.675, 12.0): [-0.64180344973512771, -0.046205907576003291, 0.015108290595438166, -0.0011207364514518488], (0.675, 13.0): [-0.64162086456823342, -0.045076099336874231, 0.0149226565346125, -0.0011126140690497352], (0.675, 14.0): [-0.64146906480198984, -0.044108523550512715, 0.014772954218646743, -0.0011069708562369386], (0.675, 15.0): [-0.64133915151966603, -0.043273370927039825, 0.014651691599222836, -0.0011032216539514398], (0.675, 16.0): [-0.64123237842752079, -0.042538925012463868, 0.014549992487506169, -0.0011005633864334021], (0.675, 17.0): [-0.64113034037536609, -0.041905699463005854, 0.014470805560767184, -0.0010995286436738471], (0.675, 18.0): [-0.64104137391561256, -0.041343885546229336, 0.014404563657113593, -0.0010991304223377683], (0.675, 19.0): [-0.64096064882827297, -0.04084569291139839, 0.014350159655133801, -0.0010993656711121901], (0.675, 20.0): [-0.64088647405089572, -0.040402175957178085, 0.014305769823654429, -0.0011001304776712105], (0.675, 24.0): [-0.64063763965937837, -0.039034716348048545, 0.014196703837251648, -0.0011061961945598175], (0.675, 30.0): [-0.64034987716294889, -0.037749651156941719, 0.014147040999127263, -0.0011188251352919833], (0.675, 40.0): [-0.6399990514713938, -0.036583307574857803, 0.014172070700846548, -0.0011391004138624943], (0.675, 60.0): [-0.63955586202430248, -0.035576938958184395, 0.014287299153378865, -0.0011675811805794236], (0.675, 120.0): [-0.63899242674778622, -0.034763757512388853, 0.014500726912982405, -0.0012028491454427466], (0.675, inf): [-0.63832682579247613, -0.034101476695520404, 0.014780921043580184, -0.0012366204114216408], (0.75, 2.0): [-0.60684073638504454, -0.096375192078057031, 0.026567529471304554, -0.0019963228971914488], (0.75, 3.0): [-0.57986144519102656, -0.078570292718034881, 0.021280637925009449, -0.0015329306898533772], (0.75, 4.0): [-0.56820771686193594, -0.0668113563896649, 0.018065284051059189, -0.0012641485481533648], (0.75, 5.0): [-0.56175292435740221, -0.058864526929603825, 0.016046735025708799, -0.0011052560286524044], (0.75, 6.0): [-0.55773449282066356, -0.053136923269827351, 0.014684258167069347, -0.0010042826823561605], (0.75, 7.0): [-0.55509524598867332, -0.048752649191139405, 0.013696566605823626, -0.00093482210003133898], (0.75, 8.0): [-0.55324993686191515, -0.045305558708724644, 0.012959681992062138, -0.00088583541601696021], (0.75, 9.0): [-0.55189259054026196, -0.042539819902381634, 0.012398791106424769, -0.00085083962241435827], (0.75, 10.0): [-0.55085384656956893, -0.040281425755686585, 0.01196442242722482, -0.00082560322161492677], (0.75, 11.0): [-0.55003198103541273, -0.038410176100193948, 0.011623294239447784, -0.00080732975034320073], (0.75, 12.0): [-0.54936541596319177, -0.036838543267887103, 0.011351822637895701, -0.0007940703654926442], (0.75, 13.0): [-0.54881015972753833, -0.035506710625568455, 0.011134691307865171, -0.0007846360016355809], (0.75, 14.0): [-0.54834094346071949, -0.034364790609906569, 0.010958873929274728, -0.00077796645357008291], (0.75, 15.0): [-0.54793602418304255, -0.033379237455748029, 0.010816140998057593, -0.00077344175064785099], (0.75, 16.0): [-0.54758347689728037, -0.032520569145898917, 0.010699240399358219, -0.00077050847328596678], (0.75, 17.0): [-0.54727115963795303, -0.031769277192927527, 0.010603749751170481, -0.0007688642392748113], (0.75, 18.0): [-0.54699351808826535, -0.031105476267880995, 0.010524669113016114, -0.00076810656837464093], (0.75, 19.0): [-0.54674357626419079, -0.030516967201954001, 0.010459478822937069, -0.00076808652582440037], (0.75, 20.0): [-0.54651728378950126, -0.029992319199769232, 0.010405694998386575, -0.0007686417223966138], (0.75, 24.0): [-0.54578309546828363, -0.028372628574010936, 0.010269939602271542, -0.00077427370647261838], (0.75, 30.0): [-0.54501246434397554, -0.026834887880579802, 0.010195603314317611, -0.00078648615954105515], (0.75, 40.0): [-0.54418127442022624, -0.025413224488871379, 0.010196455193836855, -0.00080610785749523739], (0.75, 60.0): [-0.543265189207915, -0.024141961069146383, 0.010285001019536088, -0.00083332193364294587], (0.75, 120.0): [-0.54224757817994806, -0.023039071833948214, 0.010463365295636302, -0.00086612828539477918], (0.75, inf): [-0.54114579815367159, -0.02206592527426093, 0.01070374099737127, -0.00089726564005122183], (0.8, 2.0): [-0.56895274046831146, -0.096326255190541957, 0.025815915364208686, -0.0019136561019354845], (0.8, 3.0): [-0.5336038380862278, -0.077585191014876181, 0.020184759265389905, -0.0014242746007323785], (0.8, 4.0): [-0.51780274285934258, -0.064987738443608709, 0.016713309796866204, -0.001135379856633562], (0.8, 5.0): [-0.50894361222268403, -0.056379186603362705, 0.014511270339773345, -0.00096225604117493205], (0.8, 6.0): [-0.50335153028630408, -0.050168860294790812, 0.01302807093593626, -0.00085269812692536306], (0.8, 7.0): [-0.49960934380896432, -0.045417333787806033, 0.011955593330247398, -0.00077759605604250882], (0.8, 8.0): [-0.49694518248979763, -0.041689151516021969, 0.011158986677273709, -0.00072497430103953366], (0.8, 9.0): [-0.4949559974898507, -0.038702217132906024, 0.010554360004521268, -0.0006875213117164109], (0.8, 10.0): [-0.49341407910162483, -0.036266788741325398, 0.010087354421936092, -0.00066060835062865602], (0.8, 11.0): [-0.49218129312493897, -0.034252403643273498, 0.0097218584838579536, -0.00064123459335201907], (0.8, 12.0): [-0.49117223957112183, -0.032563269730499021, 0.0094318583096021404, -0.00062725253852419032], (0.8, 13.0): [-0.49032781145131277, -0.031132495018324432, 0.0091999762562792898, -0.0006172944366003854], (0.8, 14.0): [-0.48961049628464259, -0.029906921170494854, 0.009012451847823854, -0.00061026211968669543], (0.8, 15.0): [-0.48899069793054922, -0.028849609914548158, 0.0088602820002619594, -0.00060548991575179055], (0.8, 16.0): [-0.48844921216636505, -0.027929790075266154, 0.00873599263877896, -0.00060242119796859379], (0.8, 17.0): [-0.48797119683309537, -0.027123634910159868, 0.0086338139869481887, -0.00060061821593399998], (0.8, 18.0): [-0.48754596864745836, -0.026411968723496961, 0.0085493196604705755, -0.00059977083160833624], (0.8, 19.0): [-0.48716341805691843, -0.025781422230819986, 0.0084796655915025769, -0.00059970031758323466], (0.8, 20.0): [-0.48681739197185547, -0.025219629852198749, 0.0084221844254287765, -0.00060023212822886711], (0.8, 24.0): [-0.48570639629281365, -0.023480608772518948, 0.008274490561114187, -0.000605681105792215], (0.8, 30.0): [-0.48455867067770253, -0.021824655071720423, 0.0081888502974720567, -0.00061762126933785633], (0.8, 40.0): [-0.48335478729267423, -0.020279958998363389, 0.0081765095914194709, -0.00063657117129829635], (0.8, 60.0): [-0.48207351944996679, -0.018875344346672228, 0.0082473997191472338, -0.00066242478479277243], (0.8, 120.0): [-0.48070356185330182, -0.017621686995755746, 0.0084009638803223801, -0.00069300383808949318], (0.8, inf): [-0.47926687718713606, -0.016476575352367202, 0.0086097059646591811, -0.00072160843492730911], (0.85, 2.0): [-0.53366806986381743, -0.098288178252723263, 0.026002333446289064, -0.0019567144268844896], (0.85, 3.0): [-0.48995919239619989, -0.077312722648418056, 0.019368984865418108, -0.0013449670192265796], (0.85, 4.0): [-0.46956079162382858, -0.063818518513946695, 0.015581608910696544, -0.0010264315084377606], (0.85, 5.0): [-0.45790853796153624, -0.054680511194530226, 0.013229852432203093, -0.00084248430847535898], (0.85, 6.0): [-0.4505070841695738, -0.048050936682873302, 0.011636407582714191, -0.00072491480033529815], (0.85, 7.0): [-0.44548337477336181, -0.042996612516383016, 0.010493052959891263, -0.00064528784792153239], (0.85, 8.0): [-0.44186624932664148, -0.039040005821657585, 0.0096479530794160544, -0.00058990874360967567], (0.85, 9.0): [-0.43914118689812259, -0.035875693030752713, 0.0090088804130628187, -0.00055071480339399694], (0.85, 10.0): [-0.43701255390953769, -0.033300997407157376, 0.0085172159355344848, -0.00052272770799695464], (0.85, 11.0): [-0.43530109064899053, -0.031174742038490313, 0.0081335619868386066, -0.00050268353809787927], (0.85, 12.0): [-0.43389220376610071, -0.02939618314990838, 0.007830626267772851, -0.00048836431712678222], (0.85, 13.0): [-0.43271026958463166, -0.027890759135246888, 0.0075886916668632936, -0.00047819339710596971], (0.85, 14.0): [-0.43170230265007209, -0.026604156062396189, 0.0073939099688705547, -0.00047109996854335419], (0.85, 15.0): [-0.43083160459377423, -0.025494228911600785, 0.0072358738657550868, -0.00046630677052262481], (0.85, 16.0): [-0.4300699280587239, -0.024529612608808794, 0.0071069227026219683, -0.00046323869860941791], (0.85, 17.0): [-0.42939734931902857, -0.023685025616054269, 0.0070011541609695891, -0.00046147954942994158], (0.85, 18.0): [-0.42879829041505324, -0.022940655682782165, 0.006914006369119409, -0.00046070877994711774], (0.85, 19.0): [-0.42826119448419875, -0.022280181781634649, 0.0068417746905826433, -0.00046066841214091982], (0.85, 20.0): [-0.42777654887094479, -0.021690909076747832, 0.0067817408643717969, -0.00046118620289068032], (0.85, 24.0): [-0.42622450033640852, -0.019869646711890065, 0.0066276799593494029, -0.00046668820637553747], (0.85, 30.0): [-0.42463810443233418, -0.018130114737381745, 0.0065344613060499164, -0.00047835583417510423], (0.85, 40.0): [-0.42299917804589382, -0.016498222901308417, 0.0065120558343578407, -0.00049656043685325469], (0.85, 60.0): [-0.42129387265810464, -0.014992121475265813, 0.0065657795990087635, -0.00052069705640687698], (0.85, 120.0): [-0.41951580476366368, -0.013615722489371183, 0.0066923911275726814, -0.00054846911649167492], (0.85, inf): [-0.41768751825428968, -0.012327525092266726, 0.0068664920569562592, -0.00057403720261753539], (0.9, 1.0): [-0.65851063279096722, -0.126716242078905, 0.036318801917603061, -0.002901283222928193], (0.9, 2.0): [-0.50391945369829139, -0.096996108021146235, 0.024726437623473398, -0.0017901399938303017], (0.9, 3.0): [-0.44799791843058734, -0.077180370333307199, 0.018584042055594469, -0.0012647038118363408], (0.9, 4.0): [-0.42164091756145167, -0.063427071006287514, 0.014732203755741392, -0.00094904174117957688], (0.9, 5.0): [-0.40686856251221754, -0.053361940054842398, 0.012041802076025801, -0.00072960198292410612], (0.9, 6.0): [-0.39669926026535285, -0.046951517438004242, 0.010546647213094956, -0.00062621198002366064], (0.9, 7.0): [-0.39006553675807426, -0.04169480606532109, 0.0093687546601737195, -0.00054648695713273862], (0.9, 8.0): [-0.38570205067061908, -0.037083910859179794, 0.0083233218526375836, -0.00047177586974035451], (0.9, 9.0): [-0.38190737267892938, -0.034004585655388865, 0.0077531991574119183, -0.00044306547308527872], (0.9, 10.0): [-0.37893272918125737, -0.031394677600916979, 0.0072596802503533536, -0.0004160518834299966], (0.9, 11.0): [-0.37692512492705132, -0.028780793403136471, 0.0066937909049060379, -0.00037420010136784526], (0.9, 12.0): [-0.37506345200129187, -0.026956483290567372, 0.0064147730707776523, -0.00036595383207062906], (0.9, 13.0): [-0.37339516122383209, -0.02543949524844704, 0.0061760656530197187, -0.00035678737379179527], (0.9, 14.0): [-0.37216979891087842, -0.02396347606956644, 0.0059263234465969641, -0.0003439784452550796], (0.9, 15.0): [-0.371209456600122, -0.022696132732654414, 0.0057521677184623147, -0.00033961108561770848], (0.9, 16.0): [-0.36958924377983338, -0.022227885445863002, 0.0057691706799383926, -0.00035042762538099682], (0.9, 17.0): [-0.36884224719083203, -0.021146977888668726, 0.0055957928269732716, -0.00034283810412697531], (0.9, 18.0): [-0.36803087186793326, -0.020337731477576542, 0.0054655378095212759, -0.00033452966946535248], (0.9, 19.0): [-0.3676700404163355, -0.019370115848857467, 0.0053249296207149655, -0.00032975528909580403], (0.9, 20.0): [-0.36642276267188811, -0.019344251412284838, 0.0054454968582897528, -0.00034868111677540948], (0.9, 24.0): [-0.36450650753755193, -0.017284255499990679, 0.0052337500059176749, -0.00034898202845747288], (0.9, 30.0): [-0.36251868940168608, -0.015358560437631397, 0.0050914299956134786, -0.00035574528891633978], (0.9, 40.0): [-0.36008886676510943, -0.014016835682905486, 0.0051930835959111514, -0.00038798316011984165], (0.9, 60.0): [-0.35825590690268061, -0.011991568926537646, 0.0050632208542414191, -0.00039090198974493085], (0.9, 120.0): [-0.35543612237284411, -0.011074403997811812, 0.0053504570752765162, -0.00043647137428074178], (0.9, inf): [-0.35311806343057167, -0.0096254020092145353, 0.0054548591208177181, -0.00045343916634968493], (0.95, 1.0): [-0.65330318136020071, -0.12638310760474375, 0.035987535130769424, -0.0028562665467665315], (0.95, 2.0): [-0.47225160417826934, -0.10182570362271424, 0.025846563499059158, -0.0019096769058043243], (0.95, 3.0): [-0.4056635555586528, -0.077067172693350297, 0.017789909647225533, -0.001182961668735774], (0.95, 4.0): [-0.37041675177340955, -0.063815687118939465, 0.014115210247737845, -0.00089996098435117598], (0.95, 5.0): [-0.35152398291152309, -0.052156502640669317, 0.010753738086401853, -0.0005986841939451575], (0.95, 6.0): [-0.33806730015201264, -0.045668399809578597, 0.0093168898952878162, -0.00051369719615782102], (0.95, 7.0): [-0.32924041072104465, -0.040019601775490091, 0.0080051199552865163, -0.00042054536135868043], (0.95, 8.0): [-0.32289030266989077, -0.035575345931670443, 0.0070509089344694669, -0.00035980773304803576], (0.95, 9.0): [-0.31767304201477375, -0.032464945930165703, 0.0064755950437272143, -0.0003316676253661824], (0.95, 10.0): [-0.31424318064708656, -0.029133461621153, 0.0057437449431074795, -0.00027894252261209191], (0.95, 11.0): [-0.31113589620384974, -0.02685115250591049, 0.0053517905282942889, -0.00026155954116874666], (0.95, 12.0): [-0.30848983612414582, -0.025043238019239168, 0.0050661675913488829, -0.00025017202909614005], (0.95, 13.0): [-0.3059212907410393, -0.023863874699213077, 0.0049618051135807322, -0.00025665425781125703], (0.95, 14.0): [-0.30449676902720035, -0.021983976741572344, 0.0045740513735751968, -0.00022881166323945914], (0.95, 15.0): [-0.30264908294481396, -0.02104880307520084, 0.0044866571614804382, -0.00023187587597844057], (0.95, 16.0): [-0.30118294463097917, -0.020160231061926728, 0.0044170780759056859, -0.00023733502359045826], (0.95, 17.0): [-0.30020013353427744, -0.018959271614471574, 0.0041925333038202285, -0.00022274025630789767], (0.95, 18.0): [-0.29857886556874402, -0.018664437456802001, 0.0042557787632833697, -0.00023758868868853716], (0.95, 19.0): [-0.29796289236978263, -0.017632218552317589, 0.0040792779937959866, -0.00022753271474613109], (0.95, 20.0): [-0.29681506554838077, -0.017302563243037392, 0.0041188426221428964, -0.00023913038468772782], (0.95, 24.0): [-0.29403146911167666, -0.015332330986025032, 0.0039292170319163728, -0.00024003445648641732], (0.95, 30.0): [-0.29080775563775879, -0.013844059210779323, 0.0039279165616059892, -0.00026085104496801666], (0.95, 40.0): [-0.28821583032805109, -0.011894686715666892, 0.0038202623278839982, -0.00026933325102031252], (0.95, 60.0): [-0.28525636737751447, -0.010235910558409797, 0.0038147029777580001, -0.00028598362144178959], (0.95, 120.0): [-0.28241065885026539, -0.0086103836327305026, 0.0038450612886908714, -0.00030206053671559411], (0.95, inf): [-0.27885570064169296, -0.0078122455524849222, 0.0041798538053623453, -0.0003469494881774609], (0.975, 1.0): [-0.65203598304297983, -0.12608944279227957, 0.035710038757117347, -0.0028116024425349053], (0.975, 2.0): [-0.46371891130382281, -0.096954458319996509, 0.023958312519912289, -0.0017124565391080503], (0.975, 3.0): [-0.38265282195259875, -0.076782539231612282, 0.017405078796142955, -0.0011610853687902553], (0.975, 4.0): [-0.34051193158878401, -0.063652342734671602, 0.013528310336964293, -0.00083644708934990761], (0.975, 5.0): [-0.31777655705536484, -0.051694686914334619, 0.010115807205265859, -0.00054517465344192009], (0.975, 6.0): [-0.30177149019958716, -0.044806697631189059, 0.008483551848413786, -0.00042827853925009264], (0.975, 7.0): [-0.29046972313293562, -0.039732822689098744, 0.007435356037378946, -0.00037562928283350671], (0.975, 8.0): [-0.28309484007368141, -0.034764904940713388, 0.0062932513694928518, -0.00029339243611357956], (0.975, 9.0): [-0.27711707948119785, -0.031210465194810709, 0.0055576244284178435, -0.00024663798208895803], (0.975, 10.0): [-0.27249203448553611, -0.028259756468251584, 0.00499112012528406, -0.00021535380417035389], (0.975, 11.0): [-0.26848515860011007, -0.026146703336893323, 0.0046557767110634073, -0.00020400628148271448], (0.975, 12.0): [-0.26499921540008192, -0.024522931106167097, 0.0044259624958665278, -0.00019855685376441687], (0.975, 13.0): [-0.2625023751891592, -0.022785875653297854, 0.004150277321193792, -0.00018801223218078264], (0.975, 14.0): [-0.26038552414321758, -0.021303509859738341, 0.0039195608280464681, -0.00017826200169385824], (0.975, 15.0): [-0.25801244886414665, -0.020505508012402567, 0.0038754868932712929, -0.00018588907991739744], (0.975, 16.0): [-0.25685316062360508, -0.018888418269740373, 0.0035453092842317293, -0.00016235770674204116], (0.975, 17.0): [-0.25501132271353549, -0.018362951972357794, 0.0035653933105288631, -0.00017470353354992729], (0.975, 18.0): [-0.25325045404452656, -0.017993537285026156, 0.0036035867405376691, -0.00018635492166426884], (0.975, 19.0): [-0.25236899494677928, -0.016948921372207198, 0.0034138931781330802, -0.00017462253414687881], (0.975, 20.0): [-0.25134498025027691, -0.016249564498874988, 0.0033197284005334333, -0.00017098091103245596], (0.975, 24.0): [-0.24768690797476625, -0.014668160763513996, 0.0032850791186852558, -0.00019013480716844995], (0.975, 30.0): [-0.24420834707522676, -0.012911171716272752, 0.0031977676700968051, -0.00020114907914487053], (0.975, 40.0): [-0.24105725356215926, -0.010836526056169627, 0.0030231303550754159, -0.00020128696343148667], (0.975, 60.0): [-0.23732082703955223, -0.0095442727157385391, 0.0031432904473555259, -0.00023062224109383941], (0.975, 120.0): [-0.23358581879594578, -0.0081281259918709343, 0.0031877298679120094, -0.00024496230446851501], (0.975, inf): [-0.23004105093119268, -0.0067112585174133573, 0.0032760251638919435, -0.00026244001319462992], (0.99, 1.0): [-0.65154119422706203, -0.1266603927572312, 0.03607480609672048, -0.0028668112687608113], (0.99, 2.0): [-0.45463403324378804, -0.098701236234527367, 0.024412715761684689, -0.0017613772919362193], (0.99, 3.0): [-0.36402060051035778, -0.079244959193729148, 0.017838124021360584, -0.00119080116484847], (0.99, 4.0): [-0.31903506063953818, -0.061060740682445241, 0.012093154962939612, -0.00067268347188443093], (0.99, 5.0): [-0.28917014580689182, -0.052940780099313689, 0.010231009146279354, -0.00057178339184615239], (0.99, 6.0): [-0.27283240161179012, -0.042505435573209085, 0.0072753401118264534, -0.00031314034710725922], (0.99, 7.0): [-0.25773968720546719, -0.039384214480463406, 0.0069120882597286867, -0.00032994068754356204], (0.99, 8.0): [-0.24913629282433833, -0.033831567178432859, 0.0055516244725724185, -0.00022570786249671376], (0.99, 9.0): [-0.24252380896373404, -0.029488280751457097, 0.0045215453527922998, -0.00014424552929022646], (0.99, 10.0): [-0.23654349556639986, -0.02705600214566789, 0.0041627255469343632, -0.00013804427029504753], (0.99, 11.0): [-0.23187404969432468, -0.024803662094970855, 0.0037885852786822475, -0.00012334999287725012], (0.99, 12.0): [-0.22749929386320905, -0.023655085290534145, 0.0037845051889055896, -0.00014785715789924055], (0.99, 13.0): [-0.22458989143485605, -0.021688394892771506, 0.0034075294601425251, -0.00012436961982044268], (0.99, 14.0): [-0.22197623872225777, -0.020188830700102918, 0.0031648685865587473, -0.00011320740119998819], (0.99, 15.0): [-0.2193924323730066, -0.019327469111698265, 0.0031295453754886576, -0.00012373072900083014], (0.99, 16.0): [-0.21739436875855705, -0.018215854969324128, 0.0029638341057222645, -0.00011714667871412003], (0.99, 17.0): [-0.21548926805467686, -0.017447822179412719, 0.0028994805120482812, -0.00012001887015183794], (0.99, 18.0): [-0.21365014687077843, -0.01688869353338961, 0.0028778031289216546, -0.00012591199104792711], (0.99, 19.0): [-0.21236653761262406, -0.016057151563612645, 0.0027571468998022017, -0.00012049196593780046], (0.99, 20.0): [-0.21092693178421842, -0.015641706950956638, 0.0027765989877361293, -0.00013084915163086915], (0.99, 24.0): [-0.20681960327410207, -0.013804298040271909, 0.0026308276736585674, -0.0001355061502101814], (0.99, 30.0): [-0.20271691131071576, -0.01206095288359876, 0.0025426138004198909, -0.00014589047959047533], (0.99, 40.0): [-0.19833098054449289, -0.010714533963740719, 0.0025985992420317597, -0.0001688279944262007], (0.99, 60.0): [-0.19406768821236584, -0.0093297106482013985, 0.0026521518387539584, -0.00018884874193665104], (0.99, 120.0): [-0.19010213174677365, -0.0075958207221300924, 0.0025660823297025633, -0.00018906475172834352], (0.99, inf): [-0.18602070255787137, -0.0062121155165363188, 0.0026328293420766593, -0.00020453366529867131], (0.995, 1.0): [-0.65135583544951825, -0.1266868999507193, 0.036067522182457165, -0.0028654516958844922], (0.995, 2.0): [-0.45229774013072793, -0.09869462954369547, 0.024381858599368908, -0.0017594734553033394], (0.995, 3.0): [-0.35935765236429706, -0.076650408326671915, 0.016823026893528978, -0.0010835134496404637], (0.995, 4.0): [-0.30704474720931169, -0.063093047731613019, 0.012771683306774929, -0.00075852491621809955], (0.995, 5.0): [-0.27582551740863454, -0.052533353137885791, 0.0097776009845174372, -0.00051338031756399129], (0.995, 6.0): [-0.25657971464398704, -0.043424914996692286, 0.0074324147435969991, -0.00034105188850494067], (0.995, 7.0): [-0.24090407819707738, -0.039591604712200287, 0.0068848429451020387, -0.00034737131709273414], (0.995, 8.0): [-0.23089540800827862, -0.034353305816361958, 0.0056009527629820111, -0.00024389336976992433], (0.995, 9.0): [-0.22322694848310584, -0.030294770709722547, 0.0046751239747245543, -0.00017437479314218922], (0.995, 10.0): [-0.21722684126671632, -0.026993563560163809, 0.0039811592710905491, -0.00013135281785826703], (0.995, 11.0): [-0.21171635822852911, -0.025156193618212551, 0.0037507759652964205, -0.00012959836685175671], (0.995, 12.0): [-0.20745332165849167, -0.023318819535607219, 0.0034935020002058903, -0.00012642826898405916], (0.995, 13.0): [-0.20426054591612508, -0.021189796175249527, 0.003031472176128759, -9.0497733877531618e-05], (0.995, 14.0): [-0.20113536905578902, -0.020011536696623061, 0.0029215880889956729, -9.571527213951222e-05], (0.995, 15.0): [-0.19855601561006403, -0.018808533734002542, 0.0027608859956002344, -9.2472995256929217e-05], (0.995, 16.0): [-0.19619157579534008, -0.017970461530551096, 0.0027113719105000371, -9.9864874982890861e-05], (0.995, 17.0): [-0.19428015140726104, -0.017009762497670704, 0.0025833389598201345, -9.6137545738061124e-05], (0.995, 18.0): [-0.19243180236773033, -0.01631617252107519, 0.0025227443561618621, -9.8067580523432881e-05], (0.995, 19.0): [-0.19061294393069844, -0.01586226613672222, 0.0025207005902641781, -0.00010466151274918466], (0.995, 20.0): [-0.18946302696580328, -0.014975796567260896, 0.0023700506576419867, -9.5507779057884629e-05], (0.995, 24.0): [-0.18444251428695257, -0.013770955893918012, 0.0024579445553339903, -0.00012688402863358003], (0.995, 30.0): [-0.18009742499570078, -0.011831341846559026, 0.0022801125189390046, -0.00012536249967254906], (0.995, 40.0): [-0.17562721880943261, -0.010157142650455463, 0.0022121943861923474, -0.000134542652873434], (0.995, 60.0): [-0.17084630673594547, -0.0090224965852754805, 0.0023435529965815565, -0.00016240306777440115], (0.995, 120.0): [-0.16648414081054147, -0.0074792163241677225, 0.0023284585524533607, -0.00017116464012147041], (0.995, inf): [-0.16213921875452461, -0.0058985998630496144, 0.0022605819363689093, -0.00016896211491119114], (0.999, 1.0): [-0.65233994072089363, -0.12579427445444219, 0.035830577995679271, -0.0028470555202945564], (0.999, 2.0): [-0.45050164311326341, -0.098294804380698292, 0.024134463919493736, -0.0017269603956852841], (0.999, 3.0): [-0.35161741499307819, -0.076801152272374273, 0.016695693063138672, -0.0010661121974071864], (0.999, 4.0): [-0.29398448788574133, -0.06277319725219685, 0.012454220010543127, -0.00072644165723402445], (0.999, 5.0): [-0.25725364564365477, -0.053463787584337355, 0.0099664236557431545, -0.00054866039388980659], (0.999, 6.0): [-0.23674225795168574, -0.040973155890031254, 0.0062599481191736696, -0.00021565734226586692], (0.999, 7.0): [-0.21840108878983297, -0.037037020271877719, 0.0055908063671900703, -0.00020238790479809623], (0.999, 8.0): [-0.2057964743918449, -0.032500885103194356, 0.0046441644585661756, -0.00014769592268680274], (0.999, 9.0): [-0.19604592954882674, -0.029166922919677936, 0.0040644333111949814, -0.00012854052861297006], (0.999, 10.0): [-0.18857328935948367, -0.026316705703161091, 0.0035897350868809275, -0.00011572282691335702], (0.999, 11.0): [-0.18207431428535406, -0.024201081944369412, 0.0031647372098056077, -8.1145935982296439e-05], (0.999, 12.0): [-0.17796358148991101, -0.021054306118620879, 0.0023968085939602055, -1.5907156771296993e-05], (0.999, 13.0): [-0.17371965962745489, -0.019577162950177709, 0.0022391783473999739, -2.0613023472812558e-05], (0.999, 14.0): [-0.16905298116759873, -0.01967115985443986, 0.0026495208325889269, -9.1074275220634073e-05], (0.999, 15.0): [-0.16635662558214312, -0.017903767183469876, 0.0022301322677100496, -5.1956773935885426e-05], (0.999, 16.0): [-0.16388776549525449, -0.016671918839902419, 0.0020365289602744382, -4.3592447599724942e-05], (0.999, 17.0): [-0.16131934177990759, -0.015998918405126326, 0.0019990454743285904, -4.8176277491327653e-05], (0.999, 18.0): [-0.15880633110376571, -0.015830715141055916, 0.0021688405343832091, -8.061825248932771e-05], (0.999, 19.0): [-0.15644841913314136, -0.015729364721105681, 0.0022981443610378136, -0.00010093672643417343], (0.999, 20.0): [-0.15516596606222705, -0.014725095968258637, 0.0021117117014292155, -8.8806880297328484e-05], (0.999, 24.0): [-0.14997437768645827, -0.012755323295476786, 0.0018871651510496939, -8.0896370662414938e-05], (0.999, 30.0): [-0.14459974882323703, -0.011247323832877647, 0.0018637400643826279, -9.6415323191606741e-05], (0.999, 40.0): [-0.13933285919392555, -0.0097151769692496587, 0.0018131251876208683, -0.00010452598991994023], (0.999, 60.0): [-0.13424555343804143, -0.0082163027951669444, 0.0017883427892173382, -0.00011415865110808405], (0.999, 120.0): [-0.12896119523040372, -0.0070426701112581112, 0.0018472364154226955, -0.00012862202979478294], (0.999, inf): [-0.12397213562666673, -0.0056901201604149998, 0.0018260689406957129, -0.00013263452567995485]} # p values that are defined in the A table p_keys = [.1,.5,.675,.75,.8,.85,.9,.95,.975,.99,.995,.999] # v values that are defined in the A table v_keys = lrange(2, 21) + [24, 30, 40, 60, 120, inf] def _isfloat(x): """ returns True if x is a float, returns False otherwise """ try: float(x) except: return False return True ##def _phi(p): ## """returns the pth quantile inverse norm""" ## return scipy.stats.norm.isf(p) def _phi( p ): # this function is faster than using scipy.stats.norm.isf(p) # but the permissity of the license is not explicitly listed. # using scipy.stats.norm.isf(p) is an acceptable alternative """ Modified from the author's original perl code (original comments follow below) by dfield@yahoo-inc.com. May 3, 2004. Lower tail quantile for standard normal distribution function. This function returns an approximation of the inverse cumulative standard normal distribution function. I.e., given P, it returns an approximation to the X satisfying P = Pr{Z <= X} where Z is a random variable from the standard normal distribution. The algorithm uses a minimax approximation by rational functions and the result has a relative error whose absolute value is less than 1.15e-9. Author: Peter John Acklam Time-stamp: 2000-07-19 18:26:14 E-mail: pjacklam@online.no WWW URL: http://home.online.no/~pjacklam """ if p <= 0 or p >= 1: # The original perl code exits here, we'll throw an exception instead raise ValueError( "Argument to ltqnorm %f must be in open interval (0,1)" % p ) # Coefficients in rational approximations. a = (-3.969683028665376e+01, 2.209460984245205e+02, \ -2.759285104469687e+02, 1.383577518672690e+02, \ -3.066479806614716e+01, 2.506628277459239e+00) b = (-5.447609879822406e+01, 1.615858368580409e+02, \ -1.556989798598866e+02, 6.680131188771972e+01, \ -1.328068155288572e+01 ) c = (-7.784894002430293e-03, -3.223964580411365e-01, \ -2.400758277161838e+00, -2.549732539343734e+00, \ 4.374664141464968e+00, 2.938163982698783e+00) d = ( 7.784695709041462e-03, 3.224671290700398e-01, \ 2.445134137142996e+00, 3.754408661907416e+00) # Define break-points. plow = 0.02425 phigh = 1 - plow # Rational approximation for lower region: if p < plow: q = math.sqrt(-2*math.log(p)) return -(((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) / \ ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1) # Rational approximation for upper region: if phigh < p: q = math.sqrt(-2*math.log(1-p)) return (((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) / \ ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1) # Rational approximation for central region: q = p - 0.5 r = q*q return -(((((a[0]*r+a[1])*r+a[2])*r+a[3])*r+a[4])*r+a[5])*q / \ (((((b[0]*r+b[1])*r+b[2])*r+b[3])*r+b[4])*r+1) def _ptransform(p): """function for p-value abcissa transformation""" return -1. / (1. + 1.5 * _phi((1. + p)/2.)) def _func(a, p, r, v): """ calculates f-hat for the coefficients in a, probability p, sample mean difference r, and degrees of freedom v. """ # eq. 2.3 f = a[0]*math.log(r-1.) + \ a[1]*math.log(r-1.)**2 + \ a[2]*math.log(r-1.)**3 + \ a[3]*math.log(r-1.)**4 # eq. 2.7 and 2.8 corrections if r == 3: f += -0.002 / (1. + 12. * _phi(p)**2) if v <= 4.364: v = v if not np.isinf(v) else 1e38 f += 1. / 517. - 1. / (312. * v) else: v = v if not np.isinf(v) else 1e38 f += 1. / (191. * v) return -f def _select_ps(p): # There are more generic ways of doing this but profiling # revealed that selecting these points is one of the slow # things that is easy to change. This is about 11 times # faster than the generic algorithm it is replacing. # # it is possible that different break points could yield # better estimates, but the function this is refactoring # just used linear distance. """returns the points to use for interpolating p""" if p >= .99: return .990, .995, .999 elif p >= .975: return .975, .990, .995 elif p >= .95: return .950, .975, .990 elif p >= .9125: return .900, .950, .975 elif p >= .875: return .850, .900, .950 elif p >= .825: return .800, .850, .900 elif p >= .7625: return .750, .800, .850 elif p >= .675: return .675, .750, .800 elif p >= .500: return .500, .675, .750 else: return .100, .500, .675 def _interpolate_p(p, r, v): """ interpolates p based on the values in the A table for the scalar value of r and the scalar value of v """ # interpolate p (v should be in table) # if .5 < p < .75 use linear interpolation in q # if p > .75 use quadratic interpolation in log(y + r/v) # by -1. / (1. + 1.5 * _phi((1. + p)/2.)) # find the 3 closest v values p0, p1, p2 = _select_ps(p) try: y0 = _func(A[(p0, v)], p0, r, v) + 1. except: print(p,r,v) raise y1 = _func(A[(p1, v)], p1, r, v) + 1. y2 = _func(A[(p2, v)], p2, r, v) + 1. y_log0 = math.log(y0 + float(r)/float(v)) y_log1 = math.log(y1 + float(r)/float(v)) y_log2 = math.log(y2 + float(r)/float(v)) # If p < .85 apply only the ordinate transformation # if p > .85 apply the ordinate and the abcissa transformation # In both cases apply quadratic interpolation if p > .85: p_t = _ptransform(p) p0_t = _ptransform(p0) p1_t = _ptransform(p1) p2_t = _ptransform(p2) # calculate derivatives for quadratic interpolation d2 = 2*((y_log2-y_log1)/(p2_t-p1_t) - \ (y_log1-y_log0)/(p1_t-p0_t))/(p2_t-p0_t) if (p2+p0)>=(p1+p1): d1 = (y_log2-y_log1)/(p2_t-p1_t) - 0.5*d2*(p2_t-p1_t) else: d1 = (y_log1-y_log0)/(p1_t-p0_t) + 0.5*d2*(p1_t-p0_t) d0 = y_log1 # interpolate value y_log = (d2/2.) * (p_t-p1_t)**2. + d1 * (p_t-p1_t) + d0 # transform back to y y = math.exp(y_log) - float(r)/float(v) elif p > .5: # calculate derivatives for quadratic interpolation d2 = 2*((y_log2-y_log1)/(p2-p1) - \ (y_log1-y_log0)/(p1-p0))/(p2-p0) if (p2+p0)>=(p1+p1): d1 = (y_log2-y_log1)/(p2-p1) - 0.5*d2*(p2-p1) else: d1 = (y_log1-y_log0)/(p1-p0) + 0.5*d2*(p1-p0) d0 = y_log1 # interpolate values y_log = (d2/2.) * (p-p1)**2. + d1 * (p-p1) + d0 # transform back to y y = math.exp(y_log) - float(r)/float(v) else: # linear interpolation in q and p v = min(v, 1e38) q0 = math.sqrt(2) * -y0 * \ scipy.stats.t.isf((1.+p0)/2., v) q1 = math.sqrt(2) * -y1 * \ scipy.stats.t.isf((1.+p1)/2., v) d1 = (q1-q0)/(p1-p0) d0 = q0 # interpolate values q = d1 * (p-p0) + d0 # transform back to y y = -q / (math.sqrt(2) * scipy.stats.t.isf((1.+p)/2., v)) return y def _select_vs(v, p): # This one is is about 30 times faster than # the generic algorithm it is replacing. """returns the points to use for interpolating v""" if v >= 120.: return 60, 120, inf elif v >= 60.: return 40, 60, 120 elif v >= 40.: return 30, 40, 60 elif v >= 30.: return 24, 30, 40 elif v >= 24.: return 20, 24, 30 elif v >= 19.5: return 19, 20, 24 if p >= .9: if v < 2.5: return 1, 2, 3 else: if v < 3.5: return 2, 3, 4 vi = int(round(v)) return vi - 1, vi, vi + 1 def _interpolate_v(p, r, v): """ interpolates v based on the values in the A table for the scalar value of r and th """ # interpolate v (p should be in table) # ordinate: y**2 # abcissa: 1./v # find the 3 closest v values # only p >= .9 have table values for 1 degree of freedom. # The boolean is used to index the tuple and append 1 when # p >= .9 v0, v1, v2 = _select_vs(v, p) # y = f - 1. y0_sq = (_func(A[(p,v0)], p, r, v0) + 1.)**2. y1_sq = (_func(A[(p,v1)], p, r, v1) + 1.)**2. y2_sq = (_func(A[(p,v2)], p, r, v2) + 1.)**2. # if v2 is inf set to a big number so interpolation # calculations will work if v2 > 1e38: v2 = 1e38 # transform v v_, v0_, v1_, v2_ = 1./v, 1./v0, 1./v1, 1./v2 # calculate derivatives for quadratic interpolation d2 = 2.*((y2_sq-y1_sq)/(v2_-v1_) - \ (y0_sq-y1_sq)/(v0_-v1_)) / (v2_-v0_) if (v2_ + v0_) >= (v1_ + v1_): d1 = (y2_sq-y1_sq) / (v2_-v1_) - 0.5*d2*(v2_-v1_) else: d1 = (y1_sq-y0_sq) / (v1_-v0_) + 0.5*d2*(v1_-v0_) d0 = y1_sq # calculate y y = math.sqrt((d2/2.)*(v_-v1_)**2. + d1*(v_-v1_)+ d0) return y def _qsturng(p, r, v): """scalar version of qsturng""" ## print 'q',p # r is interpolated through the q to y here we only need to # account for when p and/or v are not found in the table. global A, p_keys, v_keys if p < .1 or p > .999: raise ValueError('p must be between .1 and .999') if p < .9: if v < 2: raise ValueError('v must be > 2 when p < .9') else: if v < 1: raise ValueError('v must be > 1 when p >= .9') # The easy case. A tabled value is requested. #numpy 1.4.1: TypeError: unhashable type: 'numpy.ndarray' : p = float(p) if isinstance(v, np.ndarray): v = v.item() if (p,v) in A: y = _func(A[(p,v)], p, r, v) + 1. elif p not in p_keys and v not in v_keys+([],[1])[p>=.90]: # apply bilinear (quadratic) interpolation # # p0,v2 + o + p1,v2 + p2,v2 # r2 # # 1 # - (p,v) # v x # # r1 # p0,v1 + o + p1,v1 + p2,v1 # # # p0,v0 + o r0 + p1,v0 + p2,v0 # # _ptransform(p) # # (p1 and v1 may be below or above (p,v). The algorithm # works in both cases. For diagramatic simplicity it is # shown as above) # # 1. at v0, v1, and v2 use quadratic interpolation # to find r0, r1, r2 # # 2. use r0, r1, r2 and quadratic interpolaiton # to find y and (p,v) # find the 3 closest v values v0, v1, v2 = _select_vs(v, p) # find the 3 closest p values p0, p1, p2 = _select_ps(p) # calculate r0, r1, and r2 r0_sq = _interpolate_p(p, r, v0)**2 r1_sq = _interpolate_p(p, r, v1)**2 r2_sq = _interpolate_p(p, r, v2)**2 # transform v v_, v0_, v1_, v2_ = 1./v, 1./v0, 1./v1, 1./v2 # calculate derivatives for quadratic interpolation d2 = 2.*((r2_sq-r1_sq)/(v2_-v1_) - \ (r0_sq-r1_sq)/(v0_-v1_)) / (v2_-v0_) if (v2_ + v0_) >= (v1_ + v1_): d1 = (r2_sq-r1_sq) / (v2_-v1_) - 0.5*d2*(v2_-v1_) else: d1 = (r1_sq-r0_sq) / (v1_-v0_) + 0.5*d2*(v1_-v0_) d0 = r1_sq # calculate y y = math.sqrt((d2/2.)*(v_-v1_)**2. + d1*(v_-v1_)+ d0) elif v not in v_keys+([],[1])[p>=.90]: y = _interpolate_v(p, r, v) elif p not in p_keys: y = _interpolate_p(p, r, v) v = min(v, 1e38) return math.sqrt(2) * -y * scipy.stats.t.isf((1. + p) / 2., v) # make a qsturng functinon that will accept list-like objects _vqsturng = np.vectorize(_qsturng) _vqsturng.__doc__ = """vector version of qsturng""" def qsturng(p, r, v): """Approximates the quantile p for a studentized range distribution having v degrees of freedom and r samples for probability p. Parameters ---------- p : (scalar, array_like) The cumulative probability value p >= .1 and p <=.999 (values under .5 are not recommended) r : (scalar, array_like) The number of samples r >= 2 and r <= 200 (values over 200 are permitted but not recommended) v : (scalar, array_like) The sample degrees of freedom if p >= .9: v >=1 and v >= inf else: v >=2 and v >= inf Returns ------- q : (scalar, array_like) approximation of the Studentized Range """ if all(map(_isfloat, [p, r, v])): return _qsturng(p, r, v) return _vqsturng(p, r, v) ##def _qsturng0(p, r, v): #### print 'q0',p ## """ ## returns a first order approximation of q studentized range ## value. Based on Lund and Lund's 1983 based on the FORTRAN77 ## algorithm AS 190.2 Appl. Statist. (1983). ## """ ## vmax = 120. ## c = [0.8843, 0.2368, 1.214, 1.208, 1.4142] ## ## t = -_phi(.5+.5*p) ## if (v < vmax): ## t += (t**3. + t) / float(v) / 4. ## ## q = c[0] - c[1] * t ## if (v < vmax): ## q = q - c[2] / float(v) + c[3] * t / float(v) ## q = t * (q * math.log(r - 1.) + c[4]) ## ## # apply "bar napkin" correction for when p < .85 ## # this is good enough for our intended purpose ## if p < .85: ## q += math.log10(r) * 2.25 * (.85-p) ## return q def _psturng(q, r, v): """scalar version of psturng""" if q < 0.: raise ValueError('q should be >= 0') def opt_func(p, r, v): return np.squeeze(abs(_qsturng(p, r, v) - q)) if v == 1: if q < _qsturng(.9, r, 1): return .1 elif q > _qsturng(.999, r, 1): return .001 soln = 1. - fminbound(opt_func, .9, .999, args=(r,v)) return np.atleast_1d(soln) else: if q < _qsturng(.1, r, v): return .9 elif q > _qsturng(.999, r, v): return .001 soln = 1. - fminbound(opt_func, .1, .999, args=(r,v)) return np.atleast_1d(soln) def _psturng_scalar(q, r, v): return np.squeeze(_psturng(q, r, v)) _vpsturng = np.vectorize(_psturng_scalar) _vpsturng.__doc__ = """vector version of psturng""" def psturng(q, r, v): """Evaluates the probability from 0 to q for a studentized range having v degrees of freedom and r samples. Parameters ---------- q : (scalar, array_like) quantile value of Studentized Range q >= 0. r : (scalar, array_like) The number of samples r >= 2 and r <= 200 (values over 200 are permitted but not recommended) v : (scalar, array_like) The sample degrees of freedom if p >= .9: v >=1 and v >= inf else: v >=2 and v >= inf Returns ------- p : (scalar, array_like) 1. - area from zero to q under the Studentized Range distribution. When v == 1, p is bound between .001 and .1, when v > 1, p is bound between .001 and .9. Values between .5 and .9 are 1st order appoximations. """ if all(map(_isfloat, [q, r, v])): return _psturng(q, r, v) return _vpsturng(q, r, v) ##p, r, v = .9, 10, 20 ##print ##print 'p and v interpolation' ##print '\t20\t22\t24' ##print '.75',qsturng(.75, r, 20),qsturng(.75, r, 22),qsturng(.75, r, 24) ##print '.85',qsturng(.85, r, 20),qsturng(.85, r, 22),qsturng(.85, r, 24) ##print '.90',qsturng(.90, r, 20),qsturng(.90, r, 22),qsturng(.90, r, 24) ##print ##print 'p and v interpolation' ##print '\t120\t500\tinf' ##print '.950',qsturng(.95, r, 120),qsturng(.95, r, 500),qsturng(.95, r, inf) ##print '.960',qsturng(.96, r, 120),qsturng(.96, r, 500),qsturng(.96, r, inf) ##print '.975',qsturng(.975, r, 120),qsturng(.975, r, 500),qsturng(.975, r, inf) ##print ##print 'p and v interpolation' ##print '\t40\t50\t60' ##print '.950',qsturng(.95, r, 40),qsturng(.95, r, 50),qsturng(.95, r, 60) ##print '.960',qsturng(.96, r, 40),qsturng(.96, r, 50),qsturng(.96, r, 60) ##print '.975',qsturng(.975, r, 40),qsturng(.975, r, 50),qsturng(.975, r, 60) ##print ##print 'p and v interpolation' ##print '\t20\t22\t24' ##print '.50',qsturng(.5, r, 20),qsturng(.5, r, 22),qsturng(.5, r, 24) ##print '.60',qsturng(.6, r, 20),qsturng(.6, r, 22),qsturng(.6, r, 24) ##print '.75',qsturng(.75, r, 20),qsturng(.75, r, 22),qsturng(.75, r, 24)