    -- |
    -- Module    : Statistics.Resampling.Bootstrap
    -- Copyright : (c) 2009, 2011 Bryan O'Sullivan
    -- License   : BSD3
    --
    -- Maintainer  : bos@serpentine.com
    -- Stability   : experimental
    -- Portability : portable
    --
    -- The bootstrap method for statistical inference.

    module Statistics.Resampling.Bootstrap
        ( bootstrapBCA
        , basicBootstrap
        -- * References
        -- $references
        ) where

    import           Data.Vector.Generic ((!))
    import qualified Data.Vector.Unboxed as U
    import qualified Data.Vector.Generic as G

    import Statistics.Distribution (cumulative, quantile)
    import Statistics.Distribution.Normal
    import Statistics.Resampling (Bootstrap(..), jackknife)
    import Statistics.Sample (mean)
    import Statistics.Types (Sample, CL, Estimate, ConfInt, estimateFromInterval,
                             estimateFromErr, CL, significanceLevel)
    import Statistics.Function (gsort)

    import qualified Statistics.Resampling as R

    import Control.Parallel.Strategies (parMap, rdeepseq)

    data T = {-# UNPACK #-} !Double :< {-# UNPACK #-} !Double
    infixl 2 :<

    -- | Bias-corrected accelerated (BCA) bootstrap. This adjusts for both
    --   bias and skewness in the resampled distribution.
    --
    --   BCA algorithm is described in ch. 5 of Davison, Hinkley "Confidence
    --   intervals" in section 5.3 "Percentile method"
    bootstrapBCA
      :: CL Double       -- ^ Confidence level
      -> Sample          -- ^ Full data sample
      -> [(R.Estimator, Bootstrap U.Vector Double)]
      -- ^ Estimates obtained from resampled data and estimator used for
      --   this.
      -> [Estimate ConfInt Double]
    bootstrapBCA :: CL Double
    -> Sample
    -> [(Estimator, Bootstrap Vector Double)]
    -> [Estimate ConfInt Double]
    bootstrapBCA CL Double
    confidenceLevel Sample
    sample [(Estimator, Bootstrap Vector Double)]
    resampledData
      = Strategy (Estimate ConfInt Double)
    -> ((Estimator, Bootstrap Vector Double)
        -> Estimate ConfInt Double)
    -> [(Estimator, Bootstrap Vector Double)]
    -> [Estimate ConfInt Double]
    forall b a. Strategy b -> (a -> b) -> [a] -> [b]
    parMap Strategy (Estimate ConfInt Double)
    forall a. NFData a => Strategy a
    rdeepseq (Estimator, Bootstrap Vector Double) -> Estimate ConfInt Double
    forall {a}.
    (Num a, Unbox a, Ord a) =>
    (Estimator, Bootstrap Vector a) -> Estimate ConfInt a
    e [(Estimator, Bootstrap Vector Double)]
    resampledData
      where
        e :: (Estimator, Bootstrap Vector a) -> Estimate ConfInt a
    e (Estimator
    est, Bootstrap a
    pt Vector a
    resample)
          | Sample -> Int
    forall a. Unbox a => Vector a -> Int
    U.length Sample
    sample Int -> Int -> Bool
    forall a. Eq a => a -> a -> Bool
    == Int
    1 Bool -> Bool -> Bool
    || Double -> Bool
    forall a. RealFloat a => a -> Bool
    isInfinite Double
    bias =
              a -> (a, a) -> CL Double -> Estimate ConfInt a
    forall a. a -> (a, a) -> CL Double -> Estimate ConfInt a
    estimateFromErr      a
    pt (a
    0,a
    0) CL Double
    confidenceLevel
          | Bool
    otherwise =
              a -> (a, a) -> CL Double -> Estimate ConfInt a
    forall a. Num a => a -> (a, a) -> CL Double -> Estimate ConfInt a
    estimateFromInterval a
    pt (Vector a
    resample Vector a -> Int -> a
    forall (v :: * -> *) a.
    (HasCallStack, Vector v a) =>
    v a -> Int -> a
    ! Int
    lo, Vector a
    resample Vector a -> Int -> a
    forall (v :: * -> *) a.
    (HasCallStack, Vector v a) =>
    v a -> Int -> a
    ! Int
    hi) CL Double
    confidenceLevel
          where
            -- Quantile estimates for given CL
            lo :: Int
    lo    = Int -> Int -> Int
    forall a. Ord a => a -> a -> a
    min (Int -> Int -> Int
    forall a. Ord a => a -> a -> a
    max (Double -> Int
    cumn Double
    a1) Int
    0) (Int
    ni Int -> Int -> Int
    forall a. Num a => a -> a -> a
    - Int
    1)
              where a1 :: Double
    a1 = Double
    bias Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + Double
    b1 Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / (Double
    1 Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    accel Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    b1)
                    b1 :: Double
    b1 = Double
    bias Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + Double
    z1
            hi :: Int
    hi    = Int -> Int -> Int
    forall a. Ord a => a -> a -> a
    max (Int -> Int -> Int
    forall a. Ord a => a -> a -> a
    min (Double -> Int
    cumn Double
    a2) (Int
    ni Int -> Int -> Int
    forall a. Num a => a -> a -> a
    - Int
    1)) Int
    0
              where a2 :: Double
    a2 = Double
    bias Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + Double
    b2 Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / (Double
    1 Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    accel Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    b2)
                    b2 :: Double
    b2 = Double
    bias Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    z1
            -- Number of resamples
            ni :: Int
    ni    = Vector a -> Int
    forall a. Unbox a => Vector a -> Int
    U.length Vector a
    resample
            n :: Double
    n     = Int -> Double
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral Int
    ni
            -- Corrections
            z1 :: Double
    z1    = NormalDistribution -> Double -> Double
    forall d. ContDistr d => d -> Double -> Double
    quantile NormalDistribution
    standard (CL Double -> Double
    forall a. CL a -> a
    significanceLevel CL Double
    confidenceLevel Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / Double
    2)
            cumn :: Double -> Int
    cumn  = Double -> Int
    forall b. Integral b => Double -> b
    forall a b. (RealFrac a, Integral b) => a -> b
    round (Double -> Int) -> (Double -> Double) -> Double -> Int
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (Double -> Double -> Double
    forall a. Num a => a -> a -> a
    *Double
    n) (Double -> Double) -> (Double -> Double) -> Double -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . NormalDistribution -> Double -> Double
    forall d. Distribution d => d -> Double -> Double
    cumulative NormalDistribution
    standard
            bias :: Double
    bias  = NormalDistribution -> Double -> Double
    forall d. ContDistr d => d -> Double -> Double
    quantile NormalDistribution
    standard (Double
    probN Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / Double
    n)
              where probN :: Double
    probN = Int -> Double
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral (Int -> Double) -> (Vector a -> Int) -> Vector a -> Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector a -> Int
    forall a. Unbox a => Vector a -> Int
    U.length (Vector a -> Int) -> (Vector a -> Vector a) -> Vector a -> Int
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (a -> Bool) -> Vector a -> Vector a
    forall a. Unbox a => (a -> Bool) -> Vector a -> Vector a
    U.filter (a -> a -> Bool
    forall a. Ord a => a -> a -> Bool
    <a
    pt) (Vector a -> Double) -> Vector a -> Double
    forall a b. (a -> b) -> a -> b
    $ Vector a
    resample
            accel :: Double
    accel = Double
    sumCubes Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / (Double
    6 Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * (Double
    sumSquares Double -> Double -> Double
    forall a. Floating a => a -> a -> a
    ** Double
    1.5))
              where (Double
    sumSquares :< Double
    sumCubes) = (T -> Double -> T) -> T -> Sample -> T
    forall b a. Unbox b => (a -> b -> a) -> a -> Vector b -> a
    U.foldl' T -> Double -> T
    f (Double
    0 Double -> Double -> T
    :< Double
    0) Sample
    jack
                    f :: T -> Double -> T
    f (Double
    s :< Double
    c) Double
    j = Double
    s Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + Double
    d2 Double -> Double -> T
    :< Double
    c Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + Double
    d2 Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    d
                        where d :: Double
    d  = Double
    jackMean Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    j
                              d2 :: Double
    d2 = Double
    d Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    d
                    jackMean :: Double
    jackMean     = Sample -> Double
    forall (v :: * -> *). Vector v Double => v Double -> Double
    mean Sample
    jack
            jack :: Sample
    jack  = Estimator -> Sample -> Sample
    jackknife Estimator
    est Sample
    sample

    -- | Basic bootstrap. This method simply uses empirical quantiles for
    --   confidence interval.
    basicBootstrap
      :: (G.Vector v a, Ord a, Num a)
      => CL Double       -- ^ Confidence vector
      -> Bootstrap v a   -- ^ Estimate from full sample and vector of
                         --   estimates obtained from resamples
      -> Estimate ConfInt a
    {-# INLINE basicBootstrap #-}
    basicBootstrap :: forall (v :: * -> *) a.
    (Vector v a, Ord a, Num a) =>
    CL Double -> Bootstrap v a -> Estimate ConfInt a
    basicBootstrap CL Double
    cl (Bootstrap a
    e v a
    ests)
      = a -> (a, a) -> CL Double -> Estimate ConfInt a
    forall a. Num a => a -> (a, a) -> CL Double -> Estimate ConfInt a
    estimateFromInterval a
    e (v a
    sorted v a -> Int -> a
    forall (v :: * -> *) a.
    (HasCallStack, Vector v a) =>
    v a -> Int -> a
    ! Int
    lo, v a
    sorted v a -> Int -> a
    forall (v :: * -> *) a.
    (HasCallStack, Vector v a) =>
    v a -> Int -> a
    ! Int
    hi) CL Double
    cl
      where
        sorted :: v a
    sorted = v a -> v a
    forall e (v :: * -> *). (Ord e, Vector v e) => v e -> v e
    gsort v a
    ests
        n :: Double
    n  = Int -> Double
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral (Int -> Double) -> Int -> Double
    forall a b. (a -> b) -> a -> b
    $ v a -> Int
    forall (v :: * -> *) a. Vector v a => v a -> Int
    G.length v a
    ests
        c :: Double
    c  = Double
    n Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * (CL Double -> Double
    forall a. CL a -> a
    significanceLevel CL Double
    cl Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / Double
    2)
        -- FIXME: can we have better estimates of quantiles in case when p
        --        is not multiple of 1/N
        --
        -- FIXME: we could have undercoverage here
        lo :: Int
    lo = Double -> Int
    forall b. Integral b => Double -> b
    forall a b. (RealFrac a, Integral b) => a -> b
    round Double
    c
        hi :: Int
    hi = Double -> Int
    forall b. Integral b => Double -> b
    forall a b. (RealFrac a, Integral b) => a -> b
    truncate (Double
    n Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    c)

    -- $references
    --
    -- * Davison, A.C; Hinkley, D.V. (1997) Bootstrap methods and their
    --   application. <http://statwww.epfl.ch/davison/BMA/>