diff --git a/dataframe-operations.cabal b/dataframe-operations.cabal
--- a/dataframe-operations.cabal
+++ b/dataframe-operations.cabal
@@ -1,6 +1,6 @@
-cabal-version:      2.4
+cabal-version:      3.4
 name:               dataframe-operations
-version:            1.1.1.1
+version:            2.0.0.0
 synopsis:           Column operations, expression DSL, and statistics for the dataframe ecosystem.
 description:
     Untyped column operations (select, filter, sort, join, groupBy,
@@ -27,12 +27,22 @@
         -Wunused-local-binds
         -Wunused-packages
 
+library internal
+    import:             warnings
+    visibility:         public
+    exposed-modules:    DataFrame.Internal.Statistics
+    build-depends:      base >= 4 && < 5,
+                        dataframe-core:internal >= 2.0 && < 2.1,
+                        vector >= 0.13 && < 0.15,
+                        vector-algorithms >= 0.9 && < 0.11
+    hs-source-dirs:     src-internal
+    default-language:   Haskell2010
+
 library
     import:             warnings
     exposed-modules:
                         DataFrame.Functions
                         DataFrame.Monad
-                        DataFrame.Internal.Statistics
                         DataFrame.Operations.AggregateScatter
                         DataFrame.Operations.Aggregation
                         DataFrame.Operations.Core
@@ -48,19 +58,26 @@
                         DataFrame.Operations.Typing
                         DataFrame.Typed.Access
                         DataFrame.Typed.Aggregate
+                        DataFrame.Typed.Apply
                         DataFrame.Typed.Expr
+                        DataFrame.Typed.Expr.Extra
                         DataFrame.Typed.Join
                         DataFrame.Typed.Operations
+                        DataFrame.Typed.Sampling
+                        DataFrame.Typed.Statistics
     build-depends:      base >= 4 && < 5,
-                        bytestring >= 0.11 && < 0.13,
-                        containers >= 0.6.7 && < 0.9,
-                        dataframe-core >= 1.1 && < 1.2,
-                        dataframe-parsing ^>= 1.0.2,
+                        bytestring >= 0.11 && < 0.14,
+                        containers >= 0.6.7 && < 0.10,
+                        dataframe-core >= 2.0 && < 2.1,
+                        dataframe-core:internal >= 2.0 && < 2.1,
+                        dataframe-parsing >= 2.0 && < 2.1,
+                        dataframe-parsing:internal >= 2.0 && < 2.1,
+                        dataframe-operations:internal,
                         random >= 1.2 && < 2,
                         regex-tdfa >= 1.3.0 && < 2,
                         text >= 2.1 && < 3,
                         time >= 1.12 && < 2,
-                        vector ^>= 0.13,
-                        vector-algorithms ^>= 0.9
+                        vector >= 0.13 && < 0.15,
+                        vector-algorithms >= 0.9 && < 0.11
     hs-source-dirs:     src
     default-language:   Haskell2010
diff --git a/src-internal/DataFrame/Internal/Statistics.hs b/src-internal/DataFrame/Internal/Statistics.hs
new file mode 100644
--- /dev/null
+++ b/src-internal/DataFrame/Internal/Statistics.hs
@@ -0,0 +1,285 @@
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE OverloadedStrings #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+
+module DataFrame.Internal.Statistics where
+
+import qualified Data.Vector as V
+import qualified Data.Vector.Algorithms.Intro as VA
+import qualified Data.Vector.Mutable as VM
+import qualified Data.Vector.Unboxed as VU
+import qualified Data.Vector.Unboxed.Mutable as VUM
+
+import Control.Exception (throw)
+import Control.Monad.ST (runST)
+import DataFrame.Errors (DataFrameException (..))
+
+mean' :: (Real a, VU.Unbox a) => VU.Vector a -> Double
+mean' samp
+    | VU.null samp = throw $ EmptyDataSetException "mean"
+    | otherwise = rtf (VU.sum samp) / fromIntegral (VU.length samp)
+{-# INLINE [0] mean' #-}
+
+meanDouble' :: VU.Vector Double -> Double
+meanDouble' samp
+    | VU.null samp = throw $ EmptyDataSetException "mean"
+    | otherwise = VU.sum samp / fromIntegral (VU.length samp)
+{-# INLINE meanDouble' #-}
+
+meanInt' :: VU.Vector Int -> Double
+meanInt' samp
+    | VU.null samp = throw $ EmptyDataSetException "mean"
+    | otherwise = fromIntegral (VU.sum samp) / fromIntegral (VU.length samp)
+{-# INLINE meanInt' #-}
+
+{-# RULES
+"mean'/Double" [1] forall (xs :: VU.Vector Double).
+    mean' xs =
+        meanDouble' xs
+"mean'/Int" [1] forall (xs :: VU.Vector Int).
+    mean' xs =
+        meanInt' xs
+    #-}
+
+median' :: (Real a, VU.Unbox a) => VU.Vector a -> Double
+median' samp
+    | VU.null samp = throw $ EmptyDataSetException "median"
+    | otherwise = runST $ do
+        mutableSamp <- VU.thaw samp
+        VA.sort mutableSamp
+        let len = VU.length samp
+            middleIndex = len `div` 2
+        middleElement <- VUM.read mutableSamp middleIndex
+        if odd len
+            then pure (rtf middleElement)
+            else do
+                prev <- VUM.read mutableSamp (middleIndex - 1)
+                pure (rtf (middleElement + prev) / 2)
+{-# INLINE median' #-}
+
+-- accumulator: count, mean, m2
+data VarAcc
+    = VarAcc {-# UNPACK #-} !Int {-# UNPACK #-} !Double {-# UNPACK #-} !Double
+    deriving (Show)
+
+varianceStep :: VarAcc -> Double -> VarAcc
+varianceStep (VarAcc !n !meanVal !m2) !x =
+    let !n' = n + 1
+        !delta = x - meanVal
+        !meanVal' = meanVal + delta / fromIntegral n'
+        !m2' = m2 + delta * (x - meanVal')
+     in VarAcc n' meanVal' m2'
+{-# INLINE varianceStep #-}
+
+computeVariance :: VarAcc -> Double
+computeVariance (VarAcc !n _ !m2)
+    | n < 2 = 0 -- or error "variance of <2 samples"
+    | otherwise = m2 / fromIntegral (n - 1)
+{-# INLINE computeVariance #-}
+
+variance' :: (Real a, VU.Unbox a) => VU.Vector a -> Double
+variance' = computeVariance . VU.foldl' varianceStep (VarAcc 0 0 0) . VU.map rtf
+{-# INLINE variance' #-}
+
+varianceDouble' :: VU.Vector Double -> Double
+varianceDouble' = computeVariance . VU.foldl' varianceStep (VarAcc 0 0 0)
+{-# INLINE varianceDouble' #-}
+
+-- accumulator: count, mean, m2, m3
+data SkewAcc = SkewAcc !Int !Double !Double !Double deriving (Show)
+
+skewnessStep :: (VU.Unbox a, Num a, Real a) => SkewAcc -> a -> SkewAcc
+skewnessStep (SkewAcc !n !meanVal !m2 !m3) !x' =
+    let !n' = n + 1
+        x = rtf x'
+        !k = fromIntegral n'
+        !delta = x - meanVal
+        !meanVal' = meanVal + delta / k
+        !m2' = m2 + (delta ^ (2 :: Int) * (k - 1)) / k
+        !m3' =
+            m3
+                + (delta ^ (3 :: Int) * (k - 1) * (k - 2)) / k ^ (2 :: Int)
+                - (3 * delta * m2) / k
+     in SkewAcc n' meanVal' m2' m3'
+{-# INLINE skewnessStep #-}
+
+computeSkewness :: SkewAcc -> Double
+computeSkewness (SkewAcc n _ m2 m3)
+    | n < 3 = 0 -- or error "skewness of <3 samples"
+    | otherwise = (sqrt (fromIntegral n - 1) * m3) / sqrt (m2 ^ (3 :: Int))
+{-# INLINE computeSkewness #-}
+
+skewness' :: (VU.Unbox a, Real a, Num a) => VU.Vector a -> Double
+skewness' = computeSkewness . VU.foldl' skewnessStep (SkewAcc 0 0 0 0)
+{-# INLINE skewness' #-}
+
+data CorrelationStats
+    = CorrelationStats
+        {-# UNPACK #-} !Double
+        {-# UNPACK #-} !Double
+        {-# UNPACK #-} !Double
+        {-# UNPACK #-} !Double
+        {-# UNPACK #-} !Double
+
+correlation' :: VU.Vector Double -> VU.Vector Double -> Maybe Double
+correlation' xs ys
+    | n < 2 = Nothing
+    | VU.length xs /= VU.length ys = Nothing
+    | otherwise =
+        let nf = fromIntegral n
+            initial = CorrelationStats 0 0 0 0 0
+            (CorrelationStats sumX sumY sumXX sumYY sumXY) = VU.ifoldl' step initial xs
+
+            !num = nf * sumXY - sumX * sumY
+            !den = sqrt ((nf * sumXX - sumX * sumX) * (nf * sumYY - sumY * sumY))
+         in Just (num / den)
+  where
+    n = VU.length xs
+    step (CorrelationStats sx sy sxx syy sxy) i x =
+        let !y = VU.unsafeIndex ys i
+         in CorrelationStats (sx + x) (sy + y) (sxx + x * x) (syy + y * y) (sxy + x * y)
+{-# INLINE correlation' #-}
+
+quantiles' ::
+    (VU.Unbox a, Num a, Real a) =>
+    VU.Vector Int -> Int -> VU.Vector a -> VU.Vector Double
+quantiles' qs q samp
+    | VU.null samp = throw $ EmptyDataSetException "quantiles"
+    | q < 2 = throw $ WrongQuantileNumberException q
+    | VU.any (\i -> i < 0 || i > q) qs = throw $ WrongQuantileIndexException qs q
+    | otherwise = runST $ do
+        let !n = VU.length samp
+        mutableSamp <- VU.thaw samp
+        VA.sort mutableSamp
+        VU.mapM
+            ( \i -> do
+                let !p = fromIntegral i / fromIntegral q
+                    !position = p * fromIntegral (n - 1) :: Double
+                    !index = floor position :: Int
+                    !f = position - fromIntegral index
+                x <- fmap rtf (VUM.read mutableSamp index)
+                if f == 0
+                    then return x
+                    else do
+                        y <- fmap rtf (VUM.read mutableSamp (index + 1))
+                        return $ (1 - f) * x + f * y
+            )
+            qs
+{-# INLINE quantiles' #-}
+
+percentile' :: (VU.Unbox a, Num a, Real a) => Int -> VU.Vector a -> Double
+percentile' n = VU.head . quantiles' (VU.fromList [n]) 100
+
+quantilesOrd' ::
+    (Ord a, Eq a) =>
+    VU.Vector Int -> Int -> V.Vector a -> V.Vector a
+quantilesOrd' qs q samp
+    | V.null samp = throw $ EmptyDataSetException "quantiles"
+    | q < 2 = throw $ WrongQuantileNumberException q
+    | VU.any (\i -> i < 0 || i > q) qs = throw $ WrongQuantileIndexException qs q
+    | otherwise = runST $ do
+        let !n = V.length samp
+        mutableSamp <- V.thaw samp
+        VA.sort mutableSamp
+        V.mapM
+            ( \i -> do
+                let !p = fromIntegral i / fromIntegral q :: Double
+                    !position = p * fromIntegral (n - 1)
+                    !index = floor position :: Int
+                -- This is not exact for Ord instances.
+                -- Figure out how to make it so.
+                VM.read mutableSamp index
+            )
+            (V.convert qs)
+
+percentileOrd' :: (Ord a, Eq a) => Int -> V.Vector a -> a
+percentileOrd' n = V.head . quantilesOrd' (VU.fromList [n]) 100
+
+interQuartileRange' :: (VU.Unbox a, Num a, Real a) => VU.Vector a -> Double
+interQuartileRange' samp =
+    let quartiles = quantiles' (VU.fromList [1, 3]) 4 samp
+     in quartiles VU.! 1 - quartiles VU.! 0
+{-# INLINE interQuartileRange' #-}
+
+meanSquaredError :: VU.Vector Double -> VU.Vector Double -> Maybe Double
+meanSquaredError target prediction =
+    let
+        squareDiff = VU.ifoldl' (\sq i e -> (e - target VU.! i) ^ (2 :: Int) + sq) 0 prediction
+     in
+        Just $ squareDiff / fromIntegral (max (VU.length target) (VU.length prediction))
+{-# INLINE meanSquaredError #-}
+
+mutualInformationBinned ::
+    Int -> VU.Vector Double -> VU.Vector Double -> Maybe Double
+mutualInformationBinned k xs ys
+    | VU.length xs /= VU.length ys = Nothing
+    | VU.null xs = Nothing
+    | k < 2 = Nothing
+    | rx <= 0 || ry <= 0 = Just 0
+    | otherwise =
+        let bx = VU.map (binIndex xmin xmax k) xs
+            by = VU.map (binIndex ymin ymax k) ys
+            n = fromIntegral (VU.length xs) :: Double
+            mx = bincount k bx
+            my = bincount k by
+            mxy = jointBincount k bx by
+         in Just $
+                sum
+                    [ let !cxy = fromIntegral c
+                          !pxy = cxy / n
+                          !px = fromIntegral (mx VU.! i) / n
+                          !py = fromIntegral (my VU.! j) / n
+                       in if c == 0 then 0 else pxy * logBase 2 (pxy / (px * py))
+                    | i <- [0 .. k - 1]
+                    , j <- [0 .. k - 1]
+                    , let !c = mxy VU.! (i * k + j)
+                    ]
+  where
+    (xmin, xmax) = (VU.minimum xs, VU.maximum xs)
+    (ymin, ymax) = (VU.minimum ys, VU.maximum ys)
+    rx = xmax - xmin
+    ry = ymax - ymin
+
+binIndex :: Double -> Double -> Int -> Double -> Int
+binIndex lo hi k x
+    | hi == lo = 0
+    | otherwise =
+        let !t = (x - lo) / (hi - lo)
+            !ix = floor (fromIntegral k * t) :: Int
+         in max 0 (min (k - 1) ix)
+{-# INLINE binIndex #-}
+
+bincount :: Int -> VU.Vector Int -> VU.Vector Int
+bincount k bs = VU.create $ do
+    mv <- VU.thaw (VU.replicate k 0)
+    VU.forM_ bs $ \b -> do
+        let i
+                | b < 0 = 0
+                | b >= k = k - 1
+                | otherwise = b
+        x <- VUM.read mv i
+        VUM.write mv i (x + 1)
+    pure mv
+{-# INLINE bincount #-}
+
+jointBincount :: Int -> VU.Vector Int -> VU.Vector Int -> VU.Vector Int
+jointBincount k bx by = VU.create $ do
+    mv <- VU.thaw (VU.replicate (k * k) 0)
+    VU.forM_ (VU.zip bx by) $ \(i, j) -> do
+        let ii = clamp i 0 (k - 1)
+            jj = clamp j 0 (k - 1)
+            ix = ii * k + jj
+        x <- VUM.read mv ix
+        VUM.write mv ix (x + 1)
+    pure mv
+  where
+    clamp z a b = max a (min b z)
+{-# INLINE jointBincount #-}
+
+rtf :: (Real a) => a -> Double
+rtf = realToFrac
+{-# NOINLINE [1] rtf #-}
+
+{-# RULES
+"rtf/Double" [2] forall (x :: Double). rtf x = x
+    #-}
diff --git a/src/DataFrame/Internal/Statistics.hs b/src/DataFrame/Internal/Statistics.hs
deleted file mode 100644
--- a/src/DataFrame/Internal/Statistics.hs
+++ /dev/null
@@ -1,285 +0,0 @@
-{-# LANGUAGE BangPatterns #-}
-{-# LANGUAGE OverloadedStrings #-}
-{-# LANGUAGE ScopedTypeVariables #-}
-
-module DataFrame.Internal.Statistics where
-
-import qualified Data.Vector as V
-import qualified Data.Vector.Algorithms.Intro as VA
-import qualified Data.Vector.Mutable as VM
-import qualified Data.Vector.Unboxed as VU
-import qualified Data.Vector.Unboxed.Mutable as VUM
-
-import Control.Exception (throw)
-import Control.Monad.ST (runST)
-import DataFrame.Errors (DataFrameException (..))
-
-mean' :: (Real a, VU.Unbox a) => VU.Vector a -> Double
-mean' samp
-    | VU.null samp = throw $ EmptyDataSetException "mean"
-    | otherwise = rtf (VU.sum samp) / fromIntegral (VU.length samp)
-{-# INLINE [0] mean' #-}
-
-meanDouble' :: VU.Vector Double -> Double
-meanDouble' samp
-    | VU.null samp = throw $ EmptyDataSetException "mean"
-    | otherwise = VU.sum samp / fromIntegral (VU.length samp)
-{-# INLINE meanDouble' #-}
-
-meanInt' :: VU.Vector Int -> Double
-meanInt' samp
-    | VU.null samp = throw $ EmptyDataSetException "mean"
-    | otherwise = fromIntegral (VU.sum samp) / fromIntegral (VU.length samp)
-{-# INLINE meanInt' #-}
-
-{-# RULES
-"mean'/Double" [1] forall (xs :: VU.Vector Double).
-    mean' xs =
-        meanDouble' xs
-"mean'/Int" [1] forall (xs :: VU.Vector Int).
-    mean' xs =
-        meanInt' xs
-    #-}
-
-median' :: (Real a, VU.Unbox a) => VU.Vector a -> Double
-median' samp
-    | VU.null samp = throw $ EmptyDataSetException "median"
-    | otherwise = runST $ do
-        mutableSamp <- VU.thaw samp
-        VA.sort mutableSamp
-        let len = VU.length samp
-            middleIndex = len `div` 2
-        middleElement <- VUM.read mutableSamp middleIndex
-        if odd len
-            then pure (rtf middleElement)
-            else do
-                prev <- VUM.read mutableSamp (middleIndex - 1)
-                pure (rtf (middleElement + prev) / 2)
-{-# INLINE median' #-}
-
--- accumulator: count, mean, m2
-data VarAcc
-    = VarAcc {-# UNPACK #-} !Int {-# UNPACK #-} !Double {-# UNPACK #-} !Double
-    deriving (Show)
-
-varianceStep :: VarAcc -> Double -> VarAcc
-varianceStep (VarAcc !n !meanVal !m2) !x =
-    let !n' = n + 1
-        !delta = x - meanVal
-        !meanVal' = meanVal + delta / fromIntegral n'
-        !m2' = m2 + delta * (x - meanVal')
-     in VarAcc n' meanVal' m2'
-{-# INLINE varianceStep #-}
-
-computeVariance :: VarAcc -> Double
-computeVariance (VarAcc !n _ !m2)
-    | n < 2 = 0 -- or error "variance of <2 samples"
-    | otherwise = m2 / fromIntegral (n - 1)
-{-# INLINE computeVariance #-}
-
-variance' :: (Real a, VU.Unbox a) => VU.Vector a -> Double
-variance' = computeVariance . VU.foldl' varianceStep (VarAcc 0 0 0) . VU.map rtf
-{-# INLINE variance' #-}
-
-varianceDouble' :: VU.Vector Double -> Double
-varianceDouble' = computeVariance . VU.foldl' varianceStep (VarAcc 0 0 0)
-{-# INLINE varianceDouble' #-}
-
--- accumulator: count, mean, m2, m3
-data SkewAcc = SkewAcc !Int !Double !Double !Double deriving (Show)
-
-skewnessStep :: (VU.Unbox a, Num a, Real a) => SkewAcc -> a -> SkewAcc
-skewnessStep (SkewAcc !n !meanVal !m2 !m3) !x' =
-    let !n' = n + 1
-        x = rtf x'
-        !k = fromIntegral n'
-        !delta = x - meanVal
-        !meanVal' = meanVal + delta / k
-        !m2' = m2 + (delta ^ (2 :: Int) * (k - 1)) / k
-        !m3' =
-            m3
-                + (delta ^ (3 :: Int) * (k - 1) * (k - 2)) / k ^ (2 :: Int)
-                - (3 * delta * m2) / k
-     in SkewAcc n' meanVal' m2' m3'
-{-# INLINE skewnessStep #-}
-
-computeSkewness :: SkewAcc -> Double
-computeSkewness (SkewAcc n _ m2 m3)
-    | n < 3 = 0 -- or error "skewness of <3 samples"
-    | otherwise = (sqrt (fromIntegral n - 1) * m3) / sqrt (m2 ^ (3 :: Int))
-{-# INLINE computeSkewness #-}
-
-skewness' :: (VU.Unbox a, Real a, Num a) => VU.Vector a -> Double
-skewness' = computeSkewness . VU.foldl' skewnessStep (SkewAcc 0 0 0 0)
-{-# INLINE skewness' #-}
-
-data CorrelationStats
-    = CorrelationStats
-        {-# UNPACK #-} !Double
-        {-# UNPACK #-} !Double
-        {-# UNPACK #-} !Double
-        {-# UNPACK #-} !Double
-        {-# UNPACK #-} !Double
-
-correlation' :: VU.Vector Double -> VU.Vector Double -> Maybe Double
-correlation' xs ys
-    | n < 2 = Nothing
-    | VU.length xs /= VU.length ys = Nothing
-    | otherwise =
-        let nf = fromIntegral n
-            initial = CorrelationStats 0 0 0 0 0
-            (CorrelationStats sumX sumY sumXX sumYY sumXY) = VU.ifoldl' step initial xs
-
-            !num = nf * sumXY - sumX * sumY
-            !den = sqrt ((nf * sumXX - sumX * sumX) * (nf * sumYY - sumY * sumY))
-         in Just (num / den)
-  where
-    n = VU.length xs
-    step (CorrelationStats sx sy sxx syy sxy) i x =
-        let !y = VU.unsafeIndex ys i
-         in CorrelationStats (sx + x) (sy + y) (sxx + x * x) (syy + y * y) (sxy + x * y)
-{-# INLINE correlation' #-}
-
-quantiles' ::
-    (VU.Unbox a, Num a, Real a) =>
-    VU.Vector Int -> Int -> VU.Vector a -> VU.Vector Double
-quantiles' qs q samp
-    | VU.null samp = throw $ EmptyDataSetException "quantiles"
-    | q < 2 = throw $ WrongQuantileNumberException q
-    | VU.any (\i -> i < 0 || i > q) qs = throw $ WrongQuantileIndexException qs q
-    | otherwise = runST $ do
-        let !n = VU.length samp
-        mutableSamp <- VU.thaw samp
-        VA.sort mutableSamp
-        VU.mapM
-            ( \i -> do
-                let !p = fromIntegral i / fromIntegral q
-                    !position = p * fromIntegral (n - 1) :: Double
-                    !index = floor position :: Int
-                    !f = position - fromIntegral index
-                x <- fmap rtf (VUM.read mutableSamp index)
-                if f == 0
-                    then return x
-                    else do
-                        y <- fmap rtf (VUM.read mutableSamp (index + 1))
-                        return $ (1 - f) * x + f * y
-            )
-            qs
-{-# INLINE quantiles' #-}
-
-percentile' :: (VU.Unbox a, Num a, Real a) => Int -> VU.Vector a -> Double
-percentile' n = VU.head . quantiles' (VU.fromList [n]) 100
-
-quantilesOrd' ::
-    (Ord a, Eq a) =>
-    VU.Vector Int -> Int -> V.Vector a -> V.Vector a
-quantilesOrd' qs q samp
-    | V.null samp = throw $ EmptyDataSetException "quantiles"
-    | q < 2 = throw $ WrongQuantileNumberException q
-    | VU.any (\i -> i < 0 || i > q) qs = throw $ WrongQuantileIndexException qs q
-    | otherwise = runST $ do
-        let !n = V.length samp
-        mutableSamp <- V.thaw samp
-        VA.sort mutableSamp
-        V.mapM
-            ( \i -> do
-                let !p = fromIntegral i / fromIntegral q :: Double
-                    !position = p * fromIntegral (n - 1)
-                    !index = floor position :: Int
-                -- This is not exact for Ord instances.
-                -- Figure out how to make it so.
-                VM.read mutableSamp index
-            )
-            (V.convert qs)
-
-percentileOrd' :: (Ord a, Eq a) => Int -> V.Vector a -> a
-percentileOrd' n = V.head . quantilesOrd' (VU.fromList [n]) 100
-
-interQuartileRange' :: (VU.Unbox a, Num a, Real a) => VU.Vector a -> Double
-interQuartileRange' samp =
-    let quartiles = quantiles' (VU.fromList [1, 3]) 4 samp
-     in quartiles VU.! 1 - quartiles VU.! 0
-{-# INLINE interQuartileRange' #-}
-
-meanSquaredError :: VU.Vector Double -> VU.Vector Double -> Maybe Double
-meanSquaredError target prediction =
-    let
-        squareDiff = VU.ifoldl' (\sq i e -> (e - target VU.! i) ^ (2 :: Int) + sq) 0 prediction
-     in
-        Just $ squareDiff / fromIntegral (max (VU.length target) (VU.length prediction))
-{-# INLINE meanSquaredError #-}
-
-mutualInformationBinned ::
-    Int -> VU.Vector Double -> VU.Vector Double -> Maybe Double
-mutualInformationBinned k xs ys
-    | VU.length xs /= VU.length ys = Nothing
-    | VU.null xs = Nothing
-    | k < 2 = Nothing
-    | rx <= 0 || ry <= 0 = Just 0
-    | otherwise =
-        let bx = VU.map (binIndex xmin xmax k) xs
-            by = VU.map (binIndex ymin ymax k) ys
-            n = fromIntegral (VU.length xs) :: Double
-            mx = bincount k bx
-            my = bincount k by
-            mxy = jointBincount k bx by
-         in Just $
-                sum
-                    [ let !cxy = fromIntegral c
-                          !pxy = cxy / n
-                          !px = fromIntegral (mx VU.! i) / n
-                          !py = fromIntegral (my VU.! j) / n
-                       in if c == 0 then 0 else pxy * logBase 2 (pxy / (px * py))
-                    | i <- [0 .. k - 1]
-                    , j <- [0 .. k - 1]
-                    , let !c = mxy VU.! (i * k + j)
-                    ]
-  where
-    (xmin, xmax) = (VU.minimum xs, VU.maximum xs)
-    (ymin, ymax) = (VU.minimum ys, VU.maximum ys)
-    rx = xmax - xmin
-    ry = ymax - ymin
-
-binIndex :: Double -> Double -> Int -> Double -> Int
-binIndex lo hi k x
-    | hi == lo = 0
-    | otherwise =
-        let !t = (x - lo) / (hi - lo)
-            !ix = floor (fromIntegral k * t) :: Int
-         in max 0 (min (k - 1) ix)
-{-# INLINE binIndex #-}
-
-bincount :: Int -> VU.Vector Int -> VU.Vector Int
-bincount k bs = VU.create $ do
-    mv <- VU.thaw (VU.replicate k 0)
-    VU.forM_ bs $ \b -> do
-        let i
-                | b < 0 = 0
-                | b >= k = k - 1
-                | otherwise = b
-        x <- VUM.read mv i
-        VUM.write mv i (x + 1)
-    pure mv
-{-# INLINE bincount #-}
-
-jointBincount :: Int -> VU.Vector Int -> VU.Vector Int -> VU.Vector Int
-jointBincount k bx by = VU.create $ do
-    mv <- VU.thaw (VU.replicate (k * k) 0)
-    VU.forM_ (VU.zip bx by) $ \(i, j) -> do
-        let ii = clamp i 0 (k - 1)
-            jj = clamp j 0 (k - 1)
-            ix = ii * k + jj
-        x <- VUM.read mv ix
-        VUM.write mv ix (x + 1)
-    pure mv
-  where
-    clamp z a b = max a (min b z)
-{-# INLINE jointBincount #-}
-
-rtf :: (Real a) => a -> Double
-rtf = realToFrac
-{-# NOINLINE [1] rtf #-}
-
-{-# RULES
-"rtf/Double" [2] forall (x :: Double). rtf x = x
-    #-}
diff --git a/src/DataFrame/Monad.hs b/src/DataFrame/Monad.hs
--- a/src/DataFrame/Monad.hs
+++ b/src/DataFrame/Monad.hs
@@ -1,16 +1,46 @@
 {-# LANGUAGE ExplicitNamespaces #-}
 {-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE GADTs #-}
 {-# LANGUAGE InstanceSigs #-}
 {-# LANGUAGE RankNTypes #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TupleSections #-}
 
-module DataFrame.Monad where
+module DataFrame.Monad (
+    -- * The frame-state monad
+    FrameM,
+    runFrameM,
+    evalFrameM,
+    execFrameM,
+    modifyM,
+    inspectM,
 
+    -- * Frame verbs
+    deriveM,
+    insertM,
+    renameM,
+    filterWhereM,
+    sampleM,
+    takeM,
+    dropM,
+    columnAsListM,
+    filterJustM,
+    imputeM,
+
+    -- * Column-derivation pipelines
+    Pipeline,
+    letAs,
+    letExpr,
+    pipelineSteps,
+    toFrameM,
+    runPipeline,
+) where
+
+import Control.Monad (void)
 import DataFrame.Internal.Column (Columnable)
 import DataFrame.Internal.DataFrame (DataFrame)
-import DataFrame.Internal.Expression (Expr (..))
+import DataFrame.Internal.Expression (Expr (..), UExpr (..), prettyPrint)
 import DataFrame.Internal.Nullable (BaseType)
 import qualified DataFrame.Operations.Core as D
 import qualified DataFrame.Operations.Subset as D
@@ -105,3 +135,55 @@
 
 execFrameM :: DataFrame -> FrameM a -> DataFrame
 execFrameM df m = snd (runFrameM df m)
+
+newtype Pipeline a = Pipeline {unPipeline :: Int -> (Int, [(T.Text, UExpr)], a)}
+
+instance Functor Pipeline where
+    fmap f (Pipeline g) = Pipeline $ \n -> let (n', w, a) = g n in (n', w, f a)
+
+instance Applicative Pipeline where
+    pure x = Pipeline (,[],x)
+    Pipeline gf <*> Pipeline gx = Pipeline $ \n ->
+        let (n1, w1, f) = gf n
+            (n2, w2, x) = gx n1
+         in (n2, w1 ++ w2, f x)
+
+instance Monad Pipeline where
+    Pipeline g >>= f = Pipeline $ \n ->
+        let (n1, w1, a) = g n
+            (n2, w2, b) = unPipeline (f a) n1
+         in (n2, w1 ++ w2, b)
+
+-- | Derive a column under an explicit name, returning a reference to it.
+letAs :: (Columnable a) => T.Text -> Expr a -> Pipeline (Expr a)
+letAs nm e = Pipeline (,[(nm, UExpr e)],Col nm)
+
+-- | Derive a column under a fresh generated name.
+letExpr :: (Columnable a) => Expr a -> Pipeline (Expr a)
+letExpr e = Pipeline $ \n ->
+    let nm = T.pack ('_' : 'v' : show n) in (n + 1, [(nm, UExpr e)], Col nm)
+
+instance Show (Pipeline (Expr a)) where
+    show p =
+        let (_, steps, res) = unPipeline p 0
+         in concatMap
+                (\(nm, UExpr e) -> T.unpack nm ++ " = " ++ prettyPrint e ++ "\n")
+                steps
+                ++ "return "
+                ++ prettyPrint res
+
+-- | Number of column derivations the pipeline performs.
+pipelineSteps :: Pipeline a -> Int
+pipelineSteps p = let (_, w, _) = unPipeline p 0 in length w
+
+-- | Interpret a pipeline into 'FrameM', deriving each logged column in order.
+toFrameM :: Pipeline (Expr a) -> FrameM (Expr a)
+toFrameM p =
+    let (_, steps, res) = unPipeline p 0
+     in mapM_ (\(nm, UExpr e) -> void (deriveM nm e)) steps >> pure res
+
+{- | Run a pipeline over a frame: derive its columns and return the result
+expression (a reference to the final column) and the resulting frame.
+-}
+runPipeline :: DataFrame -> Pipeline (Expr a) -> (Expr a, DataFrame)
+runPipeline df = runFrameM df . toFrameM
diff --git a/src/DataFrame/Operations/Core.hs b/src/DataFrame/Operations/Core.hs
--- a/src/DataFrame/Operations/Core.hs
+++ b/src/DataFrame/Operations/Core.hs
@@ -6,8 +6,47 @@
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeApplications #-}
 
-module DataFrame.Operations.Core where
+module DataFrame.Operations.Core (
+    -- * Dimensions
+    dimensions,
+    nRows,
+    nColumns,
 
+    -- * Construction
+    fromUnnamedColumns,
+    fromRows,
+
+    -- * Insertion
+    insert,
+    insertVector,
+    insertUnboxedVector,
+    insertWithDefault,
+    insertVectorWithDefault,
+
+    -- * Column management
+    cloneColumn,
+    rename,
+    renameMany,
+
+    -- * Inspection
+    describeColumns,
+    valueCounts,
+    valueProportions,
+    showDerivedExpressions,
+
+    -- * Folds & matrix/vector conversions
+    fold,
+    toFloatMatrix,
+    toDoubleMatrix,
+    toIntMatrix,
+    columnAsVector,
+    columnAsList,
+    columnAsIntVector,
+    columnAsDoubleVector,
+    columnAsFloatVector,
+    columnAsUnboxedVector,
+) where
+
 import qualified Data.List as L
 import qualified Data.Map as M
 import qualified Data.Map.Strict as MS
@@ -117,10 +156,8 @@
 -- 'columnNames' is now defined in "DataFrame.Internal.DataFrame" and
 -- re-exported from "DataFrame" at the top level.
 
-{- | Adds a vector to the dataframe. If the vector has less elements than the dataframe and the dataframe is not empty
-the vector is converted to type `Maybe a` filled with `Nothing` to match the size of the dataframe. Similarly,
-if the vector has more elements than what's currently in the dataframe, the other columns in the dataframe are
-change to `Maybe <Type>` and filled with `Nothing`.
+{- | Adds a vector as a named column. Size mismatches are reconciled by making
+the shorter side nullable (`Maybe a`) and padding with `Nothing`.
 
 ==== __Example__
 @
@@ -160,13 +197,9 @@
 insertVector name xs = insertColumn name (fromVector xs)
 {-# INLINE insertVector #-}
 
-{- | Adds a foldable collection to the dataframe. If the collection has less elements than the
-dataframe and the dataframe is not empty
-the collection is converted to type `Maybe a` filled with `Nothing` to match the size of the dataframe. Similarly,
-if the collection has more elements than what's currently in the dataframe, the other columns in the dataframe are
-change to `Maybe <Type>` and filled with `Nothing`.
-
-Be careful not to insert infinite collections with this function as that will crash the program.
+{- | Adds a foldable collection as a named column. Size mismatches are reconciled by
+making the shorter side nullable (`Maybe a`) and padding with `Nothing`.
+Do not pass infinite collections: they are fully forced.
 
 ==== __Example__
 @
@@ -202,7 +235,7 @@
     -- | DataFrame to add column to
     DataFrame ->
     DataFrame
-insert name xs = insertColumn name (fromList (Fold.foldr' (:) [] xs)) -- TODO: Do reflection on container type so we can sometimes avoid the list construction.
+insert name xs = insertColumn name (fromList (Fold.foldr' (:) [] xs))
 {-# INLINE insert #-}
 
 {- | Adds a vector to the dataframe and pads it with a default value if it has less elements than the number of rows.
@@ -295,10 +328,8 @@
         values = xs' ++ replicate (rows - length xs') defaultValue
      in insertColumn name (fromList values) d
 
-{- | /O(n)/ Adds an unboxed vector to the dataframe.
-
-Same as insertVector but takes an unboxed vector. If you insert a vector of numbers through insertVector it will either way be converted
-into an unboxed vector so this function saves that extra work/conversion.
+{- | /O(n)/ Like 'insertVector' but takes an already-unboxed vector,
+skipping the boxed-to-unboxed conversion 'insertVector' would do for numbers.
 -}
 insertUnboxedVector ::
     forall a.
@@ -543,7 +574,6 @@
 
 nulls :: Column -> Int
 nulls (BoxedColumn (Just bm) xs) =
-    -- count null bits in bitmap
     let n = VG.length xs
      in n - VU.foldl' (\acc b -> acc + popCount b) 0 bm
 nulls (BoxedColumn Nothing (xs :: V.Vector a)) = case testEquality (typeRep @a) (typeRep @T.Text) of
@@ -731,14 +761,8 @@
 fold :: (a -> DataFrame -> DataFrame) -> [a] -> DataFrame -> DataFrame
 fold f xs acc = L.foldl' (flip f) acc xs
 
-{- | Returns a dataframe as a two dimensional vector of floats.
-
-Converts all columns in the dataframe to float vectors and transposes them
-into a row-major matrix representation.
-
-This is useful for handing data over into ML systems.
-
-Returns 'Left' with an error if any column cannot be converted to floats.
+{- | The dataframe as a row-major matrix of floats, for handing data to ML systems.
+'Left' if any column cannot be converted to floats.
 -}
 toFloatMatrix ::
     DataFrame -> Either DataFrameException (V.Vector (VU.Vector Float))
@@ -753,14 +777,8 @@
                 (fst (dataframeDimensions df))
                 (\i -> VU.generate (V.length m) (\j -> (m VG.! j) VG.! i))
 
-{- | Returns a dataframe as a two dimensional vector of doubles.
-
-Converts all columns in the dataframe to double vectors and transposes them
-into a row-major matrix representation.
-
-This is useful for handing data over into ML systems.
-
-Returns 'Left' with an error if any column cannot be converted to doubles.
+{- | The dataframe as a row-major matrix of doubles, for handing data to ML systems.
+'Left' if any column cannot be converted to doubles.
 -}
 toDoubleMatrix ::
     DataFrame -> Either DataFrameException (V.Vector (VU.Vector Double))
@@ -775,14 +793,8 @@
                 (fst (dataframeDimensions df))
                 (\i -> VU.generate (V.length m) (\j -> (m VG.! j) VG.! i))
 
-{- | Returns a dataframe as a two dimensional vector of ints.
-
-Converts all columns in the dataframe to int vectors and transposes them
-into a row-major matrix representation.
-
-This is useful for handing data over into ML systems.
-
-Returns 'Left' with an error if any column cannot be converted to ints.
+{- | The dataframe as a row-major matrix of ints, for handing data to ML systems.
+'Left' if any column cannot be converted to ints.
 -}
 toIntMatrix :: DataFrame -> Either DataFrameException (V.Vector (VU.Vector Int))
 toIntMatrix df = case V.foldl'
@@ -833,10 +845,8 @@
             Left e -> throw e
             Right (TColumn col) -> toVector col
 
-{- | Retrieves a column as an unboxed vector of 'Int' values.
-
-Returns 'Left' with a 'DataFrameException' if the column cannot be converted to ints.
-This may occur if the column contains non-numeric data or values outside the 'Int' range.
+{- | A column as an unboxed vector of 'Int' values.
+'Left' if the column cannot be converted to ints (non-numeric or out of range).
 -}
 columnAsIntVector ::
     (Columnable a, Num a) =>
@@ -850,10 +860,8 @@
     Left e -> throw e
     Right (TColumn col) -> toIntVector col
 
-{- | Retrieves a column as an unboxed vector of 'Double' values.
-
-Returns 'Left' with a 'DataFrameException' if the column cannot be converted to doubles.
-This may occur if the column contains non-numeric data.
+{- | A column as an unboxed vector of 'Double' values.
+'Left' if the column cannot be converted to doubles (e.g. non-numeric data).
 -}
 columnAsDoubleVector ::
     (Columnable a, Num a) =>
@@ -870,10 +878,8 @@
     Left e -> throw e
     Right (TColumn col) -> toDoubleVector col
 
-{- | Retrieves a column as an unboxed vector of 'Float' values.
-
-Returns 'Left' with a 'DataFrameException' if the column cannot be converted to floats.
-This may occur if the column contains non-numeric data.
+{- | A column as an unboxed vector of 'Float' values.
+'Left' if the column cannot be converted to floats (e.g. non-numeric data).
 -}
 columnAsFloatVector ::
     (Columnable a, Num a) =>
diff --git a/src/DataFrame/Operations/Join.hs b/src/DataFrame/Operations/Join.hs
--- a/src/DataFrame/Operations/Join.hs
+++ b/src/DataFrame/Operations/Join.hs
@@ -8,8 +8,34 @@
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeApplications #-}
 
-module DataFrame.Operations.Join where
+module DataFrame.Operations.Join (
+    -- * Join types
+    JoinType (..),
 
+    -- * Joins
+    join,
+    innerJoin,
+    leftJoin,
+    rightJoin,
+    fullOuterJoin,
+
+    -- * Low-level join kernels
+
+    {- | Reused by the lazy executor and the parallel-join tests; not a
+    stable public API (candidates for a future @Join.Internal@ split).
+    -}
+    buildHashColumn,
+    buildCompactIndex,
+    hashProbeKernel,
+    hashInnerKernel,
+    hashLeftKernel,
+    innerKernel,
+    parInnerKernel,
+    parLeftKernel,
+    assembleInner,
+    assembleLeft,
+) where
+
 import Control.Applicative ((<|>))
 import Control.Exception (throw)
 import Control.Monad (when)
@@ -54,34 +80,24 @@
 join ::
     JoinType ->
     [T.Text] ->
-    DataFrame -> -- Right hand side
-    DataFrame -> -- Left hand side
+    DataFrame ->
+    DataFrame ->
     DataFrame
 join INNER xs right = innerJoin xs right
 join LEFT xs right = leftJoin xs right
 join RIGHT xs right = rightJoin xs right
 join FULL_OUTER xs right = fullOuterJoin xs right
 
-{- | Row-count threshold for the build side.
-When the build side exceeds this, sort-merge join is used instead of the
-single-threaded hash join. The parallel chunked-probe hash join
-('parInnerKernel' \/ 'parLeftKernel') is preferred above this size when more
-than one capability is available (see 'shouldParallelizeJoin'): partitioning
-the probe over cores beats sort-merge for the large build sides measured here
-(e.g. the 2M-row build in the join pipeline). Sort-merge remains the
-single-threaded fallback.
+{- | Build-side row count above which the single-threaded hash join gives way
+to sort-merge (single-threaded) or the parallel chunked-probe hash join when
+multiple capabilities are available (see 'shouldParallelizeJoin').
 -}
 joinStrategyThreshold :: Int
 joinStrategyThreshold = 500_000
 
-{- | A compact index mapping hash values to contiguous slices of
-original row indices. All indices live in a single unboxed vector
-(@ciSortedIndices@, sorted by hash). The lookup table is an open-addressing
-linear-probe hash table held in three parallel unboxed vectors keyed by hash:
-@ciKeys@ holds the hash at each slot, @ciStarts@ the run offset into
-@ciSortedIndices@ (@-1@ marks an empty slot, since real offsets are @>= 0@),
-and @ciLens@ the run length. @ciMask@ is @tableSize - 1@ (table size is a
-power of two), used to map a hash to its home slot.
+{- | Maps hash values to contiguous slices of original row indices. Indices
+live in one hash-sorted vector; an open-addressing linear-probe table (keyed
+by hash) records each run's offset and length. @-1@ starts mark empty slots.
 -}
 data CompactIndex = CompactIndex
     { ciSortedIndices :: {-# UNPACK #-} !(VU.Vector Int)
@@ -121,18 +137,14 @@
         | p > n = p
         | otherwise = go (p * 2)
 
-{- | Build a compact index from a vector of row hashes.
-Sorts @(hash, originalIndex)@ pairs by hash, scans for contiguous runs, then
-inserts each run into an open-addressing linear-probe table sized to keep the
-load factor under ~0.5.
+{- | Build a compact index from a vector of row hashes: sort by hash, scan for
+contiguous runs, insert each run into an open-addressing table. Capacity is
+sized for the worst case (every row distinct) so building never resizes.
 -}
 buildCompactIndex :: VU.Vector Int -> CompactIndex
 buildCompactIndex hashes =
     let n = VU.length hashes
         (sortedHashes, sortedIndices) = parSortByHash n hashes
-        -- Worst case every row is a distinct group; 2*n+1 keeps the table
-        -- sparse even then. Capacity is independent of the actual group count
-        -- so building stays a single pass with no resize.
         !cap = nextPow2Above (2 * n)
         !mask = cap - 1
         (keys, starts, lens) = runST $ do
@@ -169,9 +181,8 @@
     | otherwise = j
 {-# INLINE findGroupEnd #-}
 
-{- | Sort a hash vector, returning sorted hashes and corresponding original indices.
-Sorts an index array using hash values as the comparison key, avoiding the
-intermediate pair vector used by the naive zip-then-sort approach.
+{- | Sort a hash vector, returning sorted hashes and their original indices.
+Sorts an index array keyed by hash, avoiding an intermediate pair vector.
 -}
 sortWithIndices :: VU.Vector Int -> (VU.Vector Int, VU.Vector Int)
 sortWithIndices hashes = runST $ do
@@ -254,6 +265,9 @@
     | D.null right || D.null left = D.empty
     | otherwise = innerJoinNonEmpty cs left right
 
+-- Build the hash index on the smaller side and probe with the larger one;
+-- the probe side is partitioned across cores, so probing the larger maximises
+-- the parallel win.
 innerJoinNonEmpty :: [T.Text] -> DataFrame -> DataFrame -> DataFrame
 innerJoinNonEmpty cs left right =
     let
@@ -267,16 +281,7 @@
         rightHashes = D.computeRowHashes rightKeyIdxs right
 
         buildRows = min leftRows rightRows
-        -- Probe with the larger side, build on the smaller. The probe side
-        -- drives parallelism (it is partitioned across cores), so probing the
-        -- larger side maximises the win.
         probeRows = max leftRows rightRows
-        -- Parallelize the probe either when the build is large enough to spill
-        -- cache (the original sort-merge regime) or when the build is small
-        -- (cache-resident, read-only) but the probe is huge: partitioning a
-        -- ~1e7-row probe over cores wins even against a tiny hot build (the
-        -- medium-factor lever). Both routes use the same bit-identical
-        -- 'parInnerKernel'.
         useParallel =
             shouldParallelizeJoin probeRows buildRows
                 || shouldParallelizeSmallBuildProbe probeRows
@@ -284,10 +289,8 @@
             | buildRows > joinStrategyThreshold && not useParallel =
                 sortMergeInnerKernel leftHashes rightHashes
             | rightRows <= leftRows =
-                -- Build on right (smaller or equal), probe with left
                 innerKernel useParallel leftHashes rightHashes
             | otherwise =
-                -- Build on left (smaller), probe with right, swap result
                 let (!rIxs, !lIxs) = innerKernel useParallel rightHashes leftHashes
                  in (lIxs, rIxs)
      in
@@ -304,10 +307,8 @@
 {-# INLINE innerKernel #-}
 
 {- | Parallel inner-join kernel: build the 'CompactIndex' on @buildHashes@ once,
-then probe @probeHashes@ in parallel. Bit-for-bit identical output to
-'hashInnerKernel' (probe-row order preserved). Runs the IO probe via
-'unsafePerformIO'; the computation is pure (no observable effects, fixed result
-for fixed inputs).
+then probe @probeHashes@ in parallel. Output is bit-for-bit identical to
+'hashInnerKernel' (probe-row order preserved).
 -}
 parInnerKernel ::
     VU.Vector Int -> VU.Vector Int -> (VU.Vector Int, VU.Vector Int)
@@ -325,10 +326,9 @@
         keyIdxs = validatedKeyColIndices "buildHashColumn" csSet df
      in D.computeRowHashes keyIdxs df
 
-{- | Probe one batch of rows against a pre-built 'CompactIndex'.
-Returns @(probeExpandedIxs, buildExpandedIxs)@.
-Unlike 'hashInnerKernel', does not build the index (it is pre-built once)
-and has no cross-product row guard — the caller controls probe batch size.
+{- | Probe one batch of rows against a pre-built 'CompactIndex', returning
+@(probeExpandedIxs, buildExpandedIxs)@. Unlike 'hashInnerKernel' it neither
+builds the index nor guards cross-product size — the caller sizes batches.
 -}
 hashProbeKernel ::
     -- | Built once from the full right\/build side.
@@ -393,16 +393,13 @@
                 <*> VU.unsafeFreeze (VUM.slice 0 total bv)
      in (VU.force pFrozen, VU.force bFrozen)
 
-{- | Hash-based inner join kernel.
-Builds compact index on @buildHashes@ (second arg), probes with
-@probeHashes@ (first arg).
-Returns @(probeExpandedIndices, buildExpandedIndices)@.
-Uses a dynamically growing output buffer to avoid pre-allocating the full
-cross-product size (which can be astronomically large for low-cardinality keys).
+{- | Hash-based inner join kernel: build the index on @buildHashes@, probe with
+@probeHashes@, return @(probeExpandedIndices, buildExpandedIndices)@. Grows its
+output buffer dynamically rather than pre-allocating the full cross product.
 -}
 
-{- | Maximum number of output rows allowed from a join kernel.
-Exceeding this limit indicates a cross-product explosion (e.g. low-cardinality keys).
+{- | Output-row ceiling for a join kernel; exceeding it signals a cross-product
+explosion (e.g. low-cardinality keys).
 -}
 maxJoinOutputRows :: Int
 maxJoinOutputRows = 500_000_000
@@ -472,15 +469,11 @@
             (,)
                 <$> VU.unsafeFreeze (VUM.slice 0 total pv)
                 <*> VU.unsafeFreeze (VUM.slice 0 total bv)
-     in -- VU.force copies the slice into a compact array, releasing the oversized
-        -- backing buffer allocated by the doubling strategy.
-        (VU.force pFrozen, VU.force bFrozen)
+     in (VU.force pFrozen, VU.force bFrozen)
 
-{- | Sort-merge inner join kernel.
-Sorts both sides by hash, walks in lockstep.
-Returns @(leftExpandedIndices, rightExpandedIndices)@.
-Uses a dynamically growing output buffer instead of a two-pass count-then-allocate
-strategy, which OOMs when low-cardinality keys produce large cross products.
+{- | Sort-merge inner join kernel: sort both sides by hash, walk in lockstep,
+return @(leftExpandedIndices, rightExpandedIndices)@. Grows its output buffer
+dynamically so low-cardinality cross products don't OOM.
 -}
 sortMergeInnerKernel ::
     VU.Vector Int -> VU.Vector Int -> (VU.Vector Int, VU.Vector Int)
@@ -565,9 +558,7 @@
             (,)
                 <$> VU.unsafeFreeze (VUM.slice 0 total lv)
                 <*> VU.unsafeFreeze (VUM.slice 0 total rv)
-     in -- VU.force copies the slice into a compact array, releasing the oversized
-        -- backing buffer allocated by the doubling strategy.
-        (VU.force lFrozen, VU.force rFrozen)
+     in (VU.force lFrozen, VU.force rFrozen)
 
 -- | Assemble the result DataFrame for an inner join from expanded index vectors.
 assembleInner ::
@@ -582,7 +573,6 @@
         leftColSet = S.fromList (D.columnNames left)
         rightColNames = D.columnNames right
 
-        -- Pre-expand every column once
         expandedLeftCols = VB.map (D.atIndicesStable leftIxs) (D.columns left)
         expandedRightCols = VB.map (D.atIndicesStable rightIxs) (D.columns right)
 
@@ -594,7 +584,6 @@
             idx <- M.lookup name (D.columnIndices right)
             return (expandedRightCols `VB.unsafeIndex` idx)
 
-        -- Base DataFrame: all left columns, expanded
         baseDf =
             left
                 { columns = expandedLeftCols
@@ -607,15 +596,15 @@
      in D.fold
             ( \name df ->
                 if S.member name csSet
-                    then df -- Key column already present from left side
+                    then df
                     else
                         if S.member name leftColSet
-                            then -- Overlapping non-key column: merge with These
+                            then
                                 insertIfPresent
                                     name
                                     (D.mergeColumns <$> getExpandedLeft name <*> getExpandedRight name)
                                     df
-                            else -- Right-only column
+                            else
                                 insertIfPresent name (getExpandedRight name) df
             )
             rightColNames
@@ -657,6 +646,8 @@
     | D.null left || D.nRows left == 0 = D.empty
     | otherwise = leftJoinNonEmpty callPoint cs left right
 
+-- The right side is always the build side; the left is probed (and drives
+-- parallelism). Unmatched right rows are marked with a @-1@ sentinel.
 leftJoinNonEmpty :: T.Text -> [T.Text] -> DataFrame -> DataFrame -> DataFrame
 leftJoinNonEmpty callPoint cs left right =
     let
@@ -669,10 +660,6 @@
         rightHashes = D.computeRowHashes rightKeyIdxs right
 
         leftRows = fst (D.dimensions left)
-        -- Right is always the build side for left join; left is the probe side
-        -- (and drives parallelism). Parallelize either in the large-build regime
-        -- or when the build is small but the probe (left) is huge: the
-        -- read-only shared index is probed across cores with no synchronization.
         useParallel =
             shouldParallelizeJoin leftRows rightRows
                 || shouldParallelizeSmallBuildProbe leftRows
@@ -684,7 +671,6 @@
             | otherwise =
                 hashLeftKernel leftHashes rightHashes
      in
-        -- rightIxs uses -1 as sentinel for "no match"
         assembleLeft csSet left right leftIxs rightIxs
 
 {- | Parallel left-join kernel: build the 'CompactIndex' on @rightHashes@ once,
@@ -698,11 +684,9 @@
      in unsafePerformIO (parLeftProbe (ciSortedIndices ci) (ciLookup ci) leftHashes)
 {-# NOINLINE parLeftKernel #-}
 
-{- | Hash-based left join kernel.
-Returns @(leftExpandedIndices, rightExpandedIndices)@ where
-right indices use @-1@ as sentinel for unmatched rows.
-Uses a dynamically growing output buffer to avoid pre-allocating the full
-cross-product size (which can be astronomically large for low-cardinality keys).
+{- | Hash-based left join kernel, returning @(leftExpandedIndices,
+rightExpandedIndices)@ with @-1@ marking unmatched right rows. Grows its output
+buffer dynamically rather than pre-allocating the full cross product.
 -}
 hashLeftKernel ::
     VU.Vector Int -> VU.Vector Int -> (VU.Vector Int, VU.Vector Int)
@@ -769,10 +753,9 @@
         <$> VU.unsafeFreeze (VUM.slice 0 total lv)
         <*> VU.unsafeFreeze (VUM.slice 0 total rv)
 
-{- | Sort-merge left join kernel.
-Returns @(leftExpandedIndices, rightExpandedIndices)@ with @-1@ sentinel.
-Uses a dynamically growing output buffer instead of a two-pass count-then-allocate
-strategy, which OOMs when low-cardinality keys produce large cross products.
+{- | Sort-merge left join kernel, returning @(leftExpandedIndices,
+rightExpandedIndices)@ with a @-1@ sentinel. Grows its output buffer
+dynamically so low-cardinality cross products don't OOM.
 -}
 sortMergeLeftKernel ::
     VU.Vector Int -> VU.Vector Int -> (VU.Vector Int, VU.Vector Int)
@@ -967,14 +950,12 @@
         leftHashes = D.computeRowHashes leftKeyIdxs left
         rightHashes = D.computeRowHashes rightKeyIdxs right
 
-        -- Both sides can have nulls in full outer
         (leftIxs, rightIxs)
             | max leftRows rightRows > joinStrategyThreshold =
                 sortMergeFullOuterKernel leftHashes rightHashes
             | otherwise =
                 hashFullOuterKernel leftHashes rightHashes
      in
-        -- Both index vectors use -1 as sentinel
         assembleFullOuter csSet left right leftIxs rightIxs
 
 {- | Hash-based full outer join kernel.
@@ -997,8 +978,6 @@
         !leftCap = VU.length leftStarts
         !rightCap = VU.length rightStarts
 
-    -- Count: matched + left-only + right-only. Iterate over occupied slots
-    -- (ciStarts /= -1) in each table, cross-referencing the other via ciLookup.
     let countLeft !slot !acc
             | slot >= leftCap = acc
             | otherwise =
@@ -1029,7 +1008,6 @@
     rv <- VUM.unsafeNew totalCount
     posRef <- newSTRef (0 :: Int)
 
-    -- Fill matched + left-only (iterate left slots)
     let fillLeft !slot
             | slot >= leftCap = return ()
             | otherwise = do
@@ -1066,7 +1044,6 @@
                         fillLeft (slot + 1)
     fillLeft 0
 
-    -- Fill right-only (iterate right slots not in left)
     let fillRightOnly !slot
             | slot >= rightCap = return ()
             | otherwise = do
@@ -1105,7 +1082,6 @@
         !leftN = VU.length leftHashes
         !rightN = VU.length rightHashes
 
-    -- Pass 1: count
     let countLoop !li !ri !c
             | li >= leftN && ri >= rightN = c
             | li >= leftN = c + (rightN - ri)
@@ -1121,7 +1097,6 @@
             !rh = rightSH `VU.unsafeIndex` ri
         !totalRows = countLoop 0 0 0
 
-    -- Pass 2: fill
     lv <- VUM.unsafeNew totalRows
     rv <- VUM.unsafeNew totalRows
 
@@ -1164,9 +1139,8 @@
     fill 0 0 0
     (,) <$> VU.unsafeFreeze lv <*> VU.unsafeFreeze rv
 
-{- | Assemble the result DataFrame for a full outer join.
-Both index vectors use @-1@ sentinel; all columns gathered via
-'gatherWithSentinel'.  Key columns are coalesced (first non-null wins).
+{- | Assemble the result DataFrame for a full outer join. Both index vectors
+use a @-1@ sentinel; key columns are coalesced (first non-null wins).
 -}
 assembleFullOuter ::
     S.Set T.Text ->
@@ -1201,8 +1175,6 @@
         insertIfPresent _ Nothing df = df
         insertIfPresent name (Just c) df = D.insertColumn name c df
 
-        -- Coalesce two nullable columns: take first non-Nothing per row,
-        -- producing a non-optional column.
         coalesceKeyColumn :: Column -> Column -> Column
         coalesceKeyColumn l r
             | D.isPackedText l || D.isPackedText r =
@@ -1252,10 +1224,9 @@
      in D.fold
             ( \name df ->
                 if S.member name csSet
-                    then -- Key column: coalesce left and right
-                        case (getExpandedLeft name, getExpandedRight name) of
-                            (Just lc, Just rc) -> D.insertColumn name (coalesceKeyColumn lc rc) df
-                            _ -> df
+                    then case (getExpandedLeft name, getExpandedRight name) of
+                        (Just lc, Just rc) -> D.insertColumn name (coalesceKeyColumn lc rc) df
+                        _ -> df
                     else
                         if S.member name leftColSet
                             then
diff --git a/src/DataFrame/Operations/Merge.hs b/src/DataFrame/Operations/Merge.hs
--- a/src/DataFrame/Operations/Merge.hs
+++ b/src/DataFrame/Operations/Merge.hs
@@ -1,7 +1,11 @@
 {-# LANGUAGE InstanceSigs #-}
 {-# OPTIONS_GHC -Wno-orphans #-}
 
-module DataFrame.Operations.Merge where
+module DataFrame.Operations.Merge (
+    -- * Horizontal (side-by-side) merge
+    (|||),
+    -- Also exports the orphan @Semigroup@/@Monoid DataFrame@ instances.
+) where
 
 import qualified Data.List as L
 import qualified Data.Text as T
@@ -38,9 +42,6 @@
                         case optB of
                             Nothing -> case optA of
                                 Nothing ->
-                                    -- N.B. this case should never happen, because we're dealing with columns coming from
-                                    -- union of column names of both dataframes. Nothing + Nothing would mean column
-                                    -- wasn't in either dataframe, which shouldn't happen
                                     D.insertColumn name (D.fromList ([] :: [T.Text])) df
                                 Just a'' ->
                                     D.insertColumn name (D.expandColumn sumRows a'') df
diff --git a/src/DataFrame/Operations/Permutation.hs b/src/DataFrame/Operations/Permutation.hs
--- a/src/DataFrame/Operations/Permutation.hs
+++ b/src/DataFrame/Operations/Permutation.hs
@@ -5,7 +5,15 @@
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeApplications #-}
 
-module DataFrame.Operations.Permutation where
+module DataFrame.Operations.Permutation (
+    -- * Sorting
+    SortOrder (..),
+    sortBy,
+
+    -- * Shuffling
+    shuffle,
+    shuffledIndices,
+) where
 
 import qualified Data.List as L
 import qualified Data.Text as T
diff --git a/src/DataFrame/Operations/Statistics.hs b/src/DataFrame/Operations/Statistics.hs
--- a/src/DataFrame/Operations/Statistics.hs
+++ b/src/DataFrame/Operations/Statistics.hs
@@ -9,8 +9,30 @@
 {-# LANGUAGE UndecidableInstances #-}
 {-# OPTIONS_GHC -Wno-orphans #-}
 
-module DataFrame.Operations.Statistics where
+module DataFrame.Operations.Statistics (
+    -- * Summaries
+    summarize,
+    frequencies,
 
+    -- * Aggregate statistics
+    mean,
+    meanMaybe,
+    median,
+    medianMaybe,
+    percentile,
+    genericPercentile,
+    standardDeviation,
+    skewness,
+    variance,
+    interQuartileRange,
+    correlation,
+    sum,
+
+    -- * Imputation
+    imputeWith,
+    -- Also exports the orphan @ImputeOp (Maybe b)@ instance.
+) where
+
 import qualified Data.List as L
 import qualified Data.Map as M
 import qualified Data.Text as T
@@ -18,7 +40,8 @@
 import qualified Data.Vector.Generic as VG
 import qualified Data.Vector.Unboxed as VU
 
-import Prelude as P
+import Prelude hiding (sum)
+import qualified Prelude as P
 
 import Control.Exception (throw)
 import Data.Function ((&))
@@ -227,7 +250,7 @@
     Nothing ->
         throw $
             ColumnsNotFoundException [name] "_getColumnAsDouble" (M.keys $ columnIndices df)
-    _ -> Nothing -- Return a type mismatch error here.
+    _ -> Nothing
 {-# INLINE _getColumnAsDouble #-}
 
 optionalToDoubleVector :: (Real a) => V.Vector (Maybe a) -> VU.Vector Double
diff --git a/src/DataFrame/Operations/Subset.hs b/src/DataFrame/Operations/Subset.hs
--- a/src/DataFrame/Operations/Subset.hs
+++ b/src/DataFrame/Operations/Subset.hs
@@ -9,8 +9,49 @@
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeApplications #-}
 
-module DataFrame.Operations.Subset where
+module DataFrame.Operations.Subset (
+    -- * Row slicing
+    take,
+    takeLast,
+    drop,
+    dropLast,
+    range,
+    cube,
+    selectRows,
 
+    -- * Filtering
+    filter,
+    filterBy,
+    filterWhere,
+    filterJust,
+    filterNothing,
+    filterAllJust,
+    filterAllNothing,
+
+    -- * Column selection
+    select,
+    selectBy,
+    exclude,
+    SelectionCriteria,
+    byName,
+    byProperty,
+    byNameProperty,
+    byNameRange,
+    byIndexRange,
+
+    -- * Sampling & splitting
+    SplittableGen,
+    sample,
+    randomSplit,
+    kFolds,
+    stratifiedSample,
+    stratifiedSplit,
+
+    -- * Label helpers (exported for satellite packages)
+    columnToTextVec,
+    rowsAtIndices,
+) where
+
 import qualified Data.List as L
 import qualified Data.Map as M
 import qualified Data.Set as S
@@ -52,7 +93,7 @@
 import DataFrame.Operators
 import System.Random
 import Type.Reflection
-import Prelude hiding (filter, take)
+import Prelude hiding (drop, filter, take)
 
 #if MIN_VERSION_random(1,3,0)
 type SplittableGen g = (SplitGen g, RandomGen g)
@@ -138,7 +179,6 @@
     Just c@(PackedText _ _) ->
         filter e condition (insertColumn filterColumnName (materializePacked c) df)
     Just _col@(BoxedColumn bm (column :: V.Vector b)) ->
-        -- Check direct type match first, then try Maybe b match for nullable columns
         case testEquality (typeRep @a) (typeRep @b) of
             Just Refl -> filterByVector filterColumnName column condition df
             Nothing -> case (bm, typeRep @a) of
diff --git a/src/DataFrame/Operations/Transformations.hs b/src/DataFrame/Operations/Transformations.hs
--- a/src/DataFrame/Operations/Transformations.hs
+++ b/src/DataFrame/Operations/Transformations.hs
@@ -10,7 +10,26 @@
 {-# LANGUAGE UndecidableInstances #-}
 {-# LANGUAGE UndecidableSuperClasses #-}
 
-module DataFrame.Operations.Transformations where
+module DataFrame.Operations.Transformations (
+    -- * Apply
+    apply,
+    safeApply,
+    applyMany,
+    applyInt,
+    applyDouble,
+    applyWhere,
+    applyAtIndex,
+
+    -- * Derive
+    derive,
+    deriveWithExpr,
+    deriveMany,
+
+    -- * Impute
+    ImputeOp (..),
+    impute,
+    imputeCore,
+) where
 
 import qualified Data.List as L
 import qualified Data.Map as M
diff --git a/src/DataFrame/Operations/Typing.hs b/src/DataFrame/Operations/Typing.hs
--- a/src/DataFrame/Operations/Typing.hs
+++ b/src/DataFrame/Operations/Typing.hs
@@ -44,18 +44,15 @@
     unsafeGetColumn,
  )
 import DataFrame.Internal.Parsing
-import DataFrame.Internal.Schema
 import DataFrame.Operations.Core ()
 import DataFrame.Operations.Inference
+import DataFrame.Schema
 import Text.Read
 import Type.Reflection
 
-{- | How parse failures are surfaced in the resulting column.
-
-* 'NoSafeRead' — strict parsing: failures throw (via 'read').
-* 'MaybeRead' — failures become 'Nothing'; columns are wrapped as @Maybe a@.
-* 'EitherRead' — failures become @Left rawText@; columns are wrapped as
-  @Either Text a@, preserving the original input so callers can inspect it.
+{- | How parse failures are surfaced: 'NoSafeRead' throws, 'MaybeRead' yields
+@Nothing@ (column wrapped @Maybe a@), 'EitherRead' yields @Left rawText@
+(column wrapped @Either Text a@, preserving the original input).
 -}
 data SafeReadMode
     = NoSafeRead
@@ -70,16 +67,10 @@
     , sampleSize :: Int
     -- ^ Number of rows to inspect when inferring a column's type (0 = all rows).
     , parseSafe :: SafeReadMode
-    {- ^ Default 'SafeReadMode' applied to every column that does not have an
-    entry in 'parseSafeOverrides'. 'NoSafeRead' only treats empty strings as
-    missing; 'MaybeRead' additionally treats 'missingValues' and nullish
-    strings as @Nothing@; 'EitherRead' wraps the resulting column as
-    @Either Text a@ with the raw input preserved on failure.
-    -}
+    -- ^ Default 'SafeReadMode' for columns without a 'parseSafeOverrides' entry.
     , parseSafeOverrides :: [(T.Text, SafeReadMode)]
-    {- ^ Per-column overrides. When a column name is present here, its value
-    takes precedence over 'parseSafe'. Typical use: strict IDs
-    (@NoSafeRead@) alongside lenient fields (@MaybeRead@/@EitherRead@).
+    {- ^ Per-column overrides taking precedence over 'parseSafe' — e.g. strict
+    IDs (@NoSafeRead@) alongside lenient fields (@MaybeRead@/@EitherRead@).
     -}
     , parseDateFormat :: DateFormat
     -- ^ Date format string as accepted by "Data.Time.Format" (e.g. @\"%Y-%m-%d\"@).
@@ -108,7 +99,6 @@
 parseDefaults :: ParseOptions -> DataFrame -> DataFrame
 parseDefaults opts df = df{columns = V.imap forCol (columns df)}
   where
-    -- Index -> column name: reverse the columnIndices map once.
     nameAt =
         let inverted = M.fromList [(i, n) | (n, i) <- M.toList (columnIndices df)]
          in \i -> M.findWithDefault "" i inverted
@@ -144,10 +134,6 @@
     let isNull = case parseSafe opts of
             NoSafeRead -> T.null
             _ -> isNullishOrMissing (missingValues opts)
-        -- `examples` is small (≤ sampleSize, default 100), so the
-        -- Maybe-wrap allocation here is ignorable.  The full-column
-        -- equivalent (`asMaybeText = V.map ... cols`) has been removed:
-        -- handlers now walk `cols` directly with `isNull`.
         examples = V.map (classify isNull) (V.take (sampleSize opts) cols)
         dfmt = parseDateFormat opts
         assumption = makeParsingAssumption dfmt examples
@@ -165,10 +151,9 @@
   where
     classify p t = if p t then Nothing else Just t
 
-{- | For 'EitherRead' mode: take the chosen parsing assumption and produce an
-@Either Text a@ column. Successful parses become @Right@; any row that fails
-to parse as the chosen type (including null/missing cells) becomes @Left@
-carrying the raw input text verbatim.
+{- | For 'EitherRead' mode: parse under the chosen assumption into an
+@Either Text a@ column. Successful parses become @Right@; failures (including
+null/missing cells) become @Left@ carrying the raw input verbatim.
 -}
 handleEitherAssumption ::
     DateFormat -> ParsingAssumption -> V.Vector T.Text -> Column
@@ -177,9 +162,6 @@
     IntAssumption -> fromVector (V.map (toEither readInt) raw)
     DoubleAssumption -> fromVector (V.map (toEither readDouble) raw)
     DateAssumption -> fromVector (V.map (toEither (parseTimeOpt dfmt)) raw)
-    -- TextAssumption and NoAssumption degenerate to Either Text Text; treat
-    -- empty strings as Left "" so the convention (Left = missing/failure) stays
-    -- consistent across column types.
     TextAssumption -> fromVector (V.map textToEither raw)
     NoAssumption -> fromVector (V.map textToEither raw)
   where
@@ -234,10 +216,9 @@
         (parseUnboxedColumnWithPred False isNull readBool cols)
         (handleTextAssumption isNull cols)
 
-{- | Int columns: one fused pass with in-place Int -> Double promotion
-('promoteIntColumn'); a cell parsing as neither demotes to Text over
-the retained raw cells. 'readIntStrict' rejects overflow so a huge
-integer promotes to its true 'Double' value instead of wrapping.
+{- | Int columns: one fused pass with in-place Int -> Double promotion; a cell
+parsing as neither demotes the column to Text. 'readIntStrict' rejects overflow
+so a huge integer promotes to 'Double' rather than wrapping.
 -}
 handleIntAssumption :: (T.Text -> Bool) -> V.Vector T.Text -> Column
 handleIntAssumption isNull cols =
@@ -251,9 +232,8 @@
         (parseUnboxedColumnWithPred 0 isNull readDouble cols)
         (handleTextAssumption isNull cols)
 
-{- | Text columns: no parse, just null-marking.  When the whole column
-is non-null we return a plain 'V.Vector T.Text'; otherwise we emit a
-@V.Vector (Maybe T.Text)@ the same shape the old code produced.
+{- | Text columns: no parse, just null-marking. An all-non-null column stays a
+plain @V.Vector T.Text@; otherwise it becomes @V.Vector (Maybe T.Text)@.
 -}
 handleTextAssumption :: (T.Text -> Bool) -> V.Vector T.Text -> Column
 handleTextAssumption isNull cols
@@ -262,19 +242,15 @@
             (V.map (\t -> if isNull t then Nothing else Just t) cols)
     | otherwise = fromVector cols
 
-{- | Date: single parse pass, boxed because 'Day' is not unboxable.
-Bails to 'handleTextAssumption' the moment a non-null cell fails to
-parse as a 'Day'.  Still avoids the outer @V.Vector (Maybe T.Text)@
-allocation — we walk @cols@ directly with @isNull@.
+{- | Date: single boxed parse pass ('Day' is not unboxable). Bails to
+'handleTextAssumption' the moment a non-null cell fails to parse as a 'Day'.
+A column with no nulls keeps type 'Day' rather than 'Maybe Day'.
 -}
 handleDateAssumption ::
     DateFormat -> (T.Text -> Bool) -> V.Vector T.Text -> Column
 handleDateAssumption dateFormat isNull cols =
     case parseBoxedMaybeColumn isNull (parseTimeOpt dateFormat) cols of
         Just (anyNull, vec)
-            -- `vec :: V.Vector (Maybe Day)`.  If no nulls, strip the
-            -- outer 'Maybe' (every cell is guaranteed 'Just') so the
-            -- column type stays 'Day' rather than becoming 'Maybe Day'.
             | anyNull -> fromVector vec
             | otherwise -> fromVector (V.mapMaybe id vec)
         Nothing -> handleTextAssumption isNull cols
@@ -304,11 +280,11 @@
                             Nothing -> return Nothing
     loop 0 False
 
+-- Reached only when the sample was all-null: try each concrete type in turn,
+-- falling back to Text. A column with no nulls keeps type 'Day', not 'Maybe Day'.
 handleNoAssumption ::
     DateFormat -> (T.Text -> Bool) -> V.Vector T.Text -> Column
 handleNoAssumption dateFormat isNull cols
-    -- Only reached when the 100-row sample was all-null.  Try each
-    -- concrete type in turn; fall back to Text otherwise.
     | V.all isNull cols =
         fromVector (V.map (const (Nothing :: Maybe T.Text)) cols)
     | Just (mbm, vec) <- parseUnboxedColumnWithPred False isNull readBool cols =
@@ -319,19 +295,12 @@
         UnboxedColumn mbm vec
     | otherwise = case parseBoxedMaybeColumn isNull (parseTimeOpt dateFormat) cols of
         Just (anyNull, vec)
-            -- `vec :: V.Vector (Maybe Day)`.  If no nulls, strip the
-            -- outer 'Maybe' (every cell is guaranteed 'Just') so the
-            -- column type stays 'Day' rather than becoming 'Maybe Day'.
             | anyNull -> fromVector vec
             | otherwise -> fromVector (V.mapMaybe id vec)
         Nothing -> handleTextAssumption isNull cols
 
-{- | Predicate matching what 'parseSafe == NoSafeRead' previously used:
-only empty strings are treated as missing.
-
-We still expose 'convertNullish' \/ 'convertOnlyEmpty' below because
-other parts of the library reference them, but neither is used by
-'parseFromExamples' any longer.
+{- | True for nullish or explicitly-listed missing strings. ('convertNullish'
+and 'convertOnlyEmpty' below are kept only for external callers.)
 -}
 isNullishOrMissing :: [T.Text] -> T.Text -> Bool
 isNullishOrMissing missing v = isNullish v || v `elem` missing
@@ -344,7 +313,7 @@
 
 unsafeParseTime :: DateFormat -> T.Text -> Day
 unsafeParseTime dateFormat s =
-    parseTimeOrError {- Accept leading/trailing whitespace -}
+    parseTimeOrError
         True
         defaultTimeLocale
         dateFormat
@@ -361,10 +330,8 @@
     hasSameConstructor Nothing Nothing = True
     hasSameConstructor _ _ = False
 
-{- | Re-type columns of a 'DataFrame' according to the supplied schema map.
-The caller provides a @resolveMode@ function that maps a column name to its
-'SafeReadMode' — typically built from a global default plus an overrides map
-via 'effectiveSafeRead'.
+{- | Re-type columns of a 'DataFrame' according to a schema map. @resolveMode@
+maps a column name to its 'SafeReadMode' (typically via 'effectiveSafeRead').
 -}
 parseWithTypes ::
     (T.Text -> SafeReadMode) ->
@@ -379,9 +346,6 @@
             df
             ts
   where
-    -- \| Re-parse a plain (non-Maybe, non-Either) target type according to the
-    -- 'SafeReadMode'. @toStr@ converts column elements to a 'String' ready for
-    -- 'Read'.
     plainType ::
         forall a b.
         (Columnable a, Read a) =>
@@ -392,8 +356,6 @@
         EitherRead -> fromVector (V.map ((readEitherRaw @a) . toStr) col)
 
     asType :: SafeReadMode -> SchemaType -> Column -> Column
-    -- A raw CSV string column may arrive as PackedText; decode to boxed Text
-    -- so the re-parse arms below can read the cells.
     asType mode st c@(PackedText _ _) = asType mode st (materializePacked c)
     asType mode (SType (_ :: P.Proxy a)) c@(BoxedColumn _ (col :: V.Vector b)) = case typeRep @a of
         App t1 _t2 -> case eqTypeRep t1 (typeRep @Maybe) of
diff --git a/src/DataFrame/Typed/Access.hs b/src/DataFrame/Typed/Access.hs
--- a/src/DataFrame/Typed/Access.hs
+++ b/src/DataFrame/Typed/Access.hs
@@ -10,18 +10,35 @@
     -- * Typed column access
     columnAsVector,
     columnAsList,
+
+    -- * Numeric vector extraction
+    columnAsIntVector,
+    columnAsDoubleVector,
+    columnAsFloatVector,
+    columnAsUnboxedVector,
+
+    -- * Matrix extraction
+    toDoubleMatrix,
+    toFloatMatrix,
+    toIntMatrix,
 ) where
 
 import Control.Exception (throw)
 import Data.Proxy (Proxy (..))
 import qualified Data.Text as T
 import qualified Data.Vector as V
+import qualified Data.Vector.Unboxed as VU
 import GHC.TypeLits (KnownSymbol, symbolVal)
 
 import DataFrame.Internal.Column (Columnable)
 import DataFrame.Internal.Expression (Expr (Col))
 import qualified DataFrame.Operations.Core as D
-import DataFrame.Typed.Schema (AssertPresent, SafeLookup)
+import DataFrame.Typed.Schema (
+    AllColumnsReal,
+    AssertPresent,
+    AssertRealColumn,
+    SafeLookup,
+ )
 import DataFrame.Typed.Types (TypedDataFrame (..))
 
 {- | Retrieve a column as a boxed 'Vector', with the type determined by
@@ -53,3 +70,94 @@
     D.columnAsList (Col @a colName) df
   where
     colName = T.pack (symbolVal (Proxy @name))
+
+{- | Retrieve a column coerced to an unboxed 'Int' vector, named by type
+application. The column must exist and be numeric — both are compile-time
+checks via 'SafeLookup', so this is total (no 'Either', no runtime throw).
+-}
+columnAsIntVector ::
+    forall name cols a.
+    ( KnownSymbol name
+    , a ~ SafeLookup name cols
+    , Columnable a
+    , AssertRealColumn "columnAsIntVector" name a
+    , Real a
+    , VU.Unbox a
+    , AssertPresent name cols
+    ) =>
+    TypedDataFrame cols -> VU.Vector Int
+columnAsIntVector (TDF df) = either throw id (D.columnAsIntVector (Col @a colName) df)
+  where
+    colName = T.pack (symbolVal (Proxy @name))
+
+-- | Retrieve a column coerced to an unboxed 'Double' vector. See 'columnAsIntVector'.
+columnAsDoubleVector ::
+    forall name cols a.
+    ( KnownSymbol name
+    , a ~ SafeLookup name cols
+    , Columnable a
+    , AssertRealColumn "columnAsDoubleVector" name a
+    , Real a
+    , VU.Unbox a
+    , AssertPresent name cols
+    ) =>
+    TypedDataFrame cols -> VU.Vector Double
+columnAsDoubleVector (TDF df) =
+    either throw id (D.columnAsDoubleVector (Col @a colName) df)
+  where
+    colName = T.pack (symbolVal (Proxy @name))
+
+-- | Retrieve a column coerced to an unboxed 'Float' vector. See 'columnAsIntVector'.
+columnAsFloatVector ::
+    forall name cols a.
+    ( KnownSymbol name
+    , a ~ SafeLookup name cols
+    , Columnable a
+    , AssertRealColumn "columnAsFloatVector" name a
+    , Real a
+    , VU.Unbox a
+    , AssertPresent name cols
+    ) =>
+    TypedDataFrame cols -> VU.Vector Float
+columnAsFloatVector (TDF df) =
+    either throw id (D.columnAsFloatVector (Col @a colName) df)
+  where
+    colName = T.pack (symbolVal (Proxy @name))
+
+{- | Retrieve a column as an unboxed vector of its own element type. The column
+must exist and be unboxable — both compile-time checks, so this is total.
+-}
+columnAsUnboxedVector ::
+    forall name cols a.
+    ( KnownSymbol name
+    , a ~ SafeLookup name cols
+    , Columnable a
+    , VU.Unbox a
+    , AssertPresent name cols
+    ) =>
+    TypedDataFrame cols -> VU.Vector a
+columnAsUnboxedVector (TDF df) =
+    either throw id (D.columnAsUnboxedVector (Col @a colName) df)
+  where
+    colName = T.pack (symbolVal (Proxy @name))
+
+{- | Convert every column to 'Double' and transpose into a row-major matrix.
+Total: 'AllColumnsReal' proves at compile time that every column is numeric and
+unboxed, so the conversion cannot fail.
+-}
+toDoubleMatrix ::
+    (AllColumnsReal "toDoubleMatrix" cols) =>
+    TypedDataFrame cols -> V.Vector (VU.Vector Double)
+toDoubleMatrix (TDF df) = either throw id (D.toDoubleMatrix df)
+
+-- | Convert every column to 'Float' and transpose into a row-major matrix. See 'toDoubleMatrix'.
+toFloatMatrix ::
+    (AllColumnsReal "toFloatMatrix" cols) =>
+    TypedDataFrame cols -> V.Vector (VU.Vector Float)
+toFloatMatrix (TDF df) = either throw id (D.toFloatMatrix df)
+
+-- | Convert every column to 'Int' and transpose into a row-major matrix. See 'toDoubleMatrix'.
+toIntMatrix ::
+    (AllColumnsReal "toIntMatrix" cols) =>
+    TypedDataFrame cols -> V.Vector (VU.Vector Int)
+toIntMatrix (TDF df) = either throw id (D.toIntMatrix df)
diff --git a/src/DataFrame/Typed/Apply.hs b/src/DataFrame/Typed/Apply.hs
new file mode 100644
--- /dev/null
+++ b/src/DataFrame/Typed/Apply.hs
@@ -0,0 +1,218 @@
+{-# LANGUAGE AllowAmbiguousTypes #-}
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeApplications #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+
+{- | Typed column transformations: the @apply@\/@derive@ family and
+default-valued inserts, plus the horizontal merge @('|||')@. Schema changes are
+tracked at the type level (e.g. 'applyColumn' rewrites a column's element type
+via 'SetColumnType').
+-}
+module DataFrame.Typed.Apply (
+    applyColumn,
+    applyMany,
+    applyWhere,
+    applyAtIndex,
+    safeApply,
+    deriveWithExpr,
+    insertWithDefault,
+    insertVectorWithDefault,
+    insertUnboxedVector,
+    (|||),
+) where
+
+import Data.Proxy (Proxy (..))
+import qualified Data.Text as T
+import qualified Data.Vector as V
+import qualified Data.Vector.Unboxed as VU
+import GHC.TypeLits (KnownSymbol, Symbol, symbolVal)
+
+import DataFrame.Errors (DataFrameException)
+import DataFrame.Internal.Column (Columnable)
+import qualified DataFrame.Operations.Core as D
+import qualified DataFrame.Operations.Merge as D
+import qualified DataFrame.Operations.Transformations as D
+import DataFrame.Typed.Freeze (unsafeFreeze)
+import DataFrame.Typed.Schema (
+    AllKnownSymbol,
+    Append,
+    AssertAbsent,
+    AssertAllColumnsHaveType,
+    AssertDisjoint,
+    AssertPresent,
+    SafeLookup,
+    SetColumnType,
+    Snoc,
+    symbolVals,
+ )
+import DataFrame.Typed.Types (Column, TExpr (..), TypedDataFrame (..))
+
+{- | Map a function over a column, rewriting its element type from @a@ to @b@.
+The schema's entry for @name@ is updated via 'SetColumnType'.
+
+@
+df' = applyColumn \@\"age\" (show :: Int -> String) df
+-- the \"age\" column is now String-typed
+@
+-}
+applyColumn ::
+    forall name a b cols.
+    ( KnownSymbol name
+    , a ~ SafeLookup name cols
+    , Columnable a
+    , Columnable b
+    , AssertPresent name cols
+    ) =>
+    (a -> b) ->
+    TypedDataFrame cols ->
+    TypedDataFrame (SetColumnType name b cols)
+applyColumn f (TDF df) = unsafeFreeze (D.apply f colName df)
+  where
+    colName = T.pack (symbolVal (Proxy @name))
+
+-- | Like 'applyColumn' but returns the error instead of throwing.
+safeApply ::
+    forall name a b cols.
+    ( KnownSymbol name
+    , a ~ SafeLookup name cols
+    , Columnable a
+    , Columnable b
+    , AssertPresent name cols
+    ) =>
+    (a -> b) ->
+    TypedDataFrame cols ->
+    Either DataFrameException (TypedDataFrame (SetColumnType name b cols))
+safeApply f (TDF df) = fmap unsafeFreeze (D.safeApply f colName df)
+  where
+    colName = T.pack (symbolVal (Proxy @name))
+
+{- | Apply a type-preserving function to several columns at once. Every named
+column must already share the element type @a@ (enforced by
+'AssertAllColumnsHaveType').
+-}
+applyMany ::
+    forall (names :: [Symbol]) a cols.
+    (AllKnownSymbol names, Columnable a, AssertAllColumnsHaveType names a cols) =>
+    (a -> a) ->
+    TypedDataFrame cols ->
+    TypedDataFrame cols
+applyMany f (TDF df) = TDF (D.applyMany f (symbolVals @names) df)
+
+{- | Apply a function to a target column only on rows where a condition holds on
+a filter column. Both columns are named by type application; the target keeps
+its type.
+
+@
+applyWhere \@\"flagged\" \@\"score\" id (* 2) df
+@
+-}
+applyWhere ::
+    forall filterName targetName a b cols.
+    ( KnownSymbol filterName
+    , KnownSymbol targetName
+    , a ~ SafeLookup filterName cols
+    , b ~ SafeLookup targetName cols
+    , Columnable a
+    , Columnable b
+    , AssertPresent filterName cols
+    , AssertPresent targetName cols
+    ) =>
+    (a -> Bool) ->
+    (b -> b) ->
+    TypedDataFrame cols ->
+    TypedDataFrame cols
+applyWhere cond f (TDF df) = TDF (D.applyWhere cond filterName f targetName df)
+  where
+    filterName = T.pack (symbolVal (Proxy @filterName))
+    targetName = T.pack (symbolVal (Proxy @targetName))
+
+-- | Apply a type-preserving function to a single row of a column.
+applyAtIndex ::
+    forall name a cols.
+    ( KnownSymbol name
+    , a ~ SafeLookup name cols
+    , Columnable a
+    , AssertPresent name cols
+    ) =>
+    Int ->
+    (a -> a) ->
+    TypedDataFrame cols ->
+    TypedDataFrame cols
+applyAtIndex i f (TDF df) = TDF (D.applyAtIndex i f colName df)
+  where
+    colName = T.pack (symbolVal (Proxy @name))
+
+{- | Derive a new column and also return a typed reference to it. The returned
+expression lives in the extended schema, so it can feed later operations.
+-}
+deriveWithExpr ::
+    forall name a cols.
+    ( KnownSymbol name
+    , Columnable a
+    , AssertAbsent name cols
+    ) =>
+    TExpr cols a ->
+    TypedDataFrame cols ->
+    ( TExpr (Snoc cols (Column name a)) a
+    , TypedDataFrame (Snoc cols (Column name a))
+    )
+deriveWithExpr (TExpr expr) (TDF df) =
+    let (e', df') = D.deriveWithExpr colName expr df
+     in (TExpr e', unsafeFreeze df')
+  where
+    colName = T.pack (symbolVal (Proxy @name))
+
+-- | Insert a column from a 'Foldable', padding missing rows with a default.
+insertWithDefault ::
+    forall name a cols t.
+    ( KnownSymbol name
+    , Columnable a
+    , Foldable t
+    , AssertAbsent name cols
+    ) =>
+    a -> t a -> TypedDataFrame cols -> TypedDataFrame (Column name a ': cols)
+insertWithDefault def xs (TDF df) =
+    unsafeFreeze (D.insertWithDefault def colName xs df)
+  where
+    colName = T.pack (symbolVal (Proxy @name))
+
+-- | Insert a boxed 'V.Vector', padding missing rows with a default.
+insertVectorWithDefault ::
+    forall name a cols.
+    ( KnownSymbol name
+    , Columnable a
+    , AssertAbsent name cols
+    ) =>
+    a -> V.Vector a -> TypedDataFrame cols -> TypedDataFrame (Column name a ': cols)
+insertVectorWithDefault def vec (TDF df) =
+    unsafeFreeze (D.insertVectorWithDefault def colName vec df)
+  where
+    colName = T.pack (symbolVal (Proxy @name))
+
+-- | Insert an unboxed 'VU.Vector' as a new column.
+insertUnboxedVector ::
+    forall name a cols.
+    ( KnownSymbol name
+    , Columnable a
+    , VU.Unbox a
+    , AssertAbsent name cols
+    ) =>
+    VU.Vector a -> TypedDataFrame cols -> TypedDataFrame (Column name a ': cols)
+insertUnboxedVector vec (TDF df) =
+    unsafeFreeze (D.insertUnboxedVector colName vec df)
+  where
+    colName = T.pack (symbolVal (Proxy @name))
+
+{- | Horizontal merge: place two DataFrames side by side. The schemas must be
+disjoint (no shared column names), enforced by 'AssertDisjoint'; the result
+schema is their concatenation.
+-}
+(|||) ::
+    (AssertDisjoint left right) =>
+    TypedDataFrame left ->
+    TypedDataFrame right ->
+    TypedDataFrame (Append left right)
+(TDF a) ||| (TDF b) = unsafeFreeze (a D.||| b)
diff --git a/src/DataFrame/Typed/Expr.hs b/src/DataFrame/Typed/Expr.hs
--- a/src/DataFrame/Typed/Expr.hs
+++ b/src/DataFrame/Typed/Expr.hs
@@ -123,6 +123,9 @@
     -- * Sort helpers
     asc,
     desc,
+
+    -- * Additional expression functions
+    module DataFrame.Typed.Expr.Extra,
 ) where
 
 import Data.Either (fromRight)
@@ -158,6 +161,7 @@
 import DataFrame.Internal.Types (Promote, PromoteDiv)
 
 import qualified Data.Vector.Unboxed as VU
+import DataFrame.Typed.Expr.Extra
 import DataFrame.Typed.Schema (
     AllKnownSymbol,
     AssertAllPresent,
diff --git a/src/DataFrame/Typed/Expr/Extra.hs b/src/DataFrame/Typed/Expr/Extra.hs
new file mode 100644
--- /dev/null
+++ b/src/DataFrame/Typed/Expr/Extra.hs
@@ -0,0 +1,220 @@
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+
+{- | Typed counterparts of the remaining "DataFrame.Functions" expression
+combinators not already provided by "DataFrame.Typed.Expr". Each wraps the
+untyped combinator 1:1, replacing @Expr@ with @'TExpr' cols@ so column
+references stay schema-checked. Re-exported from "DataFrame.Typed.Expr".
+-}
+module DataFrame.Typed.Expr.Extra (
+    div,
+    mod,
+    mode,
+    sumMaybe,
+    meanMaybe,
+    variance,
+    medianMaybe,
+    percentile,
+    stddev,
+    stddevMaybe,
+    zScore,
+    pow,
+    relu,
+    min,
+    max,
+    reduce,
+    toMaybe,
+    fromMaybe,
+    isJust,
+    isNothing,
+    fromJust,
+    whenPresent,
+    whenBothPresent,
+    recode,
+    recodeWithCondition,
+    recodeWithDefault,
+    firstOrNothing,
+    lastOrNothing,
+    splitOn,
+    match,
+    matchAll,
+    parseDate,
+    daysBetween,
+    bind,
+) where
+
+import qualified Data.Text as T
+import Data.Time (Day, ParseTime)
+import qualified Data.Vector.Unboxed as VU
+import Prelude hiding (div, max, min, mod)
+
+import qualified DataFrame.Functions as F
+import DataFrame.Internal.Column (Columnable)
+import DataFrame.Typed.Types (TExpr (..))
+
+-- | Integer division.
+div ::
+    (Integral a, Columnable a) => TExpr cols a -> TExpr cols a -> TExpr cols a
+div (TExpr a) (TExpr b) = TExpr (F.div a b)
+
+-- | Integer modulus.
+mod ::
+    (Integral a, Columnable a) => TExpr cols a -> TExpr cols a -> TExpr cols a
+mod (TExpr a) (TExpr b) = TExpr (F.mod a b)
+
+-- | Most frequent value (aggregation).
+mode :: (Ord a, Columnable a, Eq a) => TExpr cols a -> TExpr cols a
+mode (TExpr e) = TExpr (F.mode e)
+
+-- | Sum of a nullable column, ignoring 'Nothing' (aggregation).
+sumMaybe :: (Columnable a, Num a) => TExpr cols (Maybe a) -> TExpr cols a
+sumMaybe (TExpr e) = TExpr (F.sumMaybe e)
+
+-- | Mean of a nullable column, ignoring 'Nothing' (aggregation).
+meanMaybe :: (Columnable a, Real a) => TExpr cols (Maybe a) -> TExpr cols Double
+meanMaybe (TExpr e) = TExpr (F.meanMaybe e)
+
+-- | Variance (aggregation).
+variance ::
+    (Columnable a, Real a, VU.Unbox a) => TExpr cols a -> TExpr cols Double
+variance (TExpr e) = TExpr (F.variance e)
+
+-- | Median of a nullable column, ignoring 'Nothing' (aggregation).
+medianMaybe ::
+    (Columnable a, Real a) => TExpr cols (Maybe a) -> TExpr cols Double
+medianMaybe (TExpr e) = TExpr (F.medianMaybe e)
+
+-- | The @n@-th percentile (aggregation).
+percentile :: Int -> TExpr cols Double -> TExpr cols Double
+percentile n (TExpr e) = TExpr (F.percentile n e)
+
+-- | Standard deviation (aggregation).
+stddev ::
+    (Columnable a, Real a, VU.Unbox a) => TExpr cols a -> TExpr cols Double
+stddev (TExpr e) = TExpr (F.stddev e)
+
+-- | Standard deviation of a nullable column, ignoring 'Nothing' (aggregation).
+stddevMaybe ::
+    (Columnable a, Real a) => TExpr cols (Maybe a) -> TExpr cols Double
+stddevMaybe (TExpr e) = TExpr (F.stddevMaybe e)
+
+-- | Z-score (value minus group mean, over standard deviation).
+zScore :: TExpr cols Double -> TExpr cols Double
+zScore (TExpr e) = TExpr (F.zScore e)
+
+-- | Raise an expression to an integer power.
+pow :: (Columnable a, Num a) => TExpr cols a -> Int -> TExpr cols a
+pow (TExpr e) i = TExpr (F.pow e i)
+
+-- | Rectified linear unit: @max 0@.
+relu :: (Columnable a, Num a, Ord a) => TExpr cols a -> TExpr cols a
+relu (TExpr e) = TExpr (F.relu e)
+
+-- | Element-wise minimum of two expressions.
+min :: (Columnable a, Ord a) => TExpr cols a -> TExpr cols a -> TExpr cols a
+min (TExpr a) (TExpr b) = TExpr (F.min a b)
+
+-- | Element-wise maximum of two expressions.
+max :: (Columnable a, Ord a) => TExpr cols a -> TExpr cols a -> TExpr cols a
+max (TExpr a) (TExpr b) = TExpr (F.max a b)
+
+-- | Fold a column into a single value with a seed and step function (aggregation).
+reduce ::
+    (Columnable a, Columnable b) =>
+    TExpr cols b -> a -> (a -> b -> a) -> TExpr cols a
+reduce (TExpr e) start f = TExpr (F.reduce e start f)
+
+-- | Wrap each value in 'Just'.
+toMaybe :: (Columnable a) => TExpr cols a -> TExpr cols (Maybe a)
+toMaybe (TExpr e) = TExpr (F.toMaybe e)
+
+-- | Replace 'Nothing' with a default.
+fromMaybe :: (Columnable a) => a -> TExpr cols (Maybe a) -> TExpr cols a
+fromMaybe d (TExpr e) = TExpr (F.fromMaybe d e)
+
+-- | True where the value is 'Just'.
+isJust :: (Columnable a) => TExpr cols (Maybe a) -> TExpr cols Bool
+isJust (TExpr e) = TExpr (F.isJust e)
+
+-- | True where the value is 'Nothing'.
+isNothing :: (Columnable a) => TExpr cols (Maybe a) -> TExpr cols Bool
+isNothing (TExpr e) = TExpr (F.isNothing e)
+
+-- | Unwrap a 'Just', erroring on 'Nothing'.
+fromJust :: (Columnable a) => TExpr cols (Maybe a) -> TExpr cols a
+fromJust (TExpr e) = TExpr (F.fromJust e)
+
+-- | Apply a function only where the value is present.
+whenPresent ::
+    (Columnable a, Columnable b) =>
+    (a -> b) -> TExpr cols (Maybe a) -> TExpr cols (Maybe b)
+whenPresent f (TExpr e) = TExpr (F.whenPresent f e)
+
+-- | Apply a binary function only where both values are present.
+whenBothPresent ::
+    (Columnable a, Columnable b, Columnable c) =>
+    (a -> b -> c) ->
+    TExpr cols (Maybe a) ->
+    TExpr cols (Maybe b) ->
+    TExpr cols (Maybe c)
+whenBothPresent f (TExpr a) (TExpr b) = TExpr (F.whenBothPresent f a b)
+
+-- | Map values through a lookup table, yielding 'Nothing' for misses.
+recode ::
+    (Columnable a, Columnable b, Show a, Show b, Show (a, b)) =>
+    [(a, b)] -> TExpr cols a -> TExpr cols (Maybe b)
+recode mapping (TExpr e) = TExpr (F.recode mapping e)
+
+-- | Pick the first value whose condition holds, else a fallback.
+recodeWithCondition ::
+    (Columnable a, Columnable b) =>
+    TExpr cols b ->
+    [(TExpr cols a -> TExpr cols Bool, b)] ->
+    TExpr cols a ->
+    TExpr cols b
+recodeWithCondition (TExpr fallback) conds (TExpr e) =
+    TExpr (F.recodeWithCondition fallback (map untype conds) e)
+  where
+    untype (p, v) = (unTExpr . p . TExpr, v)
+
+-- | Map values through a lookup table, with a default for misses.
+recodeWithDefault ::
+    (Columnable a, Columnable b, Show (a, b)) =>
+    b -> [(a, b)] -> TExpr cols a -> TExpr cols b
+recodeWithDefault d mapping (TExpr e) = TExpr (F.recodeWithDefault d mapping e)
+
+-- | First element of a list column, or 'Nothing'.
+firstOrNothing :: (Columnable a) => TExpr cols [a] -> TExpr cols (Maybe a)
+firstOrNothing (TExpr e) = TExpr (F.firstOrNothing e)
+
+-- | Last element of a list column, or 'Nothing'.
+lastOrNothing :: (Columnable a) => TExpr cols [a] -> TExpr cols (Maybe a)
+lastOrNothing (TExpr e) = TExpr (F.lastOrNothing e)
+
+-- | Split text on a delimiter.
+splitOn :: T.Text -> TExpr cols T.Text -> TExpr cols [T.Text]
+splitOn delim (TExpr e) = TExpr (F.splitOn delim e)
+
+-- | First regex match, or 'Nothing'.
+match :: T.Text -> TExpr cols T.Text -> TExpr cols (Maybe T.Text)
+match regex (TExpr e) = TExpr (F.match regex e)
+
+-- | All regex matches.
+matchAll :: T.Text -> TExpr cols T.Text -> TExpr cols [T.Text]
+matchAll regex (TExpr e) = TExpr (F.matchAll regex e)
+
+-- | Parse text into a time value with the given format.
+parseDate ::
+    (ParseTime t, Columnable t) =>
+    T.Text -> TExpr cols T.Text -> TExpr cols (Maybe t)
+parseDate format (TExpr e) = TExpr (F.parseDate format e)
+
+-- | Number of days between two dates.
+daysBetween :: TExpr cols Day -> TExpr cols Day -> TExpr cols Int
+daysBetween (TExpr a) (TExpr b) = TExpr (F.daysBetween a b)
+
+-- | Monadic bind over a column of monadic values.
+bind ::
+    (Columnable a, Columnable (m a), Monad m, Columnable b, Columnable (m b)) =>
+    (a -> m b) -> TExpr cols (m a) -> TExpr cols (m b)
+bind f (TExpr e) = TExpr (F.bind f e)
diff --git a/src/DataFrame/Typed/Operations.hs b/src/DataFrame/Typed/Operations.hs
--- a/src/DataFrame/Typed/Operations.hs
+++ b/src/DataFrame/Typed/Operations.hs
@@ -17,6 +17,7 @@
     filterAllJust,
     filterJust,
     filterNothing,
+    filterAllNothing,
     sortBy,
     take,
     takeLast,
@@ -42,6 +43,10 @@
     dropColumn,
     replaceColumn,
 
+    -- * Frequencies
+    valueCounts,
+    valueProportions,
+
     -- * Metadata
     dimensions,
     nRows,
@@ -145,6 +150,12 @@
   where
     colName = T.pack (symbolVal (Proxy @name))
 
+{- | Keep only rows where every nullable column has Nothing.
+Schema is preserved.
+-}
+filterAllNothing :: TypedDataFrame cols -> TypedDataFrame cols
+filterAllNothing (TDF df) = TDF (D.filterAllNothing df)
+
 {- | Sort by the given typed sort orders.
 Sort orders reference columns that are validated against the schema.
 -}
@@ -374,10 +385,6 @@
     TypedDataFrame cols -> TypedDataFrame cols -> TypedDataFrame cols
 symmetricDifference (TDF a) (TDF b) = TDF (DS.symmetricDifference a b)
 
--------------------------------------------------------------------------------
--- Metadata (pass-through)
--------------------------------------------------------------------------------
-
 dimensions :: TypedDataFrame cols -> (Int, Int)
 dimensions (TDF df) = D.dimensions df
 
@@ -390,9 +397,15 @@
 columnNames :: TypedDataFrame cols -> [T.Text]
 columnNames (TDF df) = D.columnNames df
 
--------------------------------------------------------------------------------
--- Internal helpers
--------------------------------------------------------------------------------
+-- | Count occurrences of each distinct value in a column.
+valueCounts ::
+    (Ord a, Columnable a) => TExpr cols a -> TypedDataFrame cols -> [(a, Int)]
+valueCounts (TExpr e) (TDF df) = D.valueCounts e df
+
+-- | Proportion of each distinct value in a column.
+valueProportions ::
+    (Ord a, Columnable a) => TExpr cols a -> TypedDataFrame cols -> [(a, Double)]
+valueProportions (TExpr e) (TDF df) = D.valueProportions e df
 
 -- | Helper class for extracting [(Text, Text)] from type-level pairs.
 class AllKnownPairs (pairs :: [(Symbol, Symbol)]) where
diff --git a/src/DataFrame/Typed/Sampling.hs b/src/DataFrame/Typed/Sampling.hs
new file mode 100644
--- /dev/null
+++ b/src/DataFrame/Typed/Sampling.hs
@@ -0,0 +1,55 @@
+{-# LANGUAGE FlexibleContexts #-}
+
+{- | Typed sampling and splitting. All operations are schema-preserving: they
+change which rows are present, never the columns, so every result reuses the
+input schema @cols@.
+-}
+module DataFrame.Typed.Sampling (
+    randomSplit,
+    kFolds,
+    selectRows,
+    stratifiedSample,
+    stratifiedSplit,
+) where
+
+import System.Random (RandomGen)
+
+import DataFrame.Internal.Column (Columnable)
+import DataFrame.Operations.Subset (SplittableGen)
+import qualified DataFrame.Operations.Subset as D
+import DataFrame.Typed.Types (TExpr (..), TypedDataFrame (..))
+
+-- | Split rows into two DataFrames by a fraction.
+randomSplit ::
+    (RandomGen g) =>
+    g -> Double -> TypedDataFrame cols -> (TypedDataFrame cols, TypedDataFrame cols)
+randomSplit g p (TDF df) = let (a, b) = D.randomSplit g p df in (TDF a, TDF b)
+
+-- | Partition rows into @k@ folds.
+kFolds ::
+    (RandomGen g) => g -> Int -> TypedDataFrame cols -> [TypedDataFrame cols]
+kFolds g k (TDF df) = map TDF (D.kFolds g k df)
+
+{- | Select rows by index.
+| This may fail if the indices are out of bounds;
+| use with caution or use 'filter' to select rows by a predicate instead.
+-}
+selectRows :: [Int] -> TypedDataFrame cols -> TypedDataFrame cols
+selectRows ixs (TDF df) = TDF (D.selectRows ixs df)
+
+-- | Sample a fraction of rows, preserving the distribution of a strata column.
+stratifiedSample ::
+    (SplittableGen g, Columnable a) =>
+    g -> Double -> TExpr cols a -> TypedDataFrame cols -> TypedDataFrame cols
+stratifiedSample g p (TExpr e) (TDF df) = TDF (D.stratifiedSample g p e df)
+
+-- | Split rows by a fraction, preserving the distribution of a strata column.
+stratifiedSplit ::
+    (SplittableGen g, Columnable a) =>
+    g ->
+    Double ->
+    TExpr cols a ->
+    TypedDataFrame cols ->
+    (TypedDataFrame cols, TypedDataFrame cols)
+stratifiedSplit g p (TExpr e) (TDF df) =
+    let (a, b) = D.stratifiedSplit g p e df in (TDF a, TDF b)
diff --git a/src/DataFrame/Typed/Statistics.hs b/src/DataFrame/Typed/Statistics.hs
new file mode 100644
--- /dev/null
+++ b/src/DataFrame/Typed/Statistics.hs
@@ -0,0 +1,180 @@
+{-# LANGUAGE AllowAmbiguousTypes #-}
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeApplications #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+
+{- | Typed statistical reducers over a 'TypedDataFrame'.
+
+These mirror the untyped reducers in "DataFrame.Operations.Statistics", taking
+a schema-checked 'TExpr' instead of a raw @Expr@. The names (@sum@, @mean@,
+@median@, …) deliberately collide with the aggregation-expression combinators
+in "DataFrame.Typed.Expr", so this module is meant to be imported qualified:
+
+@
+import qualified DataFrame.Typed.Statistics as TS
+
+avg = TS.mean (col \@\"salary\") employees
+@
+-}
+module DataFrame.Typed.Statistics (
+    mean,
+    meanMaybe,
+    median,
+    medianMaybe,
+    percentile,
+    genericPercentile,
+    standardDeviation,
+    skewness,
+    variance,
+    interQuartileRange,
+    sum,
+    correlation,
+    frequencies,
+    imputeWith,
+    summarize,
+    describeColumns,
+) where
+
+import Data.Proxy (Proxy (..))
+import qualified Data.Text as T
+import qualified Data.Vector.Unboxed as VU
+import GHC.TypeLits (KnownSymbol, symbolVal)
+import Prelude hiding (sum)
+
+import DataFrame.Internal.Column (Columnable)
+import qualified DataFrame.Internal.DataFrame as D
+import DataFrame.Internal.Nullable (BaseType)
+import qualified DataFrame.Operations.Core as Core
+import qualified DataFrame.Operations.Statistics as Stats
+import DataFrame.Operations.Transformations (ImputeOp)
+import DataFrame.Typed.Schema (AssertPresent, SafeLookup)
+import DataFrame.Typed.Types (TExpr (..), TypedDataFrame (..))
+
+-- | Mean of a column.
+mean ::
+    (Columnable a, Real a, VU.Unbox a) =>
+    TExpr cols a -> TypedDataFrame cols -> Double
+mean (TExpr e) (TDF df) = Stats.mean e df
+
+-- | Mean of a nullable column, ignoring 'Nothing'.
+meanMaybe ::
+    (Columnable a, Real a) =>
+    TExpr cols (Maybe a) -> TypedDataFrame cols -> Double
+meanMaybe (TExpr e) (TDF df) = Stats.meanMaybe e df
+
+-- | Median of a column.
+median ::
+    (Columnable a, Real a, VU.Unbox a) =>
+    TExpr cols a -> TypedDataFrame cols -> Double
+median (TExpr e) (TDF df) = Stats.median e df
+
+-- | Median of a nullable column, ignoring 'Nothing'.
+medianMaybe ::
+    (Columnable a, Real a) =>
+    TExpr cols (Maybe a) -> TypedDataFrame cols -> Double
+medianMaybe (TExpr e) (TDF df) = Stats.medianMaybe e df
+
+-- | The @n@-th percentile of a column.
+percentile ::
+    (Columnable a, Real a, VU.Unbox a) =>
+    Int -> TExpr cols a -> TypedDataFrame cols -> Double
+percentile n (TExpr e) (TDF df) = Stats.percentile n e df
+
+-- | The @n@-th percentile of a column of any 'Ord' type.
+genericPercentile ::
+    (Columnable a, Ord a) =>
+    Int -> TExpr cols a -> TypedDataFrame cols -> a
+genericPercentile n (TExpr e) (TDF df) = Stats.genericPercentile n e df
+
+-- | Standard deviation of a column.
+standardDeviation ::
+    (Columnable a, Real a, VU.Unbox a) =>
+    TExpr cols a -> TypedDataFrame cols -> Double
+standardDeviation (TExpr e) (TDF df) = Stats.standardDeviation e df
+
+-- | Skewness of a column.
+skewness ::
+    (Columnable a, Real a, VU.Unbox a) =>
+    TExpr cols a -> TypedDataFrame cols -> Double
+skewness (TExpr e) (TDF df) = Stats.skewness e df
+
+-- | Variance of a column.
+variance ::
+    (Columnable a, Real a, VU.Unbox a) =>
+    TExpr cols a -> TypedDataFrame cols -> Double
+variance (TExpr e) (TDF df) = Stats.variance e df
+
+-- | Inter-quartile range of a column.
+interQuartileRange ::
+    (Columnable a, Real a, VU.Unbox a) =>
+    TExpr cols a -> TypedDataFrame cols -> Double
+interQuartileRange (TExpr e) (TDF df) = Stats.interQuartileRange e df
+
+-- | Sum of a column.
+sum :: (Columnable a, Num a) => TExpr cols a -> TypedDataFrame cols -> a
+sum (TExpr e) (TDF df) = Stats.sum e df
+
+{- | Pearson's correlation coefficient between two columns, named by type
+application. Both columns must exist in the schema and be numeric — these are
+checked at compile time via 'SafeLookup' on each name.
+
+@
+TS.correlation \@\"height\" \@\"weight\" people
+@
+-}
+correlation ::
+    forall c1 c2 a b cols.
+    ( KnownSymbol c1
+    , KnownSymbol c2
+    , a ~ SafeLookup c1 cols
+    , b ~ SafeLookup c2 cols
+    , Columnable a
+    , Columnable b
+    , Real a
+    , Real b
+    , VU.Unbox a
+    , VU.Unbox b
+    , AssertPresent c1 cols
+    , AssertPresent c2 cols
+    ) =>
+    TypedDataFrame cols -> Maybe Double
+correlation (TDF df) =
+    Stats.correlation
+        (T.pack (symbolVal (Proxy @c1)))
+        (T.pack (symbolVal (Proxy @c2)))
+        df
+
+{- | Frequency table for a column. The result schema is data-dependent
+(one column per distinct value), so an untyped 'D.DataFrame' is returned.
+-}
+frequencies ::
+    (Columnable a, Ord a) => TExpr cols a -> TypedDataFrame cols -> D.DataFrame
+frequencies (TExpr e) (TDF df) = Stats.frequencies e df
+
+{- | Impute missing values in a column using a derived scalar (e.g. the mean).
+Schema-preserving: the imputed column keeps its type-level @Maybe@ even though
+its runtime values are now fully populated.
+-}
+imputeWith ::
+    (ImputeOp a, Columnable (BaseType a)) =>
+    (TExpr cols (BaseType a) -> TExpr cols (BaseType a)) ->
+    TExpr cols a ->
+    TypedDataFrame cols ->
+    TypedDataFrame cols
+imputeWith f (TExpr e) (TDF df) =
+    TDF (Stats.imputeWith (unTExpr . f . TExpr) e df)
+
+{- | Descriptive statistics of the numeric columns. Returns an untyped
+'D.DataFrame' (the result is a fixed set of statistic rows, not the input schema).
+-}
+summarize :: TypedDataFrame cols -> D.DataFrame
+summarize (TDF df) = Stats.summarize df
+
+{- | Per-column summary (non-null\/null counts, unique values, type). Returns an
+untyped 'D.DataFrame'.
+-}
+describeColumns :: TypedDataFrame cols -> D.DataFrame
+describeColumns (TDF df) = Core.describeColumns df