module Statistics.Resampling.Bootstrap
(
Estimate(..)
, bootstrapBCA
) 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)
data Estimate = Estimate {
estPoint :: !Double
, estLowerBound :: !Double
, estUpperBound :: !Double
, estConfidenceLevel :: !Double
} 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 = !Double :< !Double
infixl 2 :<
bootstrapBCA :: Double
-> Sample
-> [Estimator]
-> [Resample]
-> [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 (<pt) $ resample
ni = lengthU resample
n = fromIntegral ni
accel = sumCubes / (6 * (sumSquares ** 1.5))
where (sumSquares :< sumCubes) = foldlU 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