1# Copyright (c) 2011, Roger Lew [see LICENSE.txt] 2# This software is funded in part by NIH Grant P20 RR016454. 3 4""" 5Implementation of Gleason's (1999) non-iterative upper quantile 6studentized range approximation. 7 8According to Gleason this method should be more accurate than the 9AS190 FORTRAN algorithm of Lund and Lund (1983) and works from .5 10<= p <= .999 (The AS190 only works from .9 <= p <= .99). 11 12It is more efficient then the Copenhaver & Holland (1988) algorithm 13(used by the _qtukey_ R function) although it requires storing the A 14table in memory. (q distribution) approximations in Python. 15 16see: 17 Gleason, J. R. (1999). An accurate, non-iterative approximation 18 for studentized range quantiles. Computational Statistics & 19 Data Analysis, (31), 147-158. 20 21 Gleason, J. R. (1998). A table of quantile points of the 22 Studentized range distribution. 23 http://www.stata.com/stb/stb46/dm64/sturng.pdf 24""" 25from statsmodels.compat.python import lrange 26import math 27import scipy.stats 28import numpy as np 29 30from scipy.optimize import fminbound 31 32inf = np.inf 33 34__version__ = '0.2.3' 35 36# changelog 37# 0.1 - initial release 38# 0.1.1 - vectorized 39# 0.2 - psturng added 40# 0.2.1 - T, R generation script relegated to make_tbls.py 41# 0.2.2 42# - select_points refactored for performance to select_ps and 43# select_vs 44# - pysturng tester added. 45# 0.2.3 - uses np.inf and np.isinf 46 47# Gleason's table was derived using least square estimation on the tabled 48# r values for combinations of p and v. In total there are 206 49# estimates over p-values of .5, .75, .9, .95, .975, .99, .995, 50# and .999, and over v (degrees of freedom) of (1) - 20, 24, 30, 40, 51# 60, 120, and inf. combinations with p < .95 do not have coefficients 52# for v = 1. Hence the parentheses. These coefficients allow us to 53# form f-hat. f-hat with the inverse t transform of tinv(p,v) yields 54# a fairly accurate estimate of the studentized range distribution 55# across a wide range of values. According to Gleason this method 56# should be more accurate than algorithm AS190 of Lund and Lund (1983) 57# and work across a wider range of values (The AS190 only works 58# from .9 <= p <= .99). R's qtukey algorithm was used to add tables 59# at .675, .8, and .85. These aid approximations when p < .9. 60# 61# The code that generated this table is called make_tbls.py and is 62# located in version control. 63A = {(0.1, 2.0): [-2.2485085243379075, -1.5641014278923464, 0.55942294426816752, -0.060006608853883377], 64 (0.1, 3.0): [-2.2061105943901564, -1.8415406600571855, 0.61880788039834955, -0.062217093661209831], 65 (0.1, 4.0): [-2.1686691786678178, -2.008196172372553, 0.65010084431947401, -0.06289005500114471], 66 (0.1, 5.0): [-2.145077200277393, -2.112454843879346, 0.66701240582821342, -0.062993502233654797], 67 (0.1, 6.0): [-2.0896098049743155, -2.2400004934286497, 0.70088523391700142, -0.065907568563272748], 68 (0.1, 7.0): [-2.0689296655661584, -2.3078445479584873, 0.71577374609418909, -0.067081034249350552], 69 (0.1, 8.0): [-2.0064956480711262, -2.437400413087452, 0.76297532367415266, -0.072805518121505458], 70 (0.1, 9.0): [-2.3269477513436061, -2.0469494712773089, 0.60662518717720593, -0.054887108437009016], 71 (0.1, 10.0): [-2.514024350177229, -1.8261187841127482, 0.51674358077906746, -0.044590425150963633], 72 (0.1, 11.0): [-2.5130181309130828, -1.8371718595995694, 0.51336701694862252, -0.043761825829092445], 73 (0.1, 12.0): [-2.5203508109278823, -1.8355687130611862, 0.5063486549107169, -0.042646205063108261], 74 (0.1, 13.0): [-2.5142536438310477, -1.8496969402776282, 0.50616991367764153, -0.042378379905665363], 75 (0.1, 14.0): [-2.3924634153781352, -2.013859173066078, 0.56421893251638688, -0.048716888109540266], 76 (0.1, 15.0): [-2.3573552940582574, -2.0576676976224362, 0.57424068771143233, -0.049367487649225841], 77 (0.1, 16.0): [-2.3046427483044871, -2.1295959138627993, 0.59778272657680553, -0.051864829216301617], 78 (0.1, 17.0): [-2.2230551072316125, -2.2472837435427127, 0.64255758243215211, -0.057186665209197643], 79 (0.1, 18.0): [-2.3912859179716897, -2.0350604070641269, 0.55924788749333332, -0.047729331835226464], 80 (0.1, 19.0): [-2.4169773092220623, -2.0048217969339146, 0.54493039319748915, -0.045991241346224065], 81 (0.1, 20.0): [-2.4264087194660751, -1.9916614057049267, 0.53583555139648154, -0.04463049934517662], 82 (0.1, 24.0): [-2.3969903132061869, -2.0252941869225345, 0.53428382141200137, -0.043116495567779786], 83 (0.1, 30.0): [-2.2509922780354623, -2.2309248956124894, 0.60748041324937263, -0.051427415888817322], 84 (0.1, 40.0): [-2.1310090183854946, -2.3908466074610564, 0.65844375382323217, -0.05676653804036895], 85 (0.1, 60.0): [-1.9240060179027036, -2.6685751031012233, 0.75678826647453024, -0.067938584352398995], 86 (0.1, 120.0): [-1.9814895487030182, -2.5962051736978373, 0.71793969041292693, -0.063126863201511618], 87 (0.1, inf): [-1.913410267066703, -2.6947367328724732, 0.74742335122750592, -0.06660897234304515], 88 (0.5, 2.0): [-0.88295935738770648, -0.1083576698911433, 0.035214966839394388, -0.0028576288978276461], 89 (0.5, 3.0): [-0.89085829205846834, -0.10255696422201063, 0.033613638666631696, -0.0027101699918520737], 90 (0.5, 4.0): [-0.89627345339338116, -0.099072524607668286, 0.032657774808907684, -0.0026219007698204916], 91 (0.5, 5.0): [-0.89959145511941052, -0.097272836582026817, 0.032236187675182958, -0.0025911555217019663], 92 (0.5, 6.0): [-0.89959428735702474, -0.098176292411106647, 0.032590766960226995, -0.0026319890073613164], 93 (0.5, 7.0): [-0.90131491102863937, -0.097135907620296544, 0.032304124993269533, -0.0026057965808244125], 94 (0.5, 8.0): [-0.90292500599432901, -0.096047500971337962, 0.032030946615574568, -0.0025848748659053891], 95 (0.5, 9.0): [-0.90385598607803697, -0.095390771554571888, 0.031832651111105899, -0.0025656060219315991], 96 (0.5, 10.0): [-0.90562524936125388, -0.093954488089771915, 0.031414451048323286, -0.0025257834705432031], 97 (0.5, 11.0): [-0.90420347371173826, -0.095851656370277288, 0.0321150356209743, -0.0026055056400093451], 98 (0.5, 12.0): [-0.90585973471757664, -0.094449306296728028, 0.031705945923210958, -0.0025673330195780191], 99 (0.5, 13.0): [-0.90555437067293054, -0.094792991050780248, 0.031826594964571089, -0.0025807109129488545], 100 (0.5, 14.0): [-0.90652756604388762, -0.093792156994564738, 0.031468966328889042, -0.0025395175361083741], 101 (0.5, 15.0): [-0.90642323700400085, -0.094173017520487984, 0.031657517378893905, -0.0025659271829033877], 102 (0.5, 16.0): [-0.90716338636685234, -0.093785178083820434, 0.031630091949657997, -0.0025701459247416637], 103 (0.5, 17.0): [-0.90790133816769714, -0.093001147638638884, 0.031376863944487084, -0.002545143621663892], 104 (0.5, 18.0): [-0.9077432927051563, -0.093343516378180599, 0.031518139662395313, -0.0025613906133277178], 105 (0.5, 19.0): [-0.90789499456490286, -0.09316964789456067, 0.031440782366342901, -0.0025498353345867453], 106 (0.5, 20.0): [-0.90842707861030725, -0.092696016476608592, 0.031296040311388329, -0.0025346963982742186], 107 (0.5, 24.0): [-0.9083281347135469, -0.092959308144970776, 0.031464063190077093, -0.0025611384271086285], 108 (0.5, 30.0): [-0.90857624050016828, -0.093043139391980514, 0.031578791729341332, -0.0025766595412777147], 109 (0.5, 40.0): [-0.91034085045438684, -0.091978035738914568, 0.031451631000052639, -0.0025791418103733297], 110 (0.5, 60.0): [-0.91084356681030032, -0.091452675572423425, 0.031333147984820044, -0.0025669786958144843], 111 (0.5, 120.0): [-0.90963649561463833, -0.093414563261352349, 0.032215602703677425, -0.0026704024780441257], 112 (0.5, inf): [-0.91077157500981665, -0.092899220350334571, 0.032230422399363315, -0.0026696941964372916], 113 (0.675, 2.0): [-0.67231521026565144, -0.097083624030663451, 0.027991378901661649, -0.0021425184069845558], 114 (0.675, 3.0): [-0.65661724764645824, -0.08147195494632696, 0.02345732427073333, -0.0017448570400999351], 115 (0.675, 4.0): [-0.65045677697461124, -0.071419073399450431, 0.020741962576852499, -0.0015171262565892491], 116 (0.675, 5.0): [-0.64718875357808325, -0.064720611425218344, 0.019053450246546449, -0.0013836232986228711], 117 (0.675, 6.0): [-0.64523003702018655, -0.059926313672731824, 0.017918997181483924, -0.0012992250285556828], 118 (0.675, 7.0): [-0.64403313148478836, -0.056248191513784476, 0.017091446791293721, -0.0012406558789511822], 119 (0.675, 8.0): [-0.64325095865764359, -0.053352543126426684, 0.016471879286491072, -0.0011991839050964099], 120 (0.675, 9.0): [-0.64271152754911653, -0.051023769620449078, 0.01599799600547195, -0.0011693637984597086], 121 (0.675, 10.0): [-0.64232244408502626, -0.049118327462884373, 0.015629704966568955, -0.0011477775513952285], 122 (0.675, 11.0): [-0.64203897854353564, -0.047524627960277892, 0.015334801262767227, -0.0011315057284007177], 123 (0.675, 12.0): [-0.64180344973512771, -0.046205907576003291, 0.015108290595438166, -0.0011207364514518488], 124 (0.675, 13.0): [-0.64162086456823342, -0.045076099336874231, 0.0149226565346125, -0.0011126140690497352], 125 (0.675, 14.0): [-0.64146906480198984, -0.044108523550512715, 0.014772954218646743, -0.0011069708562369386], 126 (0.675, 15.0): [-0.64133915151966603, -0.043273370927039825, 0.014651691599222836, -0.0011032216539514398], 127 (0.675, 16.0): [-0.64123237842752079, -0.042538925012463868, 0.014549992487506169, -0.0011005633864334021], 128 (0.675, 17.0): [-0.64113034037536609, -0.041905699463005854, 0.014470805560767184, -0.0010995286436738471], 129 (0.675, 18.0): [-0.64104137391561256, -0.041343885546229336, 0.014404563657113593, -0.0010991304223377683], 130 (0.675, 19.0): [-0.64096064882827297, -0.04084569291139839, 0.014350159655133801, -0.0010993656711121901], 131 (0.675, 20.0): [-0.64088647405089572, -0.040402175957178085, 0.014305769823654429, -0.0011001304776712105], 132 (0.675, 24.0): [-0.64063763965937837, -0.039034716348048545, 0.014196703837251648, -0.0011061961945598175], 133 (0.675, 30.0): [-0.64034987716294889, -0.037749651156941719, 0.014147040999127263, -0.0011188251352919833], 134 (0.675, 40.0): [-0.6399990514713938, -0.036583307574857803, 0.014172070700846548, -0.0011391004138624943], 135 (0.675, 60.0): [-0.63955586202430248, -0.035576938958184395, 0.014287299153378865, -0.0011675811805794236], 136 (0.675, 120.0): [-0.63899242674778622, -0.034763757512388853, 0.014500726912982405, -0.0012028491454427466], 137 (0.675, inf): [-0.63832682579247613, -0.034101476695520404, 0.014780921043580184, -0.0012366204114216408], 138 (0.75, 2.0): [-0.60684073638504454, -0.096375192078057031, 0.026567529471304554, -0.0019963228971914488], 139 (0.75, 3.0): [-0.57986144519102656, -0.078570292718034881, 0.021280637925009449, -0.0015329306898533772], 140 (0.75, 4.0): [-0.56820771686193594, -0.0668113563896649, 0.018065284051059189, -0.0012641485481533648], 141 (0.75, 5.0): [-0.56175292435740221, -0.058864526929603825, 0.016046735025708799, -0.0011052560286524044], 142 (0.75, 6.0): [-0.55773449282066356, -0.053136923269827351, 0.014684258167069347, -0.0010042826823561605], 143 (0.75, 7.0): [-0.55509524598867332, -0.048752649191139405, 0.013696566605823626, -0.00093482210003133898], 144 (0.75, 8.0): [-0.55324993686191515, -0.045305558708724644, 0.012959681992062138, -0.00088583541601696021], 145 (0.75, 9.0): [-0.55189259054026196, -0.042539819902381634, 0.012398791106424769, -0.00085083962241435827], 146 (0.75, 10.0): [-0.55085384656956893, -0.040281425755686585, 0.01196442242722482, -0.00082560322161492677], 147 (0.75, 11.0): [-0.55003198103541273, -0.038410176100193948, 0.011623294239447784, -0.00080732975034320073], 148 (0.75, 12.0): [-0.54936541596319177, -0.036838543267887103, 0.011351822637895701, -0.0007940703654926442], 149 (0.75, 13.0): [-0.54881015972753833, -0.035506710625568455, 0.011134691307865171, -0.0007846360016355809], 150 (0.75, 14.0): [-0.54834094346071949, -0.034364790609906569, 0.010958873929274728, -0.00077796645357008291], 151 (0.75, 15.0): [-0.54793602418304255, -0.033379237455748029, 0.010816140998057593, -0.00077344175064785099], 152 (0.75, 16.0): [-0.54758347689728037, -0.032520569145898917, 0.010699240399358219, -0.00077050847328596678], 153 (0.75, 17.0): [-0.54727115963795303, -0.031769277192927527, 0.010603749751170481, -0.0007688642392748113], 154 (0.75, 18.0): [-0.54699351808826535, -0.031105476267880995, 0.010524669113016114, -0.00076810656837464093], 155 (0.75, 19.0): [-0.54674357626419079, -0.030516967201954001, 0.010459478822937069, -0.00076808652582440037], 156 (0.75, 20.0): [-0.54651728378950126, -0.029992319199769232, 0.010405694998386575, -0.0007686417223966138], 157 (0.75, 24.0): [-0.54578309546828363, -0.028372628574010936, 0.010269939602271542, -0.00077427370647261838], 158 (0.75, 30.0): [-0.54501246434397554, -0.026834887880579802, 0.010195603314317611, -0.00078648615954105515], 159 (0.75, 40.0): [-0.54418127442022624, -0.025413224488871379, 0.010196455193836855, -0.00080610785749523739], 160 (0.75, 60.0): [-0.543265189207915, -0.024141961069146383, 0.010285001019536088, -0.00083332193364294587], 161 (0.75, 120.0): [-0.54224757817994806, -0.023039071833948214, 0.010463365295636302, -0.00086612828539477918], 162 (0.75, inf): [-0.54114579815367159, -0.02206592527426093, 0.01070374099737127, -0.00089726564005122183], 163 (0.8, 2.0): [-0.56895274046831146, -0.096326255190541957, 0.025815915364208686, -0.0019136561019354845], 164 (0.8, 3.0): [-0.5336038380862278, -0.077585191014876181, 0.020184759265389905, -0.0014242746007323785], 165 (0.8, 4.0): [-0.51780274285934258, -0.064987738443608709, 0.016713309796866204, -0.001135379856633562], 166 (0.8, 5.0): [-0.50894361222268403, -0.056379186603362705, 0.014511270339773345, -0.00096225604117493205], 167 (0.8, 6.0): [-0.50335153028630408, -0.050168860294790812, 0.01302807093593626, -0.00085269812692536306], 168 (0.8, 7.0): [-0.49960934380896432, -0.045417333787806033, 0.011955593330247398, -0.00077759605604250882], 169 (0.8, 8.0): [-0.49694518248979763, -0.041689151516021969, 0.011158986677273709, -0.00072497430103953366], 170 (0.8, 9.0): [-0.4949559974898507, -0.038702217132906024, 0.010554360004521268, -0.0006875213117164109], 171 (0.8, 10.0): [-0.49341407910162483, -0.036266788741325398, 0.010087354421936092, -0.00066060835062865602], 172 (0.8, 11.0): [-0.49218129312493897, -0.034252403643273498, 0.0097218584838579536, -0.00064123459335201907], 173 (0.8, 12.0): [-0.49117223957112183, -0.032563269730499021, 0.0094318583096021404, -0.00062725253852419032], 174 (0.8, 13.0): [-0.49032781145131277, -0.031132495018324432, 0.0091999762562792898, -0.0006172944366003854], 175 (0.8, 14.0): [-0.48961049628464259, -0.029906921170494854, 0.009012451847823854, -0.00061026211968669543], 176 (0.8, 15.0): [-0.48899069793054922, -0.028849609914548158, 0.0088602820002619594, -0.00060548991575179055], 177 (0.8, 16.0): [-0.48844921216636505, -0.027929790075266154, 0.00873599263877896, -0.00060242119796859379], 178 (0.8, 17.0): [-0.48797119683309537, -0.027123634910159868, 0.0086338139869481887, -0.00060061821593399998], 179 (0.8, 18.0): [-0.48754596864745836, -0.026411968723496961, 0.0085493196604705755, -0.00059977083160833624], 180 (0.8, 19.0): [-0.48716341805691843, -0.025781422230819986, 0.0084796655915025769, -0.00059970031758323466], 181 (0.8, 20.0): [-0.48681739197185547, -0.025219629852198749, 0.0084221844254287765, -0.00060023212822886711], 182 (0.8, 24.0): [-0.48570639629281365, -0.023480608772518948, 0.008274490561114187, -0.000605681105792215], 183 (0.8, 30.0): [-0.48455867067770253, -0.021824655071720423, 0.0081888502974720567, -0.00061762126933785633], 184 (0.8, 40.0): [-0.48335478729267423, -0.020279958998363389, 0.0081765095914194709, -0.00063657117129829635], 185 (0.8, 60.0): [-0.48207351944996679, -0.018875344346672228, 0.0082473997191472338, -0.00066242478479277243], 186 (0.8, 120.0): [-0.48070356185330182, -0.017621686995755746, 0.0084009638803223801, -0.00069300383808949318], 187 (0.8, inf): [-0.47926687718713606, -0.016476575352367202, 0.0086097059646591811, -0.00072160843492730911], 188 (0.85, 2.0): [-0.53366806986381743, -0.098288178252723263, 0.026002333446289064, -0.0019567144268844896], 189 (0.85, 3.0): [-0.48995919239619989, -0.077312722648418056, 0.019368984865418108, -0.0013449670192265796], 190 (0.85, 4.0): [-0.46956079162382858, -0.063818518513946695, 0.015581608910696544, -0.0010264315084377606], 191 (0.85, 5.0): [-0.45790853796153624, -0.054680511194530226, 0.013229852432203093, -0.00084248430847535898], 192 (0.85, 6.0): [-0.4505070841695738, -0.048050936682873302, 0.011636407582714191, -0.00072491480033529815], 193 (0.85, 7.0): [-0.44548337477336181, -0.042996612516383016, 0.010493052959891263, -0.00064528784792153239], 194 (0.85, 8.0): [-0.44186624932664148, -0.039040005821657585, 0.0096479530794160544, -0.00058990874360967567], 195 (0.85, 9.0): [-0.43914118689812259, -0.035875693030752713, 0.0090088804130628187, -0.00055071480339399694], 196 (0.85, 10.0): [-0.43701255390953769, -0.033300997407157376, 0.0085172159355344848, -0.00052272770799695464], 197 (0.85, 11.0): [-0.43530109064899053, -0.031174742038490313, 0.0081335619868386066, -0.00050268353809787927], 198 (0.85, 12.0): [-0.43389220376610071, -0.02939618314990838, 0.007830626267772851, -0.00048836431712678222], 199 (0.85, 13.0): [-0.43271026958463166, -0.027890759135246888, 0.0075886916668632936, -0.00047819339710596971], 200 (0.85, 14.0): [-0.43170230265007209, -0.026604156062396189, 0.0073939099688705547, -0.00047109996854335419], 201 (0.85, 15.0): [-0.43083160459377423, -0.025494228911600785, 0.0072358738657550868, -0.00046630677052262481], 202 (0.85, 16.0): [-0.4300699280587239, -0.024529612608808794, 0.0071069227026219683, -0.00046323869860941791], 203 (0.85, 17.0): [-0.42939734931902857, -0.023685025616054269, 0.0070011541609695891, -0.00046147954942994158], 204 (0.85, 18.0): [-0.42879829041505324, -0.022940655682782165, 0.006914006369119409, -0.00046070877994711774], 205 (0.85, 19.0): [-0.42826119448419875, -0.022280181781634649, 0.0068417746905826433, -0.00046066841214091982], 206 (0.85, 20.0): [-0.42777654887094479, -0.021690909076747832, 0.0067817408643717969, -0.00046118620289068032], 207 (0.85, 24.0): [-0.42622450033640852, -0.019869646711890065, 0.0066276799593494029, -0.00046668820637553747], 208 (0.85, 30.0): [-0.42463810443233418, -0.018130114737381745, 0.0065344613060499164, -0.00047835583417510423], 209 (0.85, 40.0): [-0.42299917804589382, -0.016498222901308417, 0.0065120558343578407, -0.00049656043685325469], 210 (0.85, 60.0): [-0.42129387265810464, -0.014992121475265813, 0.0065657795990087635, -0.00052069705640687698], 211 (0.85, 120.0): [-0.41951580476366368, -0.013615722489371183, 0.0066923911275726814, -0.00054846911649167492], 212 (0.85, inf): [-0.41768751825428968, -0.012327525092266726, 0.0068664920569562592, -0.00057403720261753539], 213 (0.9, 1.0): [-0.65851063279096722, -0.126716242078905, 0.036318801917603061, -0.002901283222928193], 214 (0.9, 2.0): [-0.50391945369829139, -0.096996108021146235, 0.024726437623473398, -0.0017901399938303017], 215 (0.9, 3.0): [-0.44799791843058734, -0.077180370333307199, 0.018584042055594469, -0.0012647038118363408], 216 (0.9, 4.0): [-0.42164091756145167, -0.063427071006287514, 0.014732203755741392, -0.00094904174117957688], 217 (0.9, 5.0): [-0.40686856251221754, -0.053361940054842398, 0.012041802076025801, -0.00072960198292410612], 218 (0.9, 6.0): [-0.39669926026535285, -0.046951517438004242, 0.010546647213094956, -0.00062621198002366064], 219 (0.9, 7.0): [-0.39006553675807426, -0.04169480606532109, 0.0093687546601737195, -0.00054648695713273862], 220 (0.9, 8.0): [-0.38570205067061908, -0.037083910859179794, 0.0083233218526375836, -0.00047177586974035451], 221 (0.9, 9.0): [-0.38190737267892938, -0.034004585655388865, 0.0077531991574119183, -0.00044306547308527872], 222 (0.9, 10.0): [-0.37893272918125737, -0.031394677600916979, 0.0072596802503533536, -0.0004160518834299966], 223 (0.9, 11.0): [-0.37692512492705132, -0.028780793403136471, 0.0066937909049060379, -0.00037420010136784526], 224 (0.9, 12.0): [-0.37506345200129187, -0.026956483290567372, 0.0064147730707776523, -0.00036595383207062906], 225 (0.9, 13.0): [-0.37339516122383209, -0.02543949524844704, 0.0061760656530197187, -0.00035678737379179527], 226 (0.9, 14.0): [-0.37216979891087842, -0.02396347606956644, 0.0059263234465969641, -0.0003439784452550796], 227 (0.9, 15.0): [-0.371209456600122, -0.022696132732654414, 0.0057521677184623147, -0.00033961108561770848], 228 (0.9, 16.0): [-0.36958924377983338, -0.022227885445863002, 0.0057691706799383926, -0.00035042762538099682], 229 (0.9, 17.0): [-0.36884224719083203, -0.021146977888668726, 0.0055957928269732716, -0.00034283810412697531], 230 (0.9, 18.0): [-0.36803087186793326, -0.020337731477576542, 0.0054655378095212759, -0.00033452966946535248], 231 (0.9, 19.0): [-0.3676700404163355, -0.019370115848857467, 0.0053249296207149655, -0.00032975528909580403], 232 (0.9, 20.0): [-0.36642276267188811, -0.019344251412284838, 0.0054454968582897528, -0.00034868111677540948], 233 (0.9, 24.0): [-0.36450650753755193, -0.017284255499990679, 0.0052337500059176749, -0.00034898202845747288], 234 (0.9, 30.0): [-0.36251868940168608, -0.015358560437631397, 0.0050914299956134786, -0.00035574528891633978], 235 (0.9, 40.0): [-0.36008886676510943, -0.014016835682905486, 0.0051930835959111514, -0.00038798316011984165], 236 (0.9, 60.0): [-0.35825590690268061, -0.011991568926537646, 0.0050632208542414191, -0.00039090198974493085], 237 (0.9, 120.0): [-0.35543612237284411, -0.011074403997811812, 0.0053504570752765162, -0.00043647137428074178], 238 (0.9, inf): [-0.35311806343057167, -0.0096254020092145353, 0.0054548591208177181, -0.00045343916634968493], 239 (0.95, 1.0): [-0.65330318136020071, -0.12638310760474375, 0.035987535130769424, -0.0028562665467665315], 240 (0.95, 2.0): [-0.47225160417826934, -0.10182570362271424, 0.025846563499059158, -0.0019096769058043243], 241 (0.95, 3.0): [-0.4056635555586528, -0.077067172693350297, 0.017789909647225533, -0.001182961668735774], 242 (0.95, 4.0): [-0.37041675177340955, -0.063815687118939465, 0.014115210247737845, -0.00089996098435117598], 243 (0.95, 5.0): [-0.35152398291152309, -0.052156502640669317, 0.010753738086401853, -0.0005986841939451575], 244 (0.95, 6.0): [-0.33806730015201264, -0.045668399809578597, 0.0093168898952878162, -0.00051369719615782102], 245 (0.95, 7.0): [-0.32924041072104465, -0.040019601775490091, 0.0080051199552865163, -0.00042054536135868043], 246 (0.95, 8.0): [-0.32289030266989077, -0.035575345931670443, 0.0070509089344694669, -0.00035980773304803576], 247 (0.95, 9.0): [-0.31767304201477375, -0.032464945930165703, 0.0064755950437272143, -0.0003316676253661824], 248 (0.95, 10.0): [-0.31424318064708656, -0.029133461621153, 0.0057437449431074795, -0.00027894252261209191], 249 (0.95, 11.0): [-0.31113589620384974, -0.02685115250591049, 0.0053517905282942889, -0.00026155954116874666], 250 (0.95, 12.0): [-0.30848983612414582, -0.025043238019239168, 0.0050661675913488829, -0.00025017202909614005], 251 (0.95, 13.0): [-0.3059212907410393, -0.023863874699213077, 0.0049618051135807322, -0.00025665425781125703], 252 (0.95, 14.0): [-0.30449676902720035, -0.021983976741572344, 0.0045740513735751968, -0.00022881166323945914], 253 (0.95, 15.0): [-0.30264908294481396, -0.02104880307520084, 0.0044866571614804382, -0.00023187587597844057], 254 (0.95, 16.0): [-0.30118294463097917, -0.020160231061926728, 0.0044170780759056859, -0.00023733502359045826], 255 (0.95, 17.0): [-0.30020013353427744, -0.018959271614471574, 0.0041925333038202285, -0.00022274025630789767], 256 (0.95, 18.0): [-0.29857886556874402, -0.018664437456802001, 0.0042557787632833697, -0.00023758868868853716], 257 (0.95, 19.0): [-0.29796289236978263, -0.017632218552317589, 0.0040792779937959866, -0.00022753271474613109], 258 (0.95, 20.0): [-0.29681506554838077, -0.017302563243037392, 0.0041188426221428964, -0.00023913038468772782], 259 (0.95, 24.0): [-0.29403146911167666, -0.015332330986025032, 0.0039292170319163728, -0.00024003445648641732], 260 (0.95, 30.0): [-0.29080775563775879, -0.013844059210779323, 0.0039279165616059892, -0.00026085104496801666], 261 (0.95, 40.0): [-0.28821583032805109, -0.011894686715666892, 0.0038202623278839982, -0.00026933325102031252], 262 (0.95, 60.0): [-0.28525636737751447, -0.010235910558409797, 0.0038147029777580001, -0.00028598362144178959], 263 (0.95, 120.0): [-0.28241065885026539, -0.0086103836327305026, 0.0038450612886908714, -0.00030206053671559411], 264 (0.95, inf): [-0.27885570064169296, -0.0078122455524849222, 0.0041798538053623453, -0.0003469494881774609], 265 (0.975, 1.0): [-0.65203598304297983, -0.12608944279227957, 0.035710038757117347, -0.0028116024425349053], 266 (0.975, 2.0): [-0.46371891130382281, -0.096954458319996509, 0.023958312519912289, -0.0017124565391080503], 267 (0.975, 3.0): [-0.38265282195259875, -0.076782539231612282, 0.017405078796142955, -0.0011610853687902553], 268 (0.975, 4.0): [-0.34051193158878401, -0.063652342734671602, 0.013528310336964293, -0.00083644708934990761], 269 (0.975, 5.0): [-0.31777655705536484, -0.051694686914334619, 0.010115807205265859, -0.00054517465344192009], 270 (0.975, 6.0): [-0.30177149019958716, -0.044806697631189059, 0.008483551848413786, -0.00042827853925009264], 271 (0.975, 7.0): [-0.29046972313293562, -0.039732822689098744, 0.007435356037378946, -0.00037562928283350671], 272 (0.975, 8.0): [-0.28309484007368141, -0.034764904940713388, 0.0062932513694928518, -0.00029339243611357956], 273 (0.975, 9.0): [-0.27711707948119785, -0.031210465194810709, 0.0055576244284178435, -0.00024663798208895803], 274 (0.975, 10.0): [-0.27249203448553611, -0.028259756468251584, 0.00499112012528406, -0.00021535380417035389], 275 (0.975, 11.0): [-0.26848515860011007, -0.026146703336893323, 0.0046557767110634073, -0.00020400628148271448], 276 (0.975, 12.0): [-0.26499921540008192, -0.024522931106167097, 0.0044259624958665278, -0.00019855685376441687], 277 (0.975, 13.0): [-0.2625023751891592, -0.022785875653297854, 0.004150277321193792, -0.00018801223218078264], 278 (0.975, 14.0): [-0.26038552414321758, -0.021303509859738341, 0.0039195608280464681, -0.00017826200169385824], 279 (0.975, 15.0): [-0.25801244886414665, -0.020505508012402567, 0.0038754868932712929, -0.00018588907991739744], 280 (0.975, 16.0): [-0.25685316062360508, -0.018888418269740373, 0.0035453092842317293, -0.00016235770674204116], 281 (0.975, 17.0): [-0.25501132271353549, -0.018362951972357794, 0.0035653933105288631, -0.00017470353354992729], 282 (0.975, 18.0): [-0.25325045404452656, -0.017993537285026156, 0.0036035867405376691, -0.00018635492166426884], 283 (0.975, 19.0): [-0.25236899494677928, -0.016948921372207198, 0.0034138931781330802, -0.00017462253414687881], 284 (0.975, 20.0): [-0.25134498025027691, -0.016249564498874988, 0.0033197284005334333, -0.00017098091103245596], 285 (0.975, 24.0): [-0.24768690797476625, -0.014668160763513996, 0.0032850791186852558, -0.00019013480716844995], 286 (0.975, 30.0): [-0.24420834707522676, -0.012911171716272752, 0.0031977676700968051, -0.00020114907914487053], 287 (0.975, 40.0): [-0.24105725356215926, -0.010836526056169627, 0.0030231303550754159, -0.00020128696343148667], 288 (0.975, 60.0): [-0.23732082703955223, -0.0095442727157385391, 0.0031432904473555259, -0.00023062224109383941], 289 (0.975, 120.0): [-0.23358581879594578, -0.0081281259918709343, 0.0031877298679120094, -0.00024496230446851501], 290 (0.975, inf): [-0.23004105093119268, -0.0067112585174133573, 0.0032760251638919435, -0.00026244001319462992], 291 (0.99, 1.0): [-0.65154119422706203, -0.1266603927572312, 0.03607480609672048, -0.0028668112687608113], 292 (0.99, 2.0): [-0.45463403324378804, -0.098701236234527367, 0.024412715761684689, -0.0017613772919362193], 293 (0.99, 3.0): [-0.36402060051035778, -0.079244959193729148, 0.017838124021360584, -0.00119080116484847], 294 (0.99, 4.0): [-0.31903506063953818, -0.061060740682445241, 0.012093154962939612, -0.00067268347188443093], 295 (0.99, 5.0): [-0.28917014580689182, -0.052940780099313689, 0.010231009146279354, -0.00057178339184615239], 296 (0.99, 6.0): [-0.27283240161179012, -0.042505435573209085, 0.0072753401118264534, -0.00031314034710725922], 297 (0.99, 7.0): [-0.25773968720546719, -0.039384214480463406, 0.0069120882597286867, -0.00032994068754356204], 298 (0.99, 8.0): [-0.24913629282433833, -0.033831567178432859, 0.0055516244725724185, -0.00022570786249671376], 299 (0.99, 9.0): [-0.24252380896373404, -0.029488280751457097, 0.0045215453527922998, -0.00014424552929022646], 300 (0.99, 10.0): [-0.23654349556639986, -0.02705600214566789, 0.0041627255469343632, -0.00013804427029504753], 301 (0.99, 11.0): [-0.23187404969432468, -0.024803662094970855, 0.0037885852786822475, -0.00012334999287725012], 302 (0.99, 12.0): [-0.22749929386320905, -0.023655085290534145, 0.0037845051889055896, -0.00014785715789924055], 303 (0.99, 13.0): [-0.22458989143485605, -0.021688394892771506, 0.0034075294601425251, -0.00012436961982044268], 304 (0.99, 14.0): [-0.22197623872225777, -0.020188830700102918, 0.0031648685865587473, -0.00011320740119998819], 305 (0.99, 15.0): [-0.2193924323730066, -0.019327469111698265, 0.0031295453754886576, -0.00012373072900083014], 306 (0.99, 16.0): [-0.21739436875855705, -0.018215854969324128, 0.0029638341057222645, -0.00011714667871412003], 307 (0.99, 17.0): [-0.21548926805467686, -0.017447822179412719, 0.0028994805120482812, -0.00012001887015183794], 308 (0.99, 18.0): [-0.21365014687077843, -0.01688869353338961, 0.0028778031289216546, -0.00012591199104792711], 309 (0.99, 19.0): [-0.21236653761262406, -0.016057151563612645, 0.0027571468998022017, -0.00012049196593780046], 310 (0.99, 20.0): [-0.21092693178421842, -0.015641706950956638, 0.0027765989877361293, -0.00013084915163086915], 311 (0.99, 24.0): [-0.20681960327410207, -0.013804298040271909, 0.0026308276736585674, -0.0001355061502101814], 312 (0.99, 30.0): [-0.20271691131071576, -0.01206095288359876, 0.0025426138004198909, -0.00014589047959047533], 313 (0.99, 40.0): [-0.19833098054449289, -0.010714533963740719, 0.0025985992420317597, -0.0001688279944262007], 314 (0.99, 60.0): [-0.19406768821236584, -0.0093297106482013985, 0.0026521518387539584, -0.00018884874193665104], 315 (0.99, 120.0): [-0.19010213174677365, -0.0075958207221300924, 0.0025660823297025633, -0.00018906475172834352], 316 (0.99, inf): [-0.18602070255787137, -0.0062121155165363188, 0.0026328293420766593, -0.00020453366529867131], 317 (0.995, 1.0): [-0.65135583544951825, -0.1266868999507193, 0.036067522182457165, -0.0028654516958844922], 318 (0.995, 2.0): [-0.45229774013072793, -0.09869462954369547, 0.024381858599368908, -0.0017594734553033394], 319 (0.995, 3.0): [-0.35935765236429706, -0.076650408326671915, 0.016823026893528978, -0.0010835134496404637], 320 (0.995, 4.0): [-0.30704474720931169, -0.063093047731613019, 0.012771683306774929, -0.00075852491621809955], 321 (0.995, 5.0): [-0.27582551740863454, -0.052533353137885791, 0.0097776009845174372, -0.00051338031756399129], 322 (0.995, 6.0): [-0.25657971464398704, -0.043424914996692286, 0.0074324147435969991, -0.00034105188850494067], 323 (0.995, 7.0): [-0.24090407819707738, -0.039591604712200287, 0.0068848429451020387, -0.00034737131709273414], 324 (0.995, 8.0): [-0.23089540800827862, -0.034353305816361958, 0.0056009527629820111, -0.00024389336976992433], 325 (0.995, 9.0): [-0.22322694848310584, -0.030294770709722547, 0.0046751239747245543, -0.00017437479314218922], 326 (0.995, 10.0): [-0.21722684126671632, -0.026993563560163809, 0.0039811592710905491, -0.00013135281785826703], 327 (0.995, 11.0): [-0.21171635822852911, -0.025156193618212551, 0.0037507759652964205, -0.00012959836685175671], 328 (0.995, 12.0): [-0.20745332165849167, -0.023318819535607219, 0.0034935020002058903, -0.00012642826898405916], 329 (0.995, 13.0): [-0.20426054591612508, -0.021189796175249527, 0.003031472176128759, -9.0497733877531618e-05], 330 (0.995, 14.0): [-0.20113536905578902, -0.020011536696623061, 0.0029215880889956729, -9.571527213951222e-05], 331 (0.995, 15.0): [-0.19855601561006403, -0.018808533734002542, 0.0027608859956002344, -9.2472995256929217e-05], 332 (0.995, 16.0): [-0.19619157579534008, -0.017970461530551096, 0.0027113719105000371, -9.9864874982890861e-05], 333 (0.995, 17.0): [-0.19428015140726104, -0.017009762497670704, 0.0025833389598201345, -9.6137545738061124e-05], 334 (0.995, 18.0): [-0.19243180236773033, -0.01631617252107519, 0.0025227443561618621, -9.8067580523432881e-05], 335 (0.995, 19.0): [-0.19061294393069844, -0.01586226613672222, 0.0025207005902641781, -0.00010466151274918466], 336 (0.995, 20.0): [-0.18946302696580328, -0.014975796567260896, 0.0023700506576419867, -9.5507779057884629e-05], 337 (0.995, 24.0): [-0.18444251428695257, -0.013770955893918012, 0.0024579445553339903, -0.00012688402863358003], 338 (0.995, 30.0): [-0.18009742499570078, -0.011831341846559026, 0.0022801125189390046, -0.00012536249967254906], 339 (0.995, 40.0): [-0.17562721880943261, -0.010157142650455463, 0.0022121943861923474, -0.000134542652873434], 340 (0.995, 60.0): [-0.17084630673594547, -0.0090224965852754805, 0.0023435529965815565, -0.00016240306777440115], 341 (0.995, 120.0): [-0.16648414081054147, -0.0074792163241677225, 0.0023284585524533607, -0.00017116464012147041], 342 (0.995, inf): [-0.16213921875452461, -0.0058985998630496144, 0.0022605819363689093, -0.00016896211491119114], 343 (0.999, 1.0): [-0.65233994072089363, -0.12579427445444219, 0.035830577995679271, -0.0028470555202945564], 344 (0.999, 2.0): [-0.45050164311326341, -0.098294804380698292, 0.024134463919493736, -0.0017269603956852841], 345 (0.999, 3.0): [-0.35161741499307819, -0.076801152272374273, 0.016695693063138672, -0.0010661121974071864], 346 (0.999, 4.0): [-0.29398448788574133, -0.06277319725219685, 0.012454220010543127, -0.00072644165723402445], 347 (0.999, 5.0): [-0.25725364564365477, -0.053463787584337355, 0.0099664236557431545, -0.00054866039388980659], 348 (0.999, 6.0): [-0.23674225795168574, -0.040973155890031254, 0.0062599481191736696, -0.00021565734226586692], 349 (0.999, 7.0): [-0.21840108878983297, -0.037037020271877719, 0.0055908063671900703, -0.00020238790479809623], 350 (0.999, 8.0): [-0.2057964743918449, -0.032500885103194356, 0.0046441644585661756, -0.00014769592268680274], 351 (0.999, 9.0): [-0.19604592954882674, -0.029166922919677936, 0.0040644333111949814, -0.00012854052861297006], 352 (0.999, 10.0): [-0.18857328935948367, -0.026316705703161091, 0.0035897350868809275, -0.00011572282691335702], 353 (0.999, 11.0): [-0.18207431428535406, -0.024201081944369412, 0.0031647372098056077, -8.1145935982296439e-05], 354 (0.999, 12.0): [-0.17796358148991101, -0.021054306118620879, 0.0023968085939602055, -1.5907156771296993e-05], 355 (0.999, 13.0): [-0.17371965962745489, -0.019577162950177709, 0.0022391783473999739, -2.0613023472812558e-05], 356 (0.999, 14.0): [-0.16905298116759873, -0.01967115985443986, 0.0026495208325889269, -9.1074275220634073e-05], 357 (0.999, 15.0): [-0.16635662558214312, -0.017903767183469876, 0.0022301322677100496, -5.1956773935885426e-05], 358 (0.999, 16.0): [-0.16388776549525449, -0.016671918839902419, 0.0020365289602744382, -4.3592447599724942e-05], 359 (0.999, 17.0): [-0.16131934177990759, -0.015998918405126326, 0.0019990454743285904, -4.8176277491327653e-05], 360 (0.999, 18.0): [-0.15880633110376571, -0.015830715141055916, 0.0021688405343832091, -8.061825248932771e-05], 361 (0.999, 19.0): [-0.15644841913314136, -0.015729364721105681, 0.0022981443610378136, -0.00010093672643417343], 362 (0.999, 20.0): [-0.15516596606222705, -0.014725095968258637, 0.0021117117014292155, -8.8806880297328484e-05], 363 (0.999, 24.0): [-0.14997437768645827, -0.012755323295476786, 0.0018871651510496939, -8.0896370662414938e-05], 364 (0.999, 30.0): [-0.14459974882323703, -0.011247323832877647, 0.0018637400643826279, -9.6415323191606741e-05], 365 (0.999, 40.0): [-0.13933285919392555, -0.0097151769692496587, 0.0018131251876208683, -0.00010452598991994023], 366 (0.999, 60.0): [-0.13424555343804143, -0.0082163027951669444, 0.0017883427892173382, -0.00011415865110808405], 367 (0.999, 120.0): [-0.12896119523040372, -0.0070426701112581112, 0.0018472364154226955, -0.00012862202979478294], 368 (0.999, inf): [-0.12397213562666673, -0.0056901201604149998, 0.0018260689406957129, -0.00013263452567995485]} 369 370# p values that are defined in the A table 371p_keys = [.1,.5,.675,.75,.8,.85,.9,.95,.975,.99,.995,.999] 372 373# v values that are defined in the A table 374v_keys = lrange(2, 21) + [24, 30, 40, 60, 120, inf] 375 376def _isfloat(x): 377 """ 378 returns True if x is a float, 379 returns False otherwise 380 """ 381 try: 382 float(x) 383 except: 384 return False 385 386 return True 387 388##def _phi(p): 389## """returns the pth quantile inverse norm""" 390## return scipy.stats.norm.isf(p) 391 392def _phi( p ): 393 # this function is faster than using scipy.stats.norm.isf(p) 394 # but the permissity of the license is not explicitly listed. 395 # using scipy.stats.norm.isf(p) is an acceptable alternative 396 """ 397 Modified from the author's original perl code (original comments follow below) 398 by dfield@yahoo-inc.com. May 3, 2004. 399 400 Lower tail quantile for standard normal distribution function. 401 402 This function returns an approximation of the inverse cumulative 403 standard normal distribution function. I.e., given P, it returns 404 an approximation to the X satisfying P = Pr{Z <= X} where Z is a 405 random variable from the standard normal distribution. 406 407 The algorithm uses a minimax approximation by rational functions 408 and the result has a relative error whose absolute value is less 409 than 1.15e-9. 410 411 Author: Peter John Acklam 412 Time-stamp: 2000-07-19 18:26:14 413 E-mail: pjacklam@online.no 414 WWW URL: http://home.online.no/~pjacklam 415 """ 416 417 if p <= 0 or p >= 1: 418 # The original perl code exits here, we'll throw an exception instead 419 raise ValueError( "Argument to ltqnorm %f must be in open interval (0,1)" % p ) 420 421 # Coefficients in rational approximations. 422 a = (-3.969683028665376e+01, 2.209460984245205e+02, \ 423 -2.759285104469687e+02, 1.383577518672690e+02, \ 424 -3.066479806614716e+01, 2.506628277459239e+00) 425 b = (-5.447609879822406e+01, 1.615858368580409e+02, \ 426 -1.556989798598866e+02, 6.680131188771972e+01, \ 427 -1.328068155288572e+01 ) 428 c = (-7.784894002430293e-03, -3.223964580411365e-01, \ 429 -2.400758277161838e+00, -2.549732539343734e+00, \ 430 4.374664141464968e+00, 2.938163982698783e+00) 431 d = ( 7.784695709041462e-03, 3.224671290700398e-01, \ 432 2.445134137142996e+00, 3.754408661907416e+00) 433 434 # Define break-points. 435 plow = 0.02425 436 phigh = 1 - plow 437 438 # Rational approximation for lower region: 439 if p < plow: 440 q = math.sqrt(-2*math.log(p)) 441 return -(((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) / \ 442 ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1) 443 444 # Rational approximation for upper region: 445 if phigh < p: 446 q = math.sqrt(-2*math.log(1-p)) 447 return (((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) / \ 448 ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1) 449 450 # Rational approximation for central region: 451 q = p - 0.5 452 r = q*q 453 return -(((((a[0]*r+a[1])*r+a[2])*r+a[3])*r+a[4])*r+a[5])*q / \ 454 (((((b[0]*r+b[1])*r+b[2])*r+b[3])*r+b[4])*r+1) 455 456def _ptransform(p): 457 """function for p-value abcissa transformation""" 458 return -1. / (1. + 1.5 * _phi((1. + p)/2.)) 459 460def _func(a, p, r, v): 461 """ 462 calculates f-hat for the coefficients in a, probability p, 463 sample mean difference r, and degrees of freedom v. 464 """ 465 # eq. 2.3 466 f = a[0]*math.log(r-1.) + \ 467 a[1]*math.log(r-1.)**2 + \ 468 a[2]*math.log(r-1.)**3 + \ 469 a[3]*math.log(r-1.)**4 470 471 # eq. 2.7 and 2.8 corrections 472 if r == 3: 473 f += -0.002 / (1. + 12. * _phi(p)**2) 474 475 if v <= 4.364: 476 v = v if not np.isinf(v) else 1e38 477 f += 1. / 517. - 1. / (312. * v) 478 else: 479 v = v if not np.isinf(v) else 1e38 480 f += 1. / (191. * v) 481 482 return -f 483 484def _select_ps(p): 485 # There are more generic ways of doing this but profiling 486 # revealed that selecting these points is one of the slow 487 # things that is easy to change. This is about 11 times 488 # faster than the generic algorithm it is replacing. 489 # 490 # it is possible that different break points could yield 491 # better estimates, but the function this is refactoring 492 # just used linear distance. 493 """returns the points to use for interpolating p""" 494 if p >= .99: 495 return .990, .995, .999 496 elif p >= .975: 497 return .975, .990, .995 498 elif p >= .95: 499 return .950, .975, .990 500 elif p >= .9125: 501 return .900, .950, .975 502 elif p >= .875: 503 return .850, .900, .950 504 elif p >= .825: 505 return .800, .850, .900 506 elif p >= .7625: 507 return .750, .800, .850 508 elif p >= .675: 509 return .675, .750, .800 510 elif p >= .500: 511 return .500, .675, .750 512 else: 513 return .100, .500, .675 514 515def _interpolate_p(p, r, v): 516 """ 517 interpolates p based on the values in the A table for the 518 scalar value of r and the scalar value of v 519 """ 520 521 # interpolate p (v should be in table) 522 # if .5 < p < .75 use linear interpolation in q 523 # if p > .75 use quadratic interpolation in log(y + r/v) 524 # by -1. / (1. + 1.5 * _phi((1. + p)/2.)) 525 526 # find the 3 closest v values 527 p0, p1, p2 = _select_ps(p) 528 try: 529 y0 = _func(A[(p0, v)], p0, r, v) + 1. 530 except: 531 print(p,r,v) 532 raise 533 y1 = _func(A[(p1, v)], p1, r, v) + 1. 534 y2 = _func(A[(p2, v)], p2, r, v) + 1. 535 536 y_log0 = math.log(y0 + float(r)/float(v)) 537 y_log1 = math.log(y1 + float(r)/float(v)) 538 y_log2 = math.log(y2 + float(r)/float(v)) 539 540 # If p < .85 apply only the ordinate transformation 541 # if p > .85 apply the ordinate and the abcissa transformation 542 # In both cases apply quadratic interpolation 543 if p > .85: 544 p_t = _ptransform(p) 545 p0_t = _ptransform(p0) 546 p1_t = _ptransform(p1) 547 p2_t = _ptransform(p2) 548 549 # calculate derivatives for quadratic interpolation 550 d2 = 2*((y_log2-y_log1)/(p2_t-p1_t) - \ 551 (y_log1-y_log0)/(p1_t-p0_t))/(p2_t-p0_t) 552 if (p2+p0)>=(p1+p1): 553 d1 = (y_log2-y_log1)/(p2_t-p1_t) - 0.5*d2*(p2_t-p1_t) 554 else: 555 d1 = (y_log1-y_log0)/(p1_t-p0_t) + 0.5*d2*(p1_t-p0_t) 556 d0 = y_log1 557 558 # interpolate value 559 y_log = (d2/2.) * (p_t-p1_t)**2. + d1 * (p_t-p1_t) + d0 560 561 # transform back to y 562 y = math.exp(y_log) - float(r)/float(v) 563 564 elif p > .5: 565 # calculate derivatives for quadratic interpolation 566 d2 = 2*((y_log2-y_log1)/(p2-p1) - \ 567 (y_log1-y_log0)/(p1-p0))/(p2-p0) 568 if (p2+p0)>=(p1+p1): 569 d1 = (y_log2-y_log1)/(p2-p1) - 0.5*d2*(p2-p1) 570 else: 571 d1 = (y_log1-y_log0)/(p1-p0) + 0.5*d2*(p1-p0) 572 d0 = y_log1 573 574 # interpolate values 575 y_log = (d2/2.) * (p-p1)**2. + d1 * (p-p1) + d0 576 577 # transform back to y 578 y = math.exp(y_log) - float(r)/float(v) 579 580 else: 581 # linear interpolation in q and p 582 v = min(v, 1e38) 583 q0 = math.sqrt(2) * -y0 * \ 584 scipy.stats.t.isf((1.+p0)/2., v) 585 q1 = math.sqrt(2) * -y1 * \ 586 scipy.stats.t.isf((1.+p1)/2., v) 587 588 d1 = (q1-q0)/(p1-p0) 589 d0 = q0 590 591 # interpolate values 592 q = d1 * (p-p0) + d0 593 594 # transform back to y 595 y = -q / (math.sqrt(2) * scipy.stats.t.isf((1.+p)/2., v)) 596 597 return y 598 599def _select_vs(v, p): 600 # This one is is about 30 times faster than 601 # the generic algorithm it is replacing. 602 """returns the points to use for interpolating v""" 603 604 if v >= 120.: 605 return 60, 120, inf 606 elif v >= 60.: 607 return 40, 60, 120 608 elif v >= 40.: 609 return 30, 40, 60 610 elif v >= 30.: 611 return 24, 30, 40 612 elif v >= 24.: 613 return 20, 24, 30 614 elif v >= 19.5: 615 return 19, 20, 24 616 617 if p >= .9: 618 if v < 2.5: 619 return 1, 2, 3 620 else: 621 if v < 3.5: 622 return 2, 3, 4 623 624 vi = int(round(v)) 625 return vi - 1, vi, vi + 1 626 627def _interpolate_v(p, r, v): 628 """ 629 interpolates v based on the values in the A table for the 630 scalar value of r and th 631 """ 632 # interpolate v (p should be in table) 633 # ordinate: y**2 634 # abcissa: 1./v 635 636 # find the 3 closest v values 637 # only p >= .9 have table values for 1 degree of freedom. 638 # The boolean is used to index the tuple and append 1 when 639 # p >= .9 640 v0, v1, v2 = _select_vs(v, p) 641 642 # y = f - 1. 643 y0_sq = (_func(A[(p,v0)], p, r, v0) + 1.)**2. 644 y1_sq = (_func(A[(p,v1)], p, r, v1) + 1.)**2. 645 y2_sq = (_func(A[(p,v2)], p, r, v2) + 1.)**2. 646 647 # if v2 is inf set to a big number so interpolation 648 # calculations will work 649 if v2 > 1e38: 650 v2 = 1e38 651 652 # transform v 653 v_, v0_, v1_, v2_ = 1./v, 1./v0, 1./v1, 1./v2 654 655 # calculate derivatives for quadratic interpolation 656 d2 = 2.*((y2_sq-y1_sq)/(v2_-v1_) - \ 657 (y0_sq-y1_sq)/(v0_-v1_)) / (v2_-v0_) 658 if (v2_ + v0_) >= (v1_ + v1_): 659 d1 = (y2_sq-y1_sq) / (v2_-v1_) - 0.5*d2*(v2_-v1_) 660 else: 661 d1 = (y1_sq-y0_sq) / (v1_-v0_) + 0.5*d2*(v1_-v0_) 662 d0 = y1_sq 663 664 # calculate y 665 y = math.sqrt((d2/2.)*(v_-v1_)**2. + d1*(v_-v1_)+ d0) 666 667 return y 668 669def _qsturng(p, r, v): 670 """scalar version of qsturng""" 671## print 'q',p 672 # r is interpolated through the q to y here we only need to 673 # account for when p and/or v are not found in the table. 674 global A, p_keys, v_keys 675 676 if p < .1 or p > .999: 677 raise ValueError('p must be between .1 and .999') 678 679 if p < .9: 680 if v < 2: 681 raise ValueError('v must be > 2 when p < .9') 682 else: 683 if v < 1: 684 raise ValueError('v must be > 1 when p >= .9') 685 686 # The easy case. A tabled value is requested. 687 688 #numpy 1.4.1: TypeError: unhashable type: 'numpy.ndarray' : 689 p = float(p) 690 if isinstance(v, np.ndarray): 691 v = v.item() 692 if (p,v) in A: 693 y = _func(A[(p,v)], p, r, v) + 1. 694 695 elif p not in p_keys and v not in v_keys+([],[1])[p>=.90]: 696 # apply bilinear (quadratic) interpolation 697 # 698 # p0,v2 + o + p1,v2 + p2,v2 699 # r2 700 # 701 # 1 702 # - (p,v) 703 # v x 704 # 705 # r1 706 # p0,v1 + o + p1,v1 + p2,v1 707 # 708 # 709 # p0,v0 + o r0 + p1,v0 + p2,v0 710 # 711 # _ptransform(p) 712 # 713 # (p1 and v1 may be below or above (p,v). The algorithm 714 # works in both cases. For diagramatic simplicity it is 715 # shown as above) 716 # 717 # 1. at v0, v1, and v2 use quadratic interpolation 718 # to find r0, r1, r2 719 # 720 # 2. use r0, r1, r2 and quadratic interpolaiton 721 # to find y and (p,v) 722 723 # find the 3 closest v values 724 v0, v1, v2 = _select_vs(v, p) 725 726 # find the 3 closest p values 727 p0, p1, p2 = _select_ps(p) 728 729 # calculate r0, r1, and r2 730 r0_sq = _interpolate_p(p, r, v0)**2 731 r1_sq = _interpolate_p(p, r, v1)**2 732 r2_sq = _interpolate_p(p, r, v2)**2 733 734 # transform v 735 v_, v0_, v1_, v2_ = 1./v, 1./v0, 1./v1, 1./v2 736 737 # calculate derivatives for quadratic interpolation 738 d2 = 2.*((r2_sq-r1_sq)/(v2_-v1_) - \ 739 (r0_sq-r1_sq)/(v0_-v1_)) / (v2_-v0_) 740 if (v2_ + v0_) >= (v1_ + v1_): 741 d1 = (r2_sq-r1_sq) / (v2_-v1_) - 0.5*d2*(v2_-v1_) 742 else: 743 d1 = (r1_sq-r0_sq) / (v1_-v0_) + 0.5*d2*(v1_-v0_) 744 d0 = r1_sq 745 746 # calculate y 747 y = math.sqrt((d2/2.)*(v_-v1_)**2. + d1*(v_-v1_)+ d0) 748 749 elif v not in v_keys+([],[1])[p>=.90]: 750 y = _interpolate_v(p, r, v) 751 752 elif p not in p_keys: 753 y = _interpolate_p(p, r, v) 754 755 v = min(v, 1e38) 756 return math.sqrt(2) * -y * scipy.stats.t.isf((1. + p) / 2., v) 757 758# make a qsturng functinon that will accept list-like objects 759_vqsturng = np.vectorize(_qsturng) 760_vqsturng.__doc__ = """vector version of qsturng""" 761 762def qsturng(p, r, v): 763 """Approximates the quantile p for a studentized range 764 distribution having v degrees of freedom and r samples 765 for probability p. 766 767 Parameters 768 ---------- 769 p : (scalar, array_like) 770 The cumulative probability value 771 p >= .1 and p <=.999 772 (values under .5 are not recommended) 773 r : (scalar, array_like) 774 The number of samples 775 r >= 2 and r <= 200 776 (values over 200 are permitted but not recommended) 777 v : (scalar, array_like) 778 The sample degrees of freedom 779 if p >= .9: 780 v >=1 and v >= inf 781 else: 782 v >=2 and v >= inf 783 784 Returns 785 ------- 786 q : (scalar, array_like) 787 approximation of the Studentized Range 788 """ 789 790 if all(map(_isfloat, [p, r, v])): 791 return _qsturng(p, r, v) 792 return _vqsturng(p, r, v) 793 794##def _qsturng0(p, r, v): 795#### print 'q0',p 796## """ 797## returns a first order approximation of q studentized range 798## value. Based on Lund and Lund's 1983 based on the FORTRAN77 799## algorithm AS 190.2 Appl. Statist. (1983). 800## """ 801## vmax = 120. 802## c = [0.8843, 0.2368, 1.214, 1.208, 1.4142] 803## 804## t = -_phi(.5+.5*p) 805## if (v < vmax): 806## t += (t**3. + t) / float(v) / 4. 807## 808## q = c[0] - c[1] * t 809## if (v < vmax): 810## q = q - c[2] / float(v) + c[3] * t / float(v) 811## q = t * (q * math.log(r - 1.) + c[4]) 812## 813## # apply "bar napkin" correction for when p < .85 814## # this is good enough for our intended purpose 815## if p < .85: 816## q += math.log10(r) * 2.25 * (.85-p) 817## return q 818 819def _psturng(q, r, v): 820 """scalar version of psturng""" 821 if q < 0.: 822 raise ValueError('q should be >= 0') 823 824 def opt_func(p, r, v): 825 return np.squeeze(abs(_qsturng(p, r, v) - q)) 826 827 if v == 1: 828 if q < _qsturng(.9, r, 1): 829 return .1 830 elif q > _qsturng(.999, r, 1): 831 return .001 832 soln = 1. - fminbound(opt_func, .9, .999, args=(r,v)) 833 return np.atleast_1d(soln) 834 else: 835 if q < _qsturng(.1, r, v): 836 return .9 837 elif q > _qsturng(.999, r, v): 838 return .001 839 soln = 1. - fminbound(opt_func, .1, .999, args=(r,v)) 840 return np.atleast_1d(soln) 841 842_vpsturng = np.vectorize(_psturng) 843_vpsturng.__doc__ = """vector version of psturng""" 844 845def psturng(q, r, v): 846 """Evaluates the probability from 0 to q for a studentized 847 range having v degrees of freedom and r samples. 848 849 Parameters 850 ---------- 851 q : (scalar, array_like) 852 quantile value of Studentized Range 853 q >= 0. 854 r : (scalar, array_like) 855 The number of samples 856 r >= 2 and r <= 200 857 (values over 200 are permitted but not recommended) 858 v : (scalar, array_like) 859 The sample degrees of freedom 860 if p >= .9: 861 v >=1 and v >= inf 862 else: 863 v >=2 and v >= inf 864 865 Returns 866 ------- 867 p : (scalar, array_like) 868 1. - area from zero to q under the Studentized Range 869 distribution. When v == 1, p is bound between .001 870 and .1, when v > 1, p is bound between .001 and .9. 871 Values between .5 and .9 are 1st order appoximations. 872 """ 873 if all(map(_isfloat, [q, r, v])): 874 return _psturng(q, r, v) 875 return _vpsturng(q, r, v) 876 877##p, r, v = .9, 10, 20 878##print 879##print 'p and v interpolation' 880##print '\t20\t22\t24' 881##print '.75',qsturng(.75, r, 20),qsturng(.75, r, 22),qsturng(.75, r, 24) 882##print '.85',qsturng(.85, r, 20),qsturng(.85, r, 22),qsturng(.85, r, 24) 883##print '.90',qsturng(.90, r, 20),qsturng(.90, r, 22),qsturng(.90, r, 24) 884##print 885##print 'p and v interpolation' 886##print '\t120\t500\tinf' 887##print '.950',qsturng(.95, r, 120),qsturng(.95, r, 500),qsturng(.95, r, inf) 888##print '.960',qsturng(.96, r, 120),qsturng(.96, r, 500),qsturng(.96, r, inf) 889##print '.975',qsturng(.975, r, 120),qsturng(.975, r, 500),qsturng(.975, r, inf) 890##print 891##print 'p and v interpolation' 892##print '\t40\t50\t60' 893##print '.950',qsturng(.95, r, 40),qsturng(.95, r, 50),qsturng(.95, r, 60) 894##print '.960',qsturng(.96, r, 40),qsturng(.96, r, 50),qsturng(.96, r, 60) 895##print '.975',qsturng(.975, r, 40),qsturng(.975, r, 50),qsturng(.975, r, 60) 896##print 897##print 'p and v interpolation' 898##print '\t20\t22\t24' 899##print '.50',qsturng(.5, r, 20),qsturng(.5, r, 22),qsturng(.5, r, 24) 900##print '.60',qsturng(.6, r, 20),qsturng(.6, r, 22),qsturng(.6, r, 24) 901##print '.75',qsturng(.75, r, 20),qsturng(.75, r, 22),qsturng(.75, r, 24) 902