# Commits

committed ce8e845

Documentation and formatting

• Participants
• Parent commits f44d6e0

# Numeric/SpecFunctions.hs

` `
` `
` -- | Compute inverse of regularized incomplete beta function. Uses`
`--- initial approximation from AS109 and Halley method to solve equation.`
`+-- initial approximation from AS109, AS64 and Halley method to solve`
`+-- equation.`
` invIncompleteBeta :: Double     -- ^ /p/ > 0`
`                   -> Double     -- ^ /q/ > 0`
`                   -> Double     -- ^ /a/ ∈ [0,1]`
`   | a <  0 || a >  1 = error "bad a"`
`   | a == 0 || a == 1 = a`
`   | a > 0.5          = 1 - invIncompleteBetaWorker (logBeta p q) q p (1 - a)`
`-  | otherwise        = invIncompleteBetaWorker (logBeta p q) p q a`
`+  | otherwise        =     invIncompleteBetaWorker (logBeta p q) p q  a`
` `
` invIncompleteBetaWorker :: Double -> Double -> Double -> Double -> Double`
` invIncompleteBetaWorker beta p q a = loop (0::Int) guess`
`             | z > 1     = (x + 1) / 2`
`             | otherwise = z`
`             where z = x - dx`
`-    -- Calculate initial guess`
`+    -- Calculate initial guess. Approximations are described in the AS64.`
`+    --`
`+    -- Note that a <= 0.5.`
`     guess `
`       | p > 1 && q > 1 = `
`           let rr = (y*y - 3) / 6`
`       where`
`         r   = sqrt \$ - 2 * log a`
`         y   = r - ( 2.30753 + 0.27061 * r )`
`-                   / ( 1.0 + ( 0.99229 + 0.04481 * r ) * r )`
`+                  / ( 1.0 + ( 0.99229 + 0.04481 * r ) * r )`
`         t   = 1 / (9 * q)`
`         t'  = 2 * q * (1 - t + y * sqrt t) ** 3`
`         t'' = (4*p + 2*q - 2) / t'`