# 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
"""

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 don't 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 = range(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 isn't 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:
            f += 1./517. - 1./(312.*(v,1e38)[np.isinf(v)])
        else:
            f += 1./(191.*(v,1e38)[np.isinf(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
    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
        q0 = math.sqrt(2) * -y0 * \
             scipy.stats.t.isf((1.+p0)/2., (v,1e38)[v>1e38])
        q1 = math.sqrt(2) * -y1 * \
             scipy.stats.t.isf((1.+p1)/2., (v,1e38)[v>1e38])

        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,1e38)[v>1e38]))
        
    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.
    if A.has_key((p,v)):
        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)
        
    return math.sqrt(2) * -y * \
           scipy.stats.t.isf((1.+p)/2., (v,1e38)[v>1e38])

# 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')

    opt_func = lambda p, r, v : 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
        return 1. - fminbound(opt_func, .9, .999, args=(r,v))
    else:
        if q < _qsturng(.1, r, v):
            return .9
        elif q > _qsturng(.999, r, v):
            return .001
        return 1. - fminbound(opt_func, .1, .999, args=(r,v))

_vpsturng = np.vectorize(_psturng)
_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)
