Permalink
Browse files

Turns out I did need the fancier density function after all!

  • Loading branch information...
1 parent 13325b0 commit 5085ef8ffaf5472b530f646f23b52670067cbf1d @bos committed Nov 28, 2009
Showing with 11 additions and 12 deletions.
  1. +11 −12 Statistics/Distribution/Binomial.hs
@@ -86,35 +86,34 @@ floorf :: Double -> Double
floorf = fromIntegral . (floor :: Double -> Int64)
quantile :: BinomialDistribution -> Double -> Double
-quantile d@(BD n p) prob
+quantile dist@(BD n p) prob
| isNaN prob = prob
| p == 1 = n'
- | n' < 1e5 = fst (search 1 y z)
+ | n' < 1e5 = fst (search 1 y0 z0)
| otherwise = let dy = floorf (n' / 1000)
- in narrow dy (search dy y z)
+ in narrow dy (search dy y0 z0)
where q = 1 - p
n' = fromIntegral n
- y = n' `min` floorf (µ + σ * (d + γ * (d * d - 1) / 6) + 0.5)
+ y0 = n' `min` floorf (µ + σ * (d + γ * (d * d - 1) / 6) + 0.5)
where µ = n' * p
σ = sqrt (n' * p * q)
d = D.quantile standard prob
γ = (q - p) / σ
- z = cumulative d y
- search dy y z | z >= prob' = left y z
- | otherwise = right y
+ z0 = cumulative dist y0
+ search dy y1 z1 | z0 >= prob' = left y1 z1
+ | otherwise = right y1
where
prob' = prob * (1 - 64 * m_epsilon)
left y oldZ | y == 0 || z < prob' = (y, oldZ)
| otherwise = left (max 0 y') z
- where z = cumulative d y'
+ where z = cumulative dist y'
y' = y - dy
right y | y' >= n' || z >= prob' = (y', z)
| otherwise = right y'
- where z = cumulative d y'
+ where z = cumulative dist y'
y' = y + dy
- narrow dy (y,z)
- | dy <= 1 || dy' <= n'/1e15 = y
- | otherwise = narrow dy' (search dy y z)
+ narrow dy (y,z) | dy <= 1 || dy' <= n'/1e15 = y
+ | otherwise = narrow dy' (search dy y z)
where dy' = floorf (dy / 100)
mean :: BinomialDistribution -> Double

0 comments on commit 5085ef8

Please sign in to comment.