-- | -- Module : Statistics.Resampling.Bootstrap -- Copyright : (c) 2009 Bryan O'Sullivan -- License : BSD3 -- -- Maintainer : bos@serpentine.com -- Stability : experimental -- Portability : portable -- -- The bootstrap method for statistical inference. module Statistics.Resampling.Bootstrap ( Estimate(..) , bootstrapBCA -- * References -- \$references ) where import Control.Exception (assert) import Data.Array.Vector (foldlU, filterU, indexU, lengthU) import Statistics.Distribution.Normal import Statistics.Distribution (cumulative, quantile) import Statistics.Resampling (Resample(..), jackknife) import Statistics.Sample (mean) import Statistics.Types (Estimator, Sample) -- | A point and interval estimate computed via an 'Estimator'. data Estimate = Estimate { estPoint :: {-# UNPACK #-} !Double -- ^ Point estimate. , estLowerBound :: {-# UNPACK #-} !Double -- ^ Lower bound of the estimate interval (i.e. the lower bound of -- the confidence interval). , estUpperBound :: {-# UNPACK #-} !Double -- ^ Upper bound of the estimate interval (i.e. the upper bound of -- the confidence interval). , estConfidenceLevel :: {-# UNPACK #-} !Double -- ^ Confidence level of the confidence intervals. } deriving (Eq, Show) 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 = {-# UNPACK #-} !Double :< {-# UNPACK #-} !Double infixl 2 :< -- | Bias-corrected accelerated (BCA) bootstrap. This adjusts for both -- bias and skewness in the resampled distribution. bootstrapBCA :: Double -- ^ Confidence level -> Sample -- ^ Sample data -> [Estimator] -- ^ Estimators -> [Resample] -- ^ Resampled data -> [Estimate] bootstrapBCA confidenceLevel sample = assert (confidenceLevel > 0 && confidenceLevel < 1) zipWith e where e est (Resample resample) | lengthU sample == 1 = estimate pt pt pt confidenceLevel | otherwise = estimate pt (indexU resample lo) (indexU resample hi) confidenceLevel where pt = 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 . lengthU . filterU (