diff --git a/lib/Numeric/Sampling.hs b/lib/Numeric/Sampling.hs
--- a/lib/Numeric/Sampling.hs
+++ b/lib/Numeric/Sampling.hs
@@ -13,6 +13,10 @@
   , resample
   , resampleIO
 
+    -- * Unequal probability, without replacement
+  , psample
+  , psampleIO
+
     -- * Unequal probability, with replacement
   , presample
   , presampleIO
@@ -21,17 +25,19 @@
   , module System.Random.MWC
   ) where
 
-import qualified Control.Foldl               as F
-import           Control.Monad.Primitive     (PrimMonad, PrimState)
-import qualified Data.Foldable               as Foldable
+import qualified Control.Foldl as F
+import Control.Monad.Primitive (PrimMonad, PrimState)
+import qualified Data.Foldable as Foldable
 #if __GLASGOW_HASKELL__ < 710
 import Data.Foldable (Foldable)
 #endif
-import           Data.Function               (on)
-import           Data.List                   (sortBy)
-import qualified Data.Vector                 as V (toList)
-import           Numeric.Sampling.Internal
-import           System.Random.MWC
+import Data.Function (on)
+import Data.List (sortBy)
+import Data.Monoid
+import qualified Data.Sequence as S
+import qualified Data.Vector as V (toList)
+import Numeric.Sampling.Internal
+import System.Random.MWC
 
 -- | (/O(n)/) Sample uniformly, without replacement.
 --
@@ -42,7 +48,9 @@
   => Int -> f a -> Gen (PrimState m) -> m (Maybe [a])
 sample n xs gen
   | n < 0     = return Nothing
-  | otherwise = fmap (fmap V.toList) (F.foldM (randomN n gen) xs)
+  | otherwise = do
+      collected <- F.foldM (randomN n gen) xs
+      return $ fmap V.toList collected
 {-# INLINABLE sample #-}
 
 -- | (/O(n)/) 'sample' specialized to IO.
@@ -62,12 +70,52 @@
 {-# INLINABLE resample #-}
 
 -- | (/O(n log n)/) 'resample' specialized to IO.
-resampleIO :: (Foldable f) => Int -> f a -> IO [a]
+resampleIO :: Foldable f => Int -> f a -> IO [a]
 resampleIO n xs = do
   gen <- createSystemRandom
   resample n xs gen
 {-# INLINABLE resampleIO #-}
 
+-- | (/O(n log n)/) Unequal probability sampling.
+--
+--   Returns Nothing if the desired sample size is larger than the collection
+--   being sampled from.
+psample
+  :: (PrimMonad m, Foldable f)
+  => Int -> f (Double, a) -> Gen (PrimState m) -> m (Maybe [a])
+psample n weighted gen = do
+    let sorted = sortProbs weighted
+    computeSample n sorted gen
+  where
+    computeSample
+      :: PrimMonad m
+      => Int -> [(Double, a)] -> Gen (PrimState m) -> m (Maybe [a])
+    computeSample size xs g = go 1 [] size (S.fromList xs) where
+      go !mass !acc j vs
+        | j <  0    = return Nothing
+        | j <= 0    = return (Just acc)
+        | otherwise = do
+            z <- fmap (* mass) (uniform g)
+
+            let cumulative = S.drop 1 $ S.scanl (\s (pr, _) -> s + pr) 0 vs
+                midx       = S.findIndexL (>= z) cumulative
+
+                idx = case midx of
+                  Nothing -> error "psample: no index found"
+                  Just x  -> x
+
+                (p, val) = S.index vs idx
+                (l, r)   = S.splitAt idx vs
+                deleted = l <> S.drop 1 r
+
+            go (mass - p) (val:acc) (pred j) deleted
+{-# INLINABLE psample #-}
+
+-- | (/O(n log n)/) 'psample' specialized to IO.
+psampleIO :: Foldable f => Int -> f (Double, a) -> IO (Maybe [a])
+psampleIO n weighted = withSystemRandom . asGenIO $ psample n weighted
+{-# INLINABLE psampleIO #-}
+
 -- | (/O(n log n)/) Unequal probability resampling.
 presample
   :: (PrimMonad m, Foldable f)
@@ -90,15 +138,14 @@
             case F.fold (F.find ((>= z) . fst)) xs of
               Just (_, val) -> go (val:acc) (pred s)
               Nothing       -> return acc
-
-    sortProbs :: (Foldable f, Ord a) => f (a, b) -> [(a, b)]
-    sortProbs = sortBy (compare `on` fst) . Foldable.toList
 {-# INLINABLE presample #-}
 
 -- | (/O(n log n)/) 'presample' specialized to IO.
 presampleIO :: (Foldable f) => Int -> f (Double, a) -> IO [a]
-presampleIO n weighted = do
-  gen <- createSystemRandom
-  presample n weighted gen
+presampleIO n weighted = withSystemRandom . asGenIO $ presample n weighted
 {-# INLINABLE presampleIO #-}
+
+sortProbs :: (Foldable f, Ord a) => f (a, b) -> [(a, b)]
+sortProbs = sortBy (flip compare `on` fst) . Foldable.toList
+{-# INLINABLE sortProbs #-}
 
diff --git a/sampling.cabal b/sampling.cabal
--- a/sampling.cabal
+++ b/sampling.cabal
@@ -1,5 +1,5 @@
 name:                sampling
-version:             0.2.0
+version:             0.3.0
 synopsis:            Sample values from collections.
 homepage:            https://github.com/jtobin/sampling
 license:             MIT
@@ -18,6 +18,9 @@
   * 'sample', for sampling without replacement
   .
   * 'resample', for sampling with replacement (i.e., a bootstrap)
+  .
+  Each variation can be prefixed with 'p' to sample from a container of values
+  weighted by probability.
 
 Source-repository head
   Type:     git
@@ -33,18 +36,18 @@
   exposed-modules:
       Numeric.Sampling
   build-depends:
-      base        < 5
+      base        > 4 && < 6
+    , containers  >= 0.5 && < 1
     , foldl       >= 1.1 && < 2
     , mwc-random  >= 0.13 && < 0.14
-    , primitive
+    , primitive   >= 0.6  && < 1
     , vector      >= 0.11 && < 0.12
 
-executable sampling-test
-  hs-source-dirs: src
-  Main-is:        Main.hs
-  default-language:  Haskell2010
-  ghc-options:
-    -Wall -O2
+Test-suite resample
+  type:                exitcode-stdio-1.0
+  hs-source-dirs:      test
+  Main-is:             Main.hs
+  default-language:    Haskell2010
   build-depends:
       base
     , sampling
diff --git a/src/Main.hs b/src/Main.hs
deleted file mode 100644
--- a/src/Main.hs
+++ /dev/null
@@ -1,10 +0,0 @@
-{-# OPTIONS_GHC -fno-warn-type-defaults #-}
-
-module Main where
-
-import Numeric.Sampling
-
-main :: IO ()
-main = do
-  foo <- resampleIO 100 ([1..100000] :: [Int])
-  print foo
diff --git a/test/Main.hs b/test/Main.hs
new file mode 100644
--- /dev/null
+++ b/test/Main.hs
@@ -0,0 +1,10 @@
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+module Main where
+
+import Numeric.Sampling
+
+main :: IO ()
+main = do
+  foo <- resampleIO 100 ([1..100000] :: [Int])
+  print foo