Skip to content
Browse files

preparation for refactoring to support jackknifed trimean, etc

  • Loading branch information...
1 parent bc5b219 commit 600ddae98632f996f8bffd7176d44839f0440d77 @ekmett committed May 2, 2012
View
22 Data/Pass.hs
@@ -1,20 +1,29 @@
module Data.Pass
- ( Thrist(..), thrist
+ ( Thrist(..)
+ , thrist
, Accelerated(..)
, Accelerant(..)
, Prep(..)
- , Calc(..), passes
+ , Calc(..)
+ , passes
, Calculation(..)
, Pass(..)
- , Passable(..), env
+ , Passable(..)
+ , env
, Eval(..)
- , Naive(..)
- , Step(..)
+ , Naive(..), (@@@)
, Call(..)
, Named(..)
, Trans(..)
- , L(..)
+ , L(..), breakdown, Estimator(..)
, Robust(..)
+ , quantile
+ , tercile, t1, t2
+ , quartile, q1, q2, q3
+ , quintile, qu1, qu2, qu3, qu4
+ , percentile
+ , permille
+ , Step(..), midhinge, trimean, iqr, iqm
) where
import Data.Pass.Calc
@@ -23,6 +32,7 @@ import Data.Pass.Accelerated
import Data.Pass.Calculation
import Data.Pass.Call
import Data.Pass.Class
+import Data.Pass.Estimator
import Data.Pass.Eval
import Data.Pass.Eval.Naive
import Data.Pass.Named
View
15 Data/Pass/Accelerant.hs
@@ -0,0 +1,15 @@
+module Data.Pass.Accelerant
+ ( Accelerant(..)
+ ) where
+
+import Data.Pass.Type
+import Data.Pass.L
+import Data.Pass.Robust
+
+-- provide hooks to allow the user to accelerate non-robust L-estimators
+class Accelerant k where
+ meanPass :: Pass k Double Double
+ meanPass = robust LMean
+
+ totalPass :: Pass k Double Double
+ totalPass = robust LTotal
View
35 Data/Pass/Accelerated.hs
@@ -0,0 +1,35 @@
+module Data.Pass.Accelerated
+ ( Accelerated(..)
+ ) where
+
+import Data.Pass.Type
+import Data.Pass.Calc
+import Data.Pass.L
+import Data.Pass.Fun
+import Data.Pass.Thrist
+import Data.Pass.Robust
+import Data.Pass.Accelerant
+
+class Accelerated k where
+ mean :: k Double Double
+ total :: k Double Double
+
+instance Accelerated L where
+ mean = robust LMean
+ total = robust LTotal
+
+instance Accelerated k => Accelerated (Fun k) where
+ mean = Fun mean
+ total = Fun total
+
+instance Accelerated k => Accelerated (Thrist k) where
+ mean = mean :- Nil
+ total = total :- Nil
+
+instance Accelerant k => Accelerated (Calc k) where
+ mean = Stop () `Step` const meanPass
+ total = Stop () `Step` const totalPass
+
+instance Accelerant k => Accelerated (Pass k) where
+ mean = meanPass
+ total = totalPass
View
7 Data/Pass/Calc.hs
@@ -2,11 +2,13 @@
module Data.Pass.Calc
( Calc(..)
, passes
+ , (@@@)
) where
import Control.Category
import Control.Applicative
import Prelude hiding (id,(.))
+import Data.Foldable
import Data.Pass.Call
import Data.Pass.Eval
import Data.Pass.Eval.Naive
@@ -82,3 +84,8 @@ instance Call k => Naive (Calc k) where
instance Call k => Eval (Calc k) where
Stop b @@ _ = b
Step i k @@ xs = k (i @@ xs) @@ xs
+
+infixr 5 @@@
+
+(@@@) :: (Call k, Foldable f) => Calc k a b -> f a -> b
+(@@@) = naive
View
51 Data/Pass/Estimator.hs
@@ -0,0 +1,51 @@
+{-# LANGUAGE DeriveDataTypeable, PatternGuards #-}
+module Data.Pass.Estimator
+ ( Estimator(..)
+ , Estimate(..)
+ , estimateBy
+ ) where
+
+import Data.Ratio
+import Data.Data
+import qualified Data.IntMap as IM
+import Data.IntMap (IntMap)
+import Data.Pass.Util (clamp)
+import Data.Hashable
+
+data Estimator = R1 | R2 | R3 | R4 | R5 | R6 | R7 | R8 | R9 | R10
+ deriving (Eq,Ord,Enum,Bounded,Data,Typeable,Show,Read)
+
+instance Hashable Estimator where
+ hashWithSalt n e = n `hashWithSalt` fromEnum e
+
+data Estimate = Estimate {-# UNPACK #-} !Rational (IntMap Double)
+ deriving Show
+
+continuousEstimator ::
+ (Rational -> (Rational, Rational)) ->
+ (Rational -> Rational -> Rational) ->
+ Rational -> Int -> Estimate
+continuousEstimator bds f p n = Estimate h $
+ if p < lo then IM.singleton 0 1
+ else if p >= hi then IM.singleton (n - 1) 1
+ else case properFraction h of
+ (w,frac) | frac' <- fromRational frac -> IM.fromList [(w - 1, frac'), (w, 1 - frac')]
+ where
+ r = fromIntegral n
+ h = f p r
+ (lo, hi) = bds r
+
+estimateBy :: Estimator -> Rational -> Int -> Estimate
+estimateBy R1 = \p n -> let np = fromIntegral n * p in Estimate (np + 1%2) $ IM.singleton (clamp n (ceiling np - 1)) 1
+estimateBy R2 = \p n -> let np = fromIntegral n * p in Estimate (np + 1%2) $
+ if p == 0 then IM.singleton 0 1
+ else if p == 1 then IM.singleton (n - 1) 1
+ else IM.fromList [(ceiling np - 1, 0.5), (floor np, 0.5)]
+estimateBy R3 = \p n -> let np = fromIntegral n * p in Estimate np $ IM.singleton (clamp n (round np - 1)) 1
+estimateBy R4 = continuousEstimator (\n -> (recip n, 1)) (*)
+estimateBy R5 = continuousEstimator (\n -> let tn = 2 * n in (recip tn, (tn - 1) / tn)) $ \p n -> p*n + 0.5
+estimateBy R6 = continuousEstimator (\n -> (recip (n + 1), n / (n + 1))) $ \p n -> p*(n+1)
+estimateBy R7 = continuousEstimator (\_ -> (0, 1)) $ \p n -> p*(n-1) + 1
+estimateBy R8 = continuousEstimator (\n -> (2/3 / (n + 1/3), (n - 1/3)/(n + 1/3))) $ \p n -> p*(n + 1/3) + 1/3
+estimateBy R9 = continuousEstimator (\n -> (0.625 / (n + 0.25), (n - 0.375)/(n + 0.25))) $ \p n -> p*(n + 0.25) + 0.375
+estimateBy R10 = continuousEstimator (\n -> (1.5 / (n + 2), (n + 0.5)/(n + 2))) $ \p n -> p*(n + 2) - 0.5
View
1 Data/Pass/Eval/Naive.hs
@@ -6,3 +6,4 @@ import Data.Foldable
class Naive k where
naive :: Foldable f => k a b -> f a -> b
+
View
32 Data/Pass/L.hs
@@ -4,19 +4,24 @@
module Data.Pass.L
( L(..)
, callL
+ , breakdown
) where
import Data.Typeable
import Data.Hashable
import Data.Pass.Named
import Data.IntMap (IntMap)
import qualified Data.IntMap as IM
+import Data.Pass.Estimator
+import Data.Pass.Util (clamp)
+
data L a b where
LTotal :: L Double Double
LMean :: L Double Double
LScale :: L Double Double
Median :: L Double Double
+ QuantileBy :: Estimator -> Rational -> L Double Double
Winsorized :: Rational -> L a Double -> L a Double
Trimmed :: Rational -> L a Double -> L a Double
Jackknifed :: L a Double -> L a Double
@@ -37,24 +42,27 @@ instance Named L where
showsFun d (Jackknifed f) = showParen (d > 10) $ showString "Jackknifed " . showsFun 10 f
showsFun d (x :* y) = showParen (d > 7) $ showsPrec 8 x . showString " :* " . showsPrec 7 y
showsFun d (x :+ y) = showParen (d > 6) $ showsPrec 7 x . showString " :+ " . showsPrec 6 y
+ showsFun d (QuantileBy e q) = showParen (d > 10) $ showString "QuantileBy " . showsPrec 10 e . showChar ' ' . showsPrec 10 q
hashFunWithSalt n LTotal = n `hashWithSalt` 0
hashFunWithSalt n LMean = n `hashWithSalt` 1
hashFunWithSalt n LScale = n `hashWithSalt` 2
hashFunWithSalt n Median = n `hashWithSalt` 3
hashFunWithSalt n (Winsorized p f) = n `hashWithSalt` 4 `hashWithSalt` p `hashFunWithSalt` f
hashFunWithSalt n (Trimmed p f) = n `hashWithSalt` 5 `hashWithSalt` p `hashFunWithSalt` f
hashFunWithSalt n (Jackknifed f) = n `hashWithSalt` 6 `hashFunWithSalt` f
- hashFunWithSalt n (x :* y) = n `hashWithSalt` x `hashFunWithSalt` y
- hashFunWithSalt n (x :+ y) = n `hashFunWithSalt` x `hashFunWithSalt` y
+ hashFunWithSalt n (x :* y) = n `hashWithSalt` 7 `hashWithSalt` x `hashFunWithSalt` y
+ hashFunWithSalt n (x :+ y) = n `hashWithSalt` 8 `hashFunWithSalt` x `hashFunWithSalt` y
+ hashFunWithSalt n (QuantileBy e q) = n `hashWithSalt` 9 `hashWithSalt` e `hashWithSalt` q
equalFun LTotal LTotal = True
- equalFun LMean LMean = True
+ equalFun LMean LMean = True
equalFun LScale LScale = True
equalFun Median Median = True
equalFun (Winsorized p f) (Winsorized q g) = p == q && equalFun f g
- equalFun (Trimmed p f) (Trimmed q g) = p == q && equalFun f g
- equalFun (Jackknifed f) (Jackknifed g) = equalFun f g
- equalFun (a :+ b) (c :+ d) = equalFun a c && equalFun b d
- equalFun (a :* b) (c :* d) = a == c && equalFun b d
+ equalFun (Trimmed p f) (Trimmed q g) = p == q && equalFun f g
+ equalFun (Jackknifed f) (Jackknifed g) = equalFun f g
+ equalFun (a :+ b) (c :+ d) = equalFun a c && equalFun b d
+ equalFun (a :* b) (c :* d) = a == c && equalFun b d
+ equalFun (QuantileBy e p) (QuantileBy f q) = e == f && p == q
equalFun _ _ = False
instance Show (L a b) where
@@ -69,6 +77,10 @@ instance Eq (L a b) where
callL :: L a a -> Int -> IntMap a
callL (x :+ y) n = IM.unionWith (+) (callL x n) (callL y n)
callL (s :* y) n = fmap (s*) (callL y n)
+callL (QuantileBy f p) n = case estimateBy f p n of
+ Estimate h qp -> case properFraction h of
+ (w, 0) -> IM.singleton (clamp n (w - 1)) 1
+ _ -> qp
callL Median n = case quotRem n 2 of
(k, 0) -> IM.fromDistinctAscList [(k-1,0.5),(k,0.5)]
(k, _) -> IM.singleton k 1
@@ -93,3 +105,9 @@ callL (Jackknifed g) n = IM.fromAscListWith (+) $ do
(k, v) <- IM.toAscList $ callL g (n - 1)
let k' = fromIntegral k + 1
[(k, (n' - k') * v / n'), (k + 1, k' * v / n')]
+
+breakdown :: (Num a, Eq a) => L a a -> Int
+breakdown f
+ | IM.null m = 50
+ | otherwise = fst (IM.findMin m) `min` (100 - fst (IM.findMax m))
+ where m = IM.filter (/= 0) $ callL f 101
View
9 Data/Pass/Ordered.hs
@@ -1,9 +0,0 @@
-module Data.Pass.Ordered where
-
-data Ordered k a b where
- Ord b => k a b -> (Int -> (Int -> b) -> c) -> c
-
--- data Ordered k a b where
--- Ordered :: Ord b => (o -> b) -> k a o -> Ordered k a b
--- OrderedAp :: (c -> d) -> Ordered k a (b -> c) -> Ordered k a b -> Ordered k a d
--- OrderedPure :: b -> Ordered k a b
View
88 Data/Pass/Robust.hs
@@ -0,0 +1,88 @@
+module Data.Pass.Robust
+ ( Robust(..)
+ , quantile
+ , tercile, t1, t2
+ , quartile, q1, q2, q3
+ , quintile, qu1, qu2, qu3, qu4
+ , percentile
+ , permille
+ ) where
+
+import Data.Pass.Type
+import Data.Pass.Calc
+import Data.Pass.Fun
+import Data.Pass.Thrist
+import Data.Pass.L
+import Data.Pass.Estimator
+
+-- | embedding for L-estimators
+class Robust l where
+ robust :: (Fractional a, Ord a) => L a a -> l a a
+
+ winsorized :: Rational -> L Double Double -> l Double Double
+ winsorized p f = robust $ Winsorized p f
+
+ trimmed :: Rational -> L Double Double -> l Double Double
+ trimmed p f = robust $ Trimmed p f
+
+ jackknifed :: L Double Double -> l Double Double
+ jackknifed f = robust $ Jackknifed f
+
+ median :: l Double Double
+ median = robust Median
+
+ lscale :: l Double Double
+ lscale = robust LScale
+
+ quantileBy :: Estimator -> Rational -> l Double Double
+ quantileBy e q = robust $ QuantileBy e q
+
+quantile :: Robust l => Rational -> l Double Double
+quantile = quantileBy R2
+
+tercile :: Robust l => Rational -> l Double Double
+tercile n = quantile (n/3)
+
+-- | terciles 1 and 2
+--
+-- > breakdown t1 = breakdown t2 = 33%
+t1, t2 :: Robust l => l Double Double
+t1 = quantile (1/3)
+t2 = quantile (2/3)
+
+quartile :: Robust l => Rational -> l Double Double
+quartile n = quantile (n/4)
+
+-- | quantiles, with breakdown points 25%, 50%, and 25% respectively
+q1, q2, q3 :: Robust l => l Double Double
+q1 = quantile 0.25
+q2 = median
+q3 = quantile 0.75
+
+quintile :: Robust l => Rational -> l Double Double
+quintile n = quantile (n/5)
+
+-- | quintiles 1 through 4
+qu1, qu2, qu3, qu4 :: Robust l => l Double Double
+qu1 = quintile 1
+qu2 = quintile 2
+qu3 = quintile 3
+qu4 = quintile 4
+
+percentile :: Robust l => Rational -> l Double Double
+percentile n = quantile (n/100)
+
+permille :: Robust l => Rational -> l Double Double
+permille n = quantile (n/1000)
+
+instance Robust L where
+ robust = id
+
+instance Robust l => Robust (Fun l) where
+ robust = Fun . robust
+
+instance Robust (Pass k) where
+ robust l = L id l Nil
+
+instance Robust (Calc k) where
+ robust l = Stop () `Step` \_ -> robust l
View
21 Data/Pass/Step.hs
@@ -1,10 +1,16 @@
module Data.Pass.Step
( Step(..)
+ , midhinge
+ , trimean
+ , iqr
+ , iqm
) where
import Data.Pass.Type
import Data.Pass.Calc
import Data.Pass.Prep
+import Data.Pass.Robust
+import Data.Pass.Accelerated
class Prep t => Step t where
step :: Pass k a b -> t k a b
@@ -14,3 +20,18 @@ instance Step Pass where
instance Step Calc where
step k = Step (Stop ()) (const k)
+
+midhinge :: Step t => t k Double Double
+midhinge = step $ (q1 + q3) / 2
+
+-- | Tukey's trimean
+trimean :: Step t => t k Double Double
+trimean = step $ (q1 + 2 * q2 + q3) / 4
+
+-- | interquartile range
+iqr :: Step t => t k Double Double
+iqr = step $ q3 - q1
+
+-- | interquartile mean
+iqm :: Step t => t k Double Double
+iqm = step $ trimmed 0.25 mean
View
7 Data/Pass/Util.hs
@@ -0,0 +1,7 @@
+module Data.Pass.Util where
+
+clamp :: Int -> Int -> Int
+clamp n k
+ | k <= 0 = 0
+ | k >= n = n - 1
+ | otherwise = k
View
19 examples/Pass.hs
@@ -16,25 +16,25 @@ import Data.Pass
-- example calculation type
data Test a b where
Total :: Num a => Test a (Sum a)
- Count :: Num b => Test a (Sum b)
- Square :: Num a => Test a a
+ Count :: Test a (Sum Int)
+ Square :: Test Double Double
Minus :: Double -> Test Double Double
- Abs :: Num a => Test a a
+ Abs :: Test Double Double
deriving Typeable
deriving instance Typeable1 Sum -- :(
-count :: (Step t, Typeable b, Num b) => t Test a b
-count = step $ getSum <$> trans Count
+count :: (Step t, Num b) => t Test a b
+count = step $ fromIntegral . getSum <$> trans Count
-sumSq :: (Step t, Typeable a, Num a) => t Test a a
-sumSq = prep Square total
+sumSq :: Step t => t Test Double Double
+sumSq = step $ prep Square total
-- E[X^2] - E[X]^2
-var :: (Step t, Typeable a, Fractional a) => t Test a a
+var :: Step t => t Test Double Double
var = step $ sumSq / count - mean ^ 2
-stddev :: (Step t, Typeable a, Floating a) => t Test a a
+stddev :: Step t => t Test Double Double
stddev = step $ sqrt var
-- > absDev median -- median absolute deviation
@@ -61,7 +61,6 @@ instance Named Test where
hashFunWithSalt n Square = n `hashWithSalt` 2
hashFunWithSalt n (Minus m) = n `hashWithSalt` 3 `hashWithSalt` m
hashFunWithSalt n Abs = n `hashWithSalt` 4
- hashFunWithSalt n Median = n `hashWithSalt` 5
instance Call Test where
call Total a = Sum a
View
4 multipass.cabal
@@ -39,6 +39,7 @@ library
Data.Pass.Calculation
Data.Pass.Call
Data.Pass.Class
+ Data.Pass.Estimator
Data.Pass.Env
Data.Pass.Eval
Data.Pass.Eval.Naive
@@ -53,4 +54,7 @@ library
Data.Pass.Trans
Data.Pass.Type
+ other-modules:
+ Data.Pass.Util
+
ghc-options: -Wall

0 comments on commit 600ddae

Please sign in to comment.
Something went wrong with that request. Please try again.