1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
|
-- |
-- Module : Statistics.Resampling
-- Copyright : (c) 2009, 2010 Bryan O'Sullivan
-- License : BSD3
--
-- Maintainer : bos@serpentine.com
-- Stability : experimental
-- Portability : portable
--
-- Resampling statistics.
module Statistics.Resampling
(
Resample(..)
, jackknife
, resample
) where
import Control.Monad (forM_, liftM)
import Control.Monad.Primitive (PrimMonad, PrimState)
import Data.Vector.Algorithms.Intro (sort)
import Data.Vector.Generic (unsafeFreeze)
import Data.Vector.Unboxed ((!))
import Statistics.Function (create, indexed, indices)
import Statistics.Types (Estimator, Sample)
import System.Random.MWC (Gen, uniform)
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Unboxed.Mutable as MU
-- | A resample drawn randomly, with replacement, from a set of data
-- points. Distinct from a normal array to make it harder for your
-- humble author's brain to go wrong.
newtype Resample = Resample {
fromResample :: U.Vector Double
} deriving (Eq, Show)
-- | Resample a data set repeatedly, with replacement, computing each
-- estimate over the resampled data.
resample :: (PrimMonad m) => Gen (PrimState m) -> [Estimator] -> Int -> Sample -> m [Resample]
resample gen ests numResamples samples = do
results <- mapM (const (MU.new numResamples)) $ ests
loop 0 (zip ests results)
mapM_ sort results
mapM (liftM Resample . unsafeFreeze) results
where
loop k ers | k >= numResamples = return ()
| otherwise = do
re <- create 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
-- | Compute a statistical estimate repeatedly over a sample, each
-- time omitting a successive element.
jackknife :: Estimator -> Sample -> U.Vector Double
jackknife est sample = U.map f . indices $ sample
where f i = est (dropAt i sample)
{-# INLINE jackknife #-}
-- | Drop the /k/th element of a vector.
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
|
