diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c) 2012, Anthony Cowley
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the name of Anthony Cowley nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/RANSAC.cabal b/RANSAC.cabal
new file mode 100644
--- /dev/null
+++ b/RANSAC.cabal
@@ -0,0 +1,33 @@
+name:                RANSAC
+version:             0.1.0.0
+synopsis:            The RANSAC algorithm for parameter estimation.
+description:         The RANdom SAmple Consensus (RANSAC) algorithm for
+                     estimating the parameters of a mathematical model
+                     from a data set. See
+                     <http://en.wikipedia.org/wiki/RANSAC> for more
+                     information.
+                     .
+                     See @tests/LinearFit.hs@ in the package contents for 
+                     an example.
+license:             BSD3
+license-file:        LICENSE
+author:              Anthony Cowley
+maintainer:          acowley@gmail.com
+copyright:           (c) Anthony Cowley 2012
+category:            Math,Numerical
+build-type:          Simple
+cabal-version:       >=1.10
+extra-source-files:  tests/Perf.hs, tests/LinearFit.hs
+
+source-repository head
+  type: git
+  location: git://github.com/acowley/RANSAC.git
+
+library
+  exposed-modules:     Numeric.Ransac
+  build-depends:       base >= 4.6 && < 5, 
+                       vector >= 0.10, 
+                       random >= 1.0
+  hs-source-dirs:      src
+  default-language:    Haskell2010
+  ghc-options:         -Wall
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/src/Numeric/Ransac.hs b/src/Numeric/Ransac.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Ransac.hs
@@ -0,0 +1,62 @@
+{-# LANGUAGE BangPatterns #-}
+-- | The RANdom SAmple Consensus (RANSAC) algorithm for estimating the
+-- parameters of a mathematical model from a data set. See
+-- <http://en.wikipedia.org/wiki/RANSAC> for more information.
+module Numeric.Ransac (ransac) where
+import Control.Applicative
+import Control.Monad (replicateM)
+import Data.Vector.Generic ((!))
+import qualified Data.Vector.Generic as V
+import System.Random
+
+--randDistinct :: (Eq a, Random a, Monad m) => Int -> m a -> m [a]
+randDistinct :: Int -> IO Int -> IO [Int]
+randDistinct n gen = go 0 [] []
+  where go !i acc _ | i == n = return acc
+        go !i acc [] = replicateM (n-i) gen >>= go i acc
+        go !i acc (r:rs) = if r `elem` acc 
+                           then go i acc rs 
+                           else go (i+1) (r:acc) rs
+{- SPECIALIZE randDistinct :: Int -> IO Int -> IO [Int] #-}
+
+untilJust :: Monad m => m (Maybe b) -> m b
+untilJust x = go 
+  where go = x >>= maybe go return
+{-# INLINE untilJust #-}
+
+-- | @ransac iter sampleSize agreePct fit residual goodFit pts@ draws
+-- @iter@ samples of size @sampleSize@ from @pts@. The @fit@ function
+-- is used to produce a model from each of these samples. The elements
+-- of @pts@ whose residuals pass the @goodFit@ predicate with respect
+-- to this model are identified as /inliers/, and used to update the
+-- model. The model for which the size of the inliers set is at least
+-- @agreePct@ percent of the entire data set and whose error over all
+-- points is minimal among all sampled models is returned. If no
+-- acceptable model is found (i.e. no model whose inliers were at
+-- least @agreePct@ percent of the entire data set), 'Nothing' is
+-- returned.
+ransac :: (V.Vector v a, V.Vector v d, Num d, Ord d) => 
+          Int -> Int -> Float -> 
+          (v a -> Maybe c) -> (c -> a -> d) -> (d -> Bool) -> 
+          v a -> IO (Maybe (c, v a))
+ransac maxIter sampleSize agree fit residual goodFit pts = genModel >>= go 0
+  where go i r@(model, bestError, inliers)
+          | i == maxIter = if ratioInliers (V.length inliers) < agree
+                           then return Nothing
+                           else return (Just (model, inliers))
+          | otherwise = do r'@(_, err, inliers') <- genModel
+                           if ratioInliers (V.length inliers') >= agree
+                              && err < bestError
+                           then go (i+1) r'
+                           else go (i+1) r
+        sample = V.fromList . map (pts !) <$> 
+                 randDistinct sampleSize ((`rem` n) . abs <$> randomIO)
+        genModel = do model <- untilJust (fit <$> sample)
+                      let !errors = V.map (residual model) pts
+                          !inliers = V.ifilter (const . goodFit . (errors !)) pts
+                          Just model' = fit inliers
+                          err = V.sum $ V.map (residual model') pts
+                      return (model', err, inliers)
+        n = V.length pts
+        ratioInliers n' = fromIntegral n' / fromIntegral n
+{-# INLINE ransac #-}
diff --git a/tests/LinearFit.hs b/tests/LinearFit.hs
new file mode 100644
--- /dev/null
+++ b/tests/LinearFit.hs
@@ -0,0 +1,97 @@
+-- | Example use of the RANSAC algorithm to fit a line to some
+-- points. We start with points generated by a process defined by the
+-- equation of a line in 2D. These points are affected by normally
+-- distributed noise, and our data set is further corrupted by a "red
+-- herring" cluster of points that we would like to ignore. We use
+-- RANSAC to cut through the noise and fit a line to the point data
+-- set.
+-- 
+-- The important feature of RANSAC as applied here is that it manages
+-- to ignore the spurious (red herring) cluster centered at (0,8).
+--
+-- The Chart package is used to visualize the data and estimated
+-- model.
+module Main where
+import Control.Applicative
+import Control.Lens (view)
+import Data.Accessor ((^=))
+import Data.Colour (opaque)
+import Data.Colour.Names
+import qualified Data.Foldable as F
+import Data.Random.Normal (normalsIO')
+import Data.Vector.Storable (Vector)
+import qualified Data.Vector.Storable as V
+import Graphics.Rendering.Chart hiding (Vector,Point)
+import Linear
+import Numeric.Ransac
+
+type Point = V2 Float
+
+-- | Fit a 2D line to a collection of 'Point's.
+fitLine :: Vector Point -> Maybe (V2 Float)
+fitLine pts = (!* b) <$> inv22 a
+  where sx = V.sum $ V.map (view _x) pts
+        a = V2 (V2 (V.sum (V.map ((^2).view _x) pts)) sx)
+               (V2 sx (fromIntegral (V.length pts)))
+        b = V2 (V.sum (V.map F.product pts))
+               (V.sum (V.map (view _y) pts))
+
+-- | Compute the error of a 'Point' with respect to a hypothesized
+-- linear model.
+ptError :: V2 Float -> Point -> Float
+ptError (V2 m b) (V2 x y) = sq $ y - (m*x+b)
+  where sq x = x * x
+
+-- | Produce a plot of all the points we have to work with. A green
+-- dashed line indicates the ground truth linear model, the solid
+-- purple line shows the RANSAC model, and the points that are inliers
+-- for that model are circled in yellow.
+main = do noise <- v2Cast . V.fromList . take (n*2) <$> normalsIO' (0,0.3)
+          herring <- V.zipWith V2 
+                     <$> (V.fromList . take 200 <$> normalsIO' (0,0.2))
+                     <*> (V.fromList . take 200 <$> normalsIO' (8,0.6))
+          let pts' = V.zipWith (+) noise pts
+          let pts'' = pts' V.++ herring
+          res <- ransac 100 2 0.5 fitLine ptError (< 2) pts''
+          case res of
+            Nothing -> putStrLn "No model found"
+            Just (model,inliers) -> 
+              do putStrLn $ "Model "++show model++" with "++
+                            show (V.length inliers)++" inliers"
+                 let pp = PlotPoints "data" 
+                                     (filledCircles 2 (opaque blue))
+                                     (map (toTup . dub) (V.toList pts''))
+                     ppi = PlotPoints "inliers"
+                                      (hollowCircles 3 2 (opaque yellow))
+                                      (map (toTup . dub) (V.toList inliers))
+                     lp = PlotLines "truth"
+                                    (dashedLine 3 [10,10] (opaque green))
+                                    [[ toTup $ dub (mkPt 0) 
+                                     , toTup $ dub (mkPt (n-1)) ]]
+                                    []
+                     lp' = PlotLines "model"
+                                     (solidLine 5 (opaque purple))
+                                     [[ toTup $ dub (mkPt' model 0)
+                                      , toTup $ dub (mkPt' model (n-1)) ]]
+                                     []
+                     layout = layout1_title ^="2D Linear Fit"
+                            $ layout1_background ^= solidFillStyle (opaque white)
+                            $ layout1_plots ^= [ Left (toPlot pp)
+                                               , Left (toPlot ppi)
+                                               , Left (toPlot lp)
+                                               , Left (toPlot lp') ]
+                            $ setLayout1Foreground (opaque black)
+                            $ defaultLayout1
+                 renderableToPDFFile (toRenderable layout) 600 600 "foo.pdf"
+  where n = 1000
+        pts = V.generate n mkPt
+        mkPt :: Int -> V2 Float
+        mkPt i = let x = fromIntegral i / 500
+                 in V2 x (5*x + 2)
+        v2Cast :: Vector Float -> Vector Point
+        v2Cast = V.unsafeCast
+        toTup (V2 x y) = (x,y)
+        dub :: V2 Float -> V2 Double
+        dub = fmap realToFrac
+        mkPt' (V2 m b) i = let x = fromIntegral i / 500
+                           in V2 x (x * m + b)
diff --git a/tests/Perf.hs b/tests/Perf.hs
new file mode 100644
--- /dev/null
+++ b/tests/Perf.hs
@@ -0,0 +1,44 @@
+{-# LANGUAGE BangPatterns #-}
+module Main (main) where
+import Control.Applicative
+import Control.Lens (view)
+import Criterion.Main
+import Data.Accessor ((^=))
+import Data.Colour (opaque)
+import Data.Colour.Names
+import qualified Data.Foldable as F
+import Data.Random.Normal (normalsIO')
+import Data.Vector.Storable (Vector)
+import qualified Data.Vector.Storable as V
+import Linear
+import Numeric.Ransac
+
+type Point = V2 Float
+
+fitLine :: Vector Point -> Maybe (V2 Float)
+fitLine pts = (!* b) <$> inv22 a
+  where sx = V.sum $ V.map (view _x) pts
+        a = V2 (V2 (V.sum (V.map ((^2).view _x) pts)) sx)
+               (V2 sx (fromIntegral $ V.length pts))
+        b = V2 (V.sum (V.map F.product pts))
+               (V.sum (V.map (view _y) pts))
+
+ptError :: V2 Float -> Point -> Float
+ptError (V2 m b) (V2 x y) = sq $ y - (m*x+b)
+  where sq x = x * x
+
+main = do noise <- v2Cast . V.fromList . take (n*2) <$> normalsIO' (0,0.3)
+          let !pts' = V.zipWith (+) noise pts
+              ran = ransac 100 2 0.6 fitLine ptError (< 1) pts'
+          putStr $ "Sanity "
+          ran >>= putStrLn . show . fmap fst
+          defaultMain [ bench "linear fit" $ fmap (quadrance . fst) <$> ran ]
+  where n = 10000
+        !pts = V.generate n mkPt
+        mkPt :: Int -> V2 Float
+        mkPt i = let x = fromIntegral i / 500
+                 in V2 x (3*x + 2)
+        v2Cast :: Vector Float -> Vector Point
+        v2Cast = V.unsafeCast
+
+