module Statistics.Resampling
(
Resample(..)
, jackknife
, resample
) where
import Control.Monad (forM_)
import Control.Monad.ST (unsafeSTToIO)
import Data.Array.Vector
import Data.Array.Vector.Algorithms.Intro (sort)
import Statistics.Function (createU)
import Statistics.Types (Estimator, Sample)
import System.Random.Mersenne (MTGen, random)
newtype Resample = Resample {
fromResample :: UArr Double
} deriving (Eq, Show)
resample :: MTGen -> [Estimator] -> Int -> Sample -> IO [Resample]
resample gen ests numResamples samples = do
results <- unsafeSTToIO . mapM (const (newMU numResamples)) $ ests
loop 0 (zip ests results)
unsafeSTToIO $ do
mapM_ sort results
mapM (fmap Resample . unsafeFreezeAllMU) results
where
loop k ers | k >= numResamples = return ()
| otherwise = do
re <- createU n $ \_ -> do
r <- random gen
return (indexU samples (abs r `mod` n))
unsafeSTToIO . forM_ ers $ \(est,arr) ->
writeMU arr k . est $ re
loop (k+1) ers
n = lengthU samples
jackknife :: Estimator -> Sample -> UArr Double
jackknife est sample = mapU f . enumFromToU 0 . subtract 1 . lengthU $ sample
where f i = est (dropAt i sample)
dropAt :: UA e => Int -> UArr e -> UArr e
dropAt n = mapU sndT . filterU notN . indexedU
where notN (i :*: _) = i /= n
sndT (_ :*: k) = k