module Statistics.Resampling
(
Resample(..)
, jackknife
, resample
) where
import Control.Monad (forM_)
import Control.Monad.ST (ST)
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Unboxed.Mutable as MU
import Data.Vector.Unboxed ((!))
import Data.Vector.Generic (unsafeFreeze)
import Data.Vector.Algorithms.Intro (sort)
import Statistics.Function (createU, indices)
import System.Random.MWC (Gen, uniform)
import Statistics.Types (Estimator, Sample)
newtype Resample = Resample {
fromResample :: U.Vector Double
} deriving (Eq, Show)
resample :: Gen s -> [Estimator] -> Int -> Sample -> ST s [Resample]
resample gen ests numResamples samples = do
results <- mapM (const (MU.new numResamples)) $ ests
loop 0 (zip ests results)
mapM_ sort results
mapM (fmap Resample . unsafeFreeze) results
where
loop k ers | k >= numResamples = return ()
| otherwise = do
re <- createU n $ \_ -> do
r <- uniform gen
return (samples ! (abs r `mod` n))
forM_ ers $ \(est,arr) ->
MU.write arr k . est $ re
loop (k+1) ers
n = U.length samples
jackknife :: Estimator -> Sample -> U.Vector Double
jackknife est sample = U.map f . indices $ sample
where f i = est (dropAt i sample)
indexed :: U.Unbox e => U.Vector e -> U.Vector (Int,e)
indexed a = U.zip (U.enumFromN 0 (U.length a)) a
dropAt :: U.Unbox e => Int -> U.Vector e -> U.Vector e
dropAt n = U.map snd . U.filter notN . indexed
where notN (i , _) = i /= n