diff --git a/hstats.cabal b/hstats.cabal
--- a/hstats.cabal
+++ b/hstats.cabal
@@ -1,15 +1,16 @@
 Name:                hstats
-Version:             0.2
+Version:             0.3
 License:             BSD3
 License-file:        LICENSE
 Author:              Marshall Beddoe
-Copyright:           Copyright (c) 2007, SFTank
+Copyright:           Copyright (c) 2008, Marshall Beddoe
 category:            Math
 synopsis:            Statistical Computing in Haskell
 description:         A library of commonly used statistical functions
-maintainer:          mbeddoe@<nospam>sftank.net
-homepage:            http://www.sftank.net
+maintainer:          mbeddoe@<nospam>gmail.com
+homepage:            http://github.com/unmarshal/hstats/ 
 hs-source-dirs:      src
-ghc-options:         -O
 exposed-Modules:     Math.Statistics
+extensions:          BangPatterns
 build-depends:       base>=2.0, haskell98
+build-type:          Simple
diff --git a/src/Math/Statistics.hs b/src/Math/Statistics.hs
--- a/src/Math/Statistics.hs
+++ b/src/Math/Statistics.hs
@@ -1,9 +1,11 @@
+{-# OPTIONS_GHC -XBangPatterns #-}
+
 -----------------------------------------------------------------------------
 -- Module      : Math.Statistics
--- Copyright   : (c) 2007 SFTank
+-- Copyright   : (c) 2008 Marshall Beddoe
 -- License     : BSD3
 --
--- Maintainer  : mbeddoe@<nospam>sftank.net
+-- Maintainer  : mbeddoe@<nospam>gmail.com
 -- Stability   : experimental
 -- Portability : portable
 --
@@ -13,21 +15,26 @@
 
 module Math.Statistics where
 
-import List
+import Data.List
+import Data.Ord (comparing)
 
--- Numerically stable mean
+-- |Numerically stable mean
 mean :: Floating a => [a] -> a
-mean = fst.foldr (\x (m,n) -> (m+(x-m) / (fromIntegral $ n + 1), n + 1)) (0,0)
+mean x = fst $ foldl' (\(!m, !n) x -> (m+(x-m)/(n+1),n+1)) (0,0) x
 
--- Harmonic mean
-hmean :: (Floating a) => [a] -> a
-hmean xs = fromIntegral (length xs) / (sum $ map (1/) xs)
+-- |Same as 'mean' 
+average :: Floating a => [a] -> a
+average = mean
 
--- Geometric mean
-gmean :: (Floating a) => [a] -> a
-gmean xs = (foldr1 (*) xs)**(1 / fromIntegral (length xs))
+-- |Harmonic mean
+harmean :: (Floating a) => [a] -> a
+harmean xs = fromIntegral (length xs) / (sum $ map (1/) xs)
 
--- Median
+-- |Geometric mean
+geomean :: (Floating a) => [a] -> a
+geomean xs = (foldr1 (*) xs)**(1 / fromIntegral (length xs))
+
+-- |Median
 median :: (Floating a, Ord a) => [a] -> a
 median x | odd n  = head  $ drop (n `div` 2) x'
          | even n = mean $ take 2 $ drop i x'
@@ -35,40 +42,48 @@
                         x' = sort x
                         n  = length x
 
--- Modes
--- Returns a sorted list of modes in descending order
+-- |Modes returns a sorted list of modes in descending order
 modes :: (Ord a) => [a] -> [(Int, a)]
-modes xs = sortOn (negate.fst) $ map (\x->(length x, head x)) $ (group.sort) xs
-    where
-      sortOn :: Ord b => (a -> b) -> [a] -> [a]
-      sortOn f = sortBy (\x y -> compare (f x) (f y))
+modes xs = sortBy (comparing $ negate.fst) $ map (\x->(length x, head x)) $ (group.sort) xs
 
--- Central moments
+-- |Mode returns the mode of the list, otherwise Nothing
+mode :: (Ord a) => [a] -> Maybe a
+mode xs = case m of
+            [] -> Nothing
+            otherwise -> Just . snd $ head m
+    where m = filter (\(a,b) -> a > 1) (modes xs)
+
+-- |Central moments
+centralMoment :: (Floating b, Integral t) => [b] -> t -> b
 centralMoment xs 1 = 0
 centralMoment xs r = (sum (map (\x -> (x-m)^r) xs)) / n
     where
       m = mean xs
       n = fromIntegral $ length xs
 
--- Range
+-- |Range
 range :: (Num a, Ord a) => [a] -> a
 range xs = maximum xs - minimum xs
 
--- Average deviation
+-- |Average deviation
 avgdev :: (Floating a) => [a] -> a
 avgdev xs = mean $ map (\x -> abs(x - m)) xs
     where
       m = mean xs
 
--- Standard deviation
+-- |Standard deviation of sample
 stddev :: (Floating a) => [a] -> a
 stddev xs = sqrt $ var xs
 
--- Population variance
+-- |Standard deviation of population
+stddevp :: (Floating a) => [a] -> a
+stddevp xs = sqrt $ pvar xs
+
+-- |Population variance
 pvar :: (Floating a) => [a] -> a
 pvar xs = centralMoment xs 2
 
--- Numerically stable sample variance
+-- |Sample variance
 var xs = (var' 0 0 0 xs) / (fromIntegral $ length xs - 1)
     where
       var' _ _ s [] = s
@@ -77,27 +92,48 @@
            delta = x - m
            nm = m + delta/(fromIntegral $ n + 1)
 
--- Interquartile range
--- XXX: Add case that takes into account even vs odd length
+-- |Interquartile range
 iqr xs = take (length xs - 2*q) $ drop q xs
     where
       q = ((length xs) + 1) `div` 4
 
 -- Kurtosis
-kurtosis xs = ((centralMoment xs 4) / (centralMoment xs 2)^2)-3
+kurt xs = ((centralMoment xs 4) / (centralMoment xs 2)^2)-3
 
--- Skew
+-- |Arbitrary quantile q of an unsorted list.  The quantile /q/ of /N/
+-- |data points is the point whose (zero-based) index in the sorted
+-- |data set is closest to /q(N-1)/.
+quantile :: (Fractional b, Ord b) => Double -> [b] -> b
+quantile q = quantileAsc q . sort
+
+-- |As 'quantile' specialized for sorted data
+quantileAsc :: (Fractional b, Ord b) => Double -> [b] -> b
+quantileAsc _ [] = error "quantile on empty list"
+quantileAsc q xs
+    | q < 0 || q > 1 = error "quantile out of range"
+    | otherwise = xs !! (quantIndex (length xs) q)
+    where quantIndex :: Int -> Double -> Int
+          quantIndex len q = case round $ q * (fromIntegral len - 1) of
+                               idx | idx < 0    -> error "Quantile index too small"
+                                   | idx >= len -> error "Quantile index too large"
+                                   | otherwise  -> idx
+
+-- |Calculate skew
+skew :: (Floating b) => [b] -> b
 skew xs = (centralMoment xs 3) / (centralMoment xs 2)**(3/2)
 
+-- |Calculates pearson skew
+pearsonSkew1 :: (Ord a, Floating a) => [a] -> a
 pearsonSkew1 xs = 3 * (mean xs - mo) / stddev xs
     where
       mo = snd $ head $ modes xs
 
+pearsonSkew2 :: (Ord a, Floating a) => [a] -> a
 pearsonSkew2 xs = 3 * (mean xs - median xs) / stddev xs
 
--- Covariance
-cov :: (Floating a) => [a] -> [a] -> a
-cov xs ys = sum (zipWith (*) (map f1 xs) (map f2 ys)) / (n - 1)
+-- |Sample Covariance
+covar :: (Floating a) => [a] -> [a] -> a
+covar xs ys = sum (zipWith (*) (map f1 xs) (map f2 ys)) / (n-1)
     where
       n = fromIntegral $ length $ xs
       m1 = mean xs
@@ -105,14 +141,41 @@
       f1 = \x -> (x - m1)
       f2 = \x -> (x - m2)
 
--- Covariance matrix
+-- |Covariance matrix
 covMatrix :: (Floating a) => [[a]] -> [[a]]
 covMatrix xs =  split' (length xs) cs
     where
-      cs = [ cov a b | a <- xs, b <- xs]
+      cs = [ covar a b | a <- xs, b <- xs]
       split' n = unfoldr (\y -> if null y then Nothing else Just $ splitAt n y)
 
--- Pearson's product-moment correlation coefficient
-corr :: (Floating a) => [a] -> [a] -> a
-corr x y = cov x y / (stddev x * stddev y)
+-- |Pearson's product-moment correlation coefficient
+pearson :: (Floating a) => [a] -> [a] -> a
+pearson x y = covar x y / (stddev x * stddev y)
 
+-- |Same as 'pearson'
+correl :: (Floating a) => [a] -> [a] -> a
+correl = pearson
+
+-- |Least-squares linear regression of /y/ against /x/ for a
+-- |collection of (/x/, /y/) data, in the form of (/b0/, /b1/, /r/)
+-- |where the regression is /y/ = /b0/ + /b1/ * /x/ with Pearson
+-- |coefficient /r/
+linreg :: (Floating b) => [(b, b)] -> (b, b, b)
+linreg xys = let !xs = map fst xys
+                 !ys = map snd xys
+                 !n = fromIntegral $ length xys
+                 !sX = sum xs
+                 !sY = sum ys
+                 !sXX = sum $ map (^ 2) xs
+                 !sXY = sum $ map (uncurry (*)) xys
+                 !sYY = sum $ map (^ 2) ys
+                 !alpha = (sY - beta * sX) / n
+                 !beta = (n * sXY - sX * sY) / (n * sXX - sX * sX)
+                 !r = (n * sXY - sX * sY) / (sqrt $ (n * sXX - sX^2) * (n * sYY - sY ^ 2))
+             in (alpha, beta, r)
+
+
+-- |Returns the sum of square deviations from their sample mean.
+devsq :: (Floating a) => [a] -> a
+devsq xs = sum $ map (\x->(x-m)**2) xs
+    where m = mean xs