-- |
-- 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/>