module Statistics.Resampling.Bootstrap
(
Estimate(..)
, bootstrapBCA
, scale
) where
import Control.Applicative ((<$>), (<*>))
import Control.DeepSeq (NFData)
import Control.Exception (assert)
import Control.Monad.Par (parMap, runPar)
import Data.Binary (Binary)
import Data.Binary (put, get)
import Data.Data (Data)
import Data.Typeable (Typeable)
import Data.Vector.Unboxed ((!))
import GHC.Generics (Generic)
import Statistics.Distribution (cumulative, quantile)
import Statistics.Distribution.Normal
import Statistics.Resampling (Resample(..), jackknife)
import Statistics.Sample (mean)
import Statistics.Types (Estimator, Sample)
import qualified Data.Vector.Unboxed as U
import qualified Statistics.Resampling as R
data Estimate = Estimate {
estPoint :: !Double
, estLowerBound :: !Double
, estUpperBound :: !Double
, estConfidenceLevel :: !Double
} deriving (Eq, Read, Show, Typeable, Data, Generic)
instance Binary Estimate where
put (Estimate w x y z) = put w >> put x >> put y >> put z
get = Estimate <$> get <*> get <*> get <*> get
instance NFData Estimate
scale :: Double
-> Estimate -> Estimate
scale f e@Estimate{..} = e {
estPoint = f * estPoint
, estLowerBound = f * estLowerBound
, estUpperBound = f * estUpperBound
}
estimate :: Double -> Double -> Double -> Double -> Estimate
estimate pt lb ub cl =
assert (lb <= ub) .
assert (cl > 0 && cl < 1) $
Estimate { estPoint = pt
, estLowerBound = lb
, estUpperBound = ub
, estConfidenceLevel = cl
}
data T = !Double :< !Double
infixl 2 :<
bootstrapBCA :: Double
-> Sample
-> [Estimator]
-> [Resample]
-> [Estimate]
bootstrapBCA confidenceLevel sample estimators resamples
| confidenceLevel > 0 && confidenceLevel < 1
= runPar $ parMap (uncurry e) (zip estimators resamples)
| otherwise = error "Statistics.Resampling.Bootstrap.bootstrapBCA: confidence level outside (0,1) range"
where
e est (Resample resample)
| U.length sample == 1 = estimate pt pt pt confidenceLevel
| otherwise =
estimate pt (resample ! lo) (resample ! hi) confidenceLevel
where
pt = R.estimate est sample
lo = max (cumn a1) 0
where a1 = bias + b1 / (1 accel * b1)
b1 = bias + z1
hi = min (cumn a2) (ni 1)
where a2 = bias + b2 / (1 accel * b2)
b2 = bias z1
z1 = quantile standard ((1 confidenceLevel) / 2)
cumn = round . (*n) . cumulative standard
bias = quantile standard (probN / n)
where probN = fromIntegral . U.length . U.filter (<pt) $ resample
ni = U.length resample
n = fromIntegral ni
accel = sumCubes / (6 * (sumSquares ** 1.5))
where (sumSquares :< sumCubes) = U.foldl' f (0 :< 0) jack
f (s :< c) j = s + d2 :< c + d2 * d
where d = jackMean j
d2 = d * d
jackMean = mean jack
jack = jackknife est sample