-- |
-- Module    : Statistics.Regression
-- Copyright : 2014 Bryan O'Sullivan
-- License   : BSD3
--
-- Functions for regression analysis.

module Statistics.Regression
    (
      olsRegress
    , ols
    , rSquare
    , bootstrapRegress
    ) where

import Control.Concurrent.Async (forConcurrently)
import Control.DeepSeq (rnf)
import Control.Monad (when)
import Data.List (nub)
import GHC.Conc (getNumCapabilities)
import Prelude hiding (pred, sum)
import Statistics.Function as F
import Statistics.Matrix hiding (map)
import Statistics.Matrix.Algorithms (qr)
import Statistics.Resampling (splitGen)
import Statistics.Types      (Estimate(..),ConfInt,CL,estimateFromInterval,significanceLevel)
import Statistics.Sample (mean)
import Statistics.Sample.Internal (sum)
import System.Random.MWC (GenIO, uniformR)
import qualified Data.Vector as V
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Unboxed.Mutable as M

-- | Perform an ordinary least-squares regression on a set of
-- predictors, and calculate the goodness-of-fit of the regression.
--
-- The returned pair consists of:
--
-- * A vector of regression coefficients.  This vector has /one more/
--   element than the list of predictors; the last element is the
--   /y/-intercept value.
--
-- * /R²/, the coefficient of determination (see 'rSquare' for
--   details).
olsRegress :: [Vector]
              -- ^ Non-empty list of predictor vectors.  Must all have
              -- the same length.  These will become the columns of
              -- the matrix /A/ solved by 'ols'.
           -> Vector
              -- ^ Responder vector.  Must have the same length as the
              -- predictor vectors.
           -> (Vector, Double)
olsRegress :: [Vector] -> Vector -> (Vector, Double)
olsRegress preds :: [Vector]
preds@(Vector
_:[Vector]
_) Vector
resps
  | (Int -> Bool) -> [Int] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/=Int
n) [Int]
ls        = [Char] -> (Vector, Double)
forall a. HasCallStack => [Char] -> a
error ([Char] -> (Vector, Double)) -> [Char] -> (Vector, Double)
forall a b. (a -> b) -> a -> b
$ [Char]
"predictor vector length mismatch " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++
                                  [Int] -> [Char]
forall a. Show a => a -> [Char]
show [Int]
lss
  | Vector -> Int
forall (v :: * -> *) a. Vector v a => v a -> Int
G.length Vector
resps Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n = [Char] -> (Vector, Double)
forall a. HasCallStack => [Char] -> a
error ([Char] -> (Vector, Double)) -> [Char] -> (Vector, Double)
forall a b. (a -> b) -> a -> b
$ [Char]
"responder/predictor length mismatch " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++
                                  (Int, Int) -> [Char]
forall a. Show a => a -> [Char]
show (Vector -> Int
forall (v :: * -> *) a. Vector v a => v a -> Int
G.length Vector
resps, Int
n)
  | Bool
otherwise           = (Vector
coeffs, Matrix -> Vector -> Vector -> Double
rSquare Matrix
mxpreds Vector
resps Vector
coeffs)
  where
    coeffs :: Vector
coeffs    = Matrix -> Vector -> Vector
ols Matrix
mxpreds Vector
resps
    mxpreds :: Matrix
mxpreds   = Matrix -> Matrix
transpose (Matrix -> Matrix) -> ([Vector] -> Matrix) -> [Vector] -> Matrix
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
                Int -> Int -> Vector -> Matrix
fromVector ([Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
lss Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int
n (Vector -> Matrix) -> ([Vector] -> Vector) -> [Vector] -> Matrix
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
                [Vector] -> Vector
forall (v :: * -> *) a. Vector v a => [v a] -> v a
G.concat ([Vector] -> Matrix) -> [Vector] -> Matrix
forall a b. (a -> b) -> a -> b
$ [Vector]
preds [Vector] -> [Vector] -> [Vector]
forall a. [a] -> [a] -> [a]
++ [Int -> Double -> Vector
forall (v :: * -> *) a. Vector v a => Int -> a -> v a
G.replicate Int
n Double
1]
    lss :: [Int]
lss@(Int
n:[Int]
ls) = (Vector -> Int) -> [Vector] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map Vector -> Int
forall (v :: * -> *) a. Vector v a => v a -> Int
G.length [Vector]
preds
olsRegress [Vector]
_ Vector
_ = [Char] -> (Vector, Double)
forall a. HasCallStack => [Char] -> a
error [Char]
"no predictors given"

-- | Compute the ordinary least-squares solution to /A x = b/.
ols :: Matrix     -- ^ /A/ has at least as many rows as columns.
    -> Vector     -- ^ /b/ has the same length as columns in /A/.
    -> Vector
ols :: Matrix -> Vector -> Vector
ols Matrix
a Vector
b
  | Int
rs Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
cs   = [Char] -> Vector
forall a. HasCallStack => [Char] -> a
error ([Char] -> Vector) -> [Char] -> Vector
forall a b. (a -> b) -> a -> b
$ [Char]
"fewer rows than columns " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ (Int, Int) -> [Char]
forall a. Show a => a -> [Char]
show (Int, Int)
d
  | Bool
otherwise = Matrix -> Vector -> Vector
solve Matrix
r (Matrix -> Matrix
transpose Matrix
q Matrix -> Vector -> Vector
`multiplyV` Vector
b)
  where
    d :: (Int, Int)
d@(Int
rs,Int
cs) = Matrix -> (Int, Int)
dimension Matrix
a
    (Matrix
q,Matrix
r)     = Matrix -> (Matrix, Matrix)
qr Matrix
a

-- | Solve the equation /R x = b/.
solve :: Matrix     -- ^ /R/ is an upper-triangular square matrix.
      -> Vector     -- ^ /b/ is of the same length as rows\/columns in /R/.
      -> Vector
solve :: Matrix -> Vector -> Vector
solve Matrix
r Vector
b
  | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
l    = [Char] -> Vector
forall a. HasCallStack => [Char] -> a
error ([Char] -> Vector) -> [Char] -> Vector
forall a b. (a -> b) -> a -> b
$ [Char]
"row/vector mismatch " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ (Int, Int) -> [Char]
forall a. Show a => a -> [Char]
show (Int
n,Int
l)
  | Bool
otherwise = (forall s. ST s (MVector s Double)) -> Vector
forall a. Unbox a => (forall s. ST s (MVector s a)) -> Vector a
U.create ((forall s. ST s (MVector s Double)) -> Vector)
-> (forall s. ST s (MVector s Double)) -> Vector
forall a b. (a -> b) -> a -> b
$ do
  MVector s Double
s <- Vector -> ST s (MVector (PrimState (ST s)) Double)
forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
Vector a -> m (MVector (PrimState m) a)
U.thaw Vector
b
  Int -> Int -> (Int -> ST s ()) -> ST s ()
forall (m :: * -> *).
Monad m =>
Int -> Int -> (Int -> m ()) -> m ()
rfor Int
n Int
0 ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
i -> do
    Double
si <- (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Matrix -> Int -> Int -> Double
unsafeIndex Matrix
r Int
i Int
i) (Double -> Double) -> ST s Double -> ST s Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> MVector (PrimState (ST s)) Double -> Int -> ST s Double
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
M.unsafeRead MVector s Double
MVector (PrimState (ST s)) Double
s Int
i
    MVector (PrimState (ST s)) Double -> Int -> Double -> ST s ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
M.unsafeWrite MVector s Double
MVector (PrimState (ST s)) Double
s Int
i Double
si
    Int -> Int -> (Int -> ST s ()) -> ST s ()
forall (m :: * -> *).
Monad m =>
Int -> Int -> (Int -> m ()) -> m ()
F.for Int
0 Int
i ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
j -> MVector s Double -> Int -> (Double -> Double) -> ST s ()
forall s. MVector s Double -> Int -> (Double -> Double) -> ST s ()
F.unsafeModify MVector s Double
s Int
j ((Double -> Double) -> ST s ()) -> (Double -> Double) -> ST s ()
forall a b. (a -> b) -> a -> b
$ Double -> Double -> Double
forall a. Num a => a -> a -> a
subtract (Matrix -> Int -> Int -> Double
unsafeIndex Matrix
r Int
j Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
si)
  MVector s Double -> ST s (MVector s Double)
forall (m :: * -> *) a. Monad m => a -> m a
return MVector s Double
s
  where n :: Int
n = Matrix -> Int
rows Matrix
r
        l :: Int
l = Vector -> Int
forall a. Unbox a => Vector a -> Int
U.length Vector
b

-- | Compute /R&#0178;/, the coefficient of determination that
-- indicates goodness-of-fit of a regression.
--
-- This value will be 1 if the predictors fit perfectly, dropping to 0
-- if they have no explanatory power.
rSquare :: Matrix               -- ^ Predictors (regressors).
        -> Vector               -- ^ Responders.
        -> Vector               -- ^ Regression coefficients.
        -> Double
rSquare :: Matrix -> Vector -> Vector -> Double
rSquare Matrix
pred Vector
resp Vector
coeff = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
r Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
t
  where
    r :: Double
r   = Vector -> Double
forall (v :: * -> *). Vector v Double => v Double -> Double
sum (Vector -> Double) -> Vector -> Double
forall a b. (a -> b) -> a -> b
$ ((Int -> Double -> Double) -> Vector -> Vector)
-> Vector -> (Int -> Double -> Double) -> Vector
forall a b c. (a -> b -> c) -> b -> a -> c
flip (Int -> Double -> Double) -> Vector -> Vector
forall a b.
(Unbox a, Unbox b) =>
(Int -> a -> b) -> Vector a -> Vector b
U.imap Vector
resp ((Int -> Double -> Double) -> Vector)
-> (Int -> Double -> Double) -> Vector
forall a b. (a -> b) -> a -> b
$ \Int
i Double
x -> Double -> Double
square (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Int -> Double
p Int
i)
    t :: Double
t   = Vector -> Double
forall (v :: * -> *). Vector v Double => v Double -> Double
sum (Vector -> Double) -> Vector -> Double
forall a b. (a -> b) -> a -> b
$ ((Double -> Double) -> Vector -> Vector)
-> Vector -> (Double -> Double) -> Vector
forall a b c. (a -> b -> c) -> b -> a -> c
flip (Double -> Double) -> Vector -> Vector
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
U.map Vector
resp ((Double -> Double) -> Vector) -> (Double -> Double) -> Vector
forall a b. (a -> b) -> a -> b
$ \Double
x -> Double -> Double
square (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Vector -> Double
forall (v :: * -> *). Vector v Double => v Double -> Double
mean Vector
resp)
    p :: Int -> Double
p Int
i = Vector -> Double
forall (v :: * -> *). Vector v Double => v Double -> Double
sum (Vector -> Double)
-> ((Int -> Double -> Double) -> Vector)
-> (Int -> Double -> Double)
-> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Int -> Double -> Double) -> Vector -> Vector)
-> Vector -> (Int -> Double -> Double) -> Vector
forall a b c. (a -> b -> c) -> b -> a -> c
flip (Int -> Double -> Double) -> Vector -> Vector
forall a b.
(Unbox a, Unbox b) =>
(Int -> a -> b) -> Vector a -> Vector b
U.imap Vector
coeff ((Int -> Double -> Double) -> Double)
-> (Int -> Double -> Double) -> Double
forall a b. (a -> b) -> a -> b
$ \Int
j -> (Double -> Double -> Double
forall a. Num a => a -> a -> a
* Matrix -> Int -> Int -> Double
unsafeIndex Matrix
pred Int
i Int
j)

-- | Bootstrap a regression function.  Returns both the results of the
-- regression and the requested confidence interval values.
bootstrapRegress
  :: GenIO
  -> Int         -- ^ Number of resamples to compute.
  -> CL Double   -- ^ Confidence level.
  -> ([Vector] -> Vector -> (Vector, Double))
     -- ^ Regression function.
  -> [Vector]    -- ^ Predictor vectors.
  -> Vector      -- ^ Responder vector.
  -> IO (V.Vector (Estimate ConfInt Double), Estimate ConfInt Double)
bootstrapRegress :: GenIO
-> Int
-> CL Double
-> ([Vector] -> Vector -> (Vector, Double))
-> [Vector]
-> Vector
-> IO (Vector (Estimate ConfInt Double), Estimate ConfInt Double)
bootstrapRegress GenIO
gen0 Int
numResamples CL Double
cl [Vector] -> Vector -> (Vector, Double)
rgrss [Vector]
preds0 Vector
resp0
  | Int
numResamples Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
1   = [Char]
-> IO (Vector (Estimate ConfInt Double), Estimate ConfInt Double)
forall a. HasCallStack => [Char] -> a
error ([Char]
 -> IO (Vector (Estimate ConfInt Double), Estimate ConfInt Double))
-> [Char]
-> IO (Vector (Estimate ConfInt Double), Estimate ConfInt Double)
forall a b. (a -> b) -> a -> b
$ [Char]
"bootstrapRegress: number of resamples " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++
                                 [Char]
"must be positive"
  | Bool
otherwise = do

  -- some error checks so that we do not run into vector index out of bounds.
  case [Int] -> [Int]
forall a. Eq a => [a] -> [a]
nub ((Vector -> Int) -> [Vector] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map Vector -> Int
forall a. Unbox a => Vector a -> Int
U.length [Vector]
preds0) of
    [] -> [Char] -> IO ()
forall a. HasCallStack => [Char] -> a
error [Char]
"bootstrapRegress: predictor vectors must not be empty"
    [Int
plen] -> do
        let rlen :: Int
rlen = Vector -> Int
forall a. Unbox a => Vector a -> Int
U.length Vector
resp0
        Bool -> IO () -> IO ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
plen Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
rlen) (IO () -> IO ()) -> IO () -> IO ()
forall a b. (a -> b) -> a -> b
$
            [Char] -> IO ()
forall a. HasCallStack => [Char] -> a
error ([Char] -> IO ()) -> [Char] -> IO ()
forall a b. (a -> b) -> a -> b
$ [Char]
"bootstrapRegress: responder vector length ["
                [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Int -> [Char]
forall a. Show a => a -> [Char]
show Int
rlen
                [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
"] must be the same as predictor vectors' length ["
                [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Int -> [Char]
forall a. Show a => a -> [Char]
show Int
plen [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
"]"
    [Int]
xs -> [Char] -> IO ()
forall a. HasCallStack => [Char] -> a
error ([Char] -> IO ()) -> [Char] -> IO ()
forall a b. (a -> b) -> a -> b
$ [Char]
"bootstrapRegress: all predictor vectors must be of the same \
        \length, lengths provided are: " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Int] -> [Char]
forall a. Show a => a -> [Char]
show [Int]
xs

  Int
caps <- IO Int
getNumCapabilities
  [Gen RealWorld]
gens <- Int -> GenIO -> IO [GenIO]
splitGen Int
caps GenIO
gen0
  [Vector (Vector, Double)]
vs <- [(Gen RealWorld, Int)]
-> ((Gen RealWorld, Int) -> IO (Vector (Vector, Double)))
-> IO [Vector (Vector, Double)]
forall (t :: * -> *) a b.
Traversable t =>
t a -> (a -> IO b) -> IO (t b)
forConcurrently ([Gen RealWorld] -> [Int] -> [(Gen RealWorld, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Gen RealWorld]
gens (Int -> Int -> [Int]
balance Int
caps Int
numResamples)) (((Gen RealWorld, Int) -> IO (Vector (Vector, Double)))
 -> IO [Vector (Vector, Double)])
-> ((Gen RealWorld, Int) -> IO (Vector (Vector, Double)))
-> IO [Vector (Vector, Double)]
forall a b. (a -> b) -> a -> b
$ \(Gen RealWorld
gen,Int
count) -> do
      Vector (Vector, Double)
v <- Int -> IO (Vector, Double) -> IO (Vector (Vector, Double))
forall (m :: * -> *) a. Monad m => Int -> m a -> m (Vector a)
V.replicateM Int
count (IO (Vector, Double) -> IO (Vector (Vector, Double)))
-> IO (Vector, Double) -> IO (Vector (Vector, Double))
forall a b. (a -> b) -> a -> b
$ do
           let n :: Int
n = Vector -> Int
forall a. Unbox a => Vector a -> Int
U.length Vector
resp0
           Vector Int
ixs <- Int -> IO Int -> IO (Vector Int)
forall (m :: * -> *) a.
(Monad m, Unbox a) =>
Int -> m a -> m (Vector a)
U.replicateM Int
n (IO Int -> IO (Vector Int)) -> IO Int -> IO (Vector Int)
forall a b. (a -> b) -> a -> b
$ (Int, Int) -> GenIO -> IO Int
forall a (m :: * -> *).
(Variate a, PrimMonad m) =>
(a, a) -> Gen (PrimState m) -> m a
uniformR (Int
0,Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Gen RealWorld
GenIO
gen
           let resp :: Vector
resp  = Vector -> Vector Int -> Vector
forall a. Unbox a => Vector a -> Vector Int -> Vector a
U.backpermute Vector
resp0 Vector Int
ixs
               preds :: [Vector]
preds = (Vector -> Vector) -> [Vector] -> [Vector]
forall a b. (a -> b) -> [a] -> [b]
map ((Vector -> Vector Int -> Vector) -> Vector Int -> Vector -> Vector
forall a b c. (a -> b -> c) -> b -> a -> c
flip Vector -> Vector Int -> Vector
forall a. Unbox a => Vector a -> Vector Int -> Vector a
U.backpermute Vector Int
ixs) [Vector]
preds0
           (Vector, Double) -> IO (Vector, Double)
forall (m :: * -> *) a. Monad m => a -> m a
return ((Vector, Double) -> IO (Vector, Double))
-> (Vector, Double) -> IO (Vector, Double)
forall a b. (a -> b) -> a -> b
$ [Vector] -> Vector -> (Vector, Double)
rgrss [Vector]
preds Vector
resp
      Vector (Vector, Double) -> ()
forall a. NFData a => a -> ()
rnf Vector (Vector, Double)
v () -> IO (Vector (Vector, Double)) -> IO (Vector (Vector, Double))
`seq` Vector (Vector, Double) -> IO (Vector (Vector, Double))
forall (m :: * -> *) a. Monad m => a -> m a
return Vector (Vector, Double)
v
  let (Vector Vector
coeffsv, Vector Double
r2v) = Vector (Vector, Double) -> (Vector Vector, Vector Double)
forall (v :: * -> *) a b.
(Vector v a, Vector v b, Vector v (a, b)) =>
v (a, b) -> (v a, v b)
G.unzip ([Vector (Vector, Double)] -> Vector (Vector, Double)
forall a. [Vector a] -> Vector a
V.concat [Vector (Vector, Double)]
vs)
  let coeffs :: Vector (Estimate ConfInt Double)
coeffs  = ((Int -> Double -> Estimate ConfInt Double)
 -> Vector Double -> Vector (Estimate ConfInt Double))
-> Vector Double
-> (Int -> Double -> Estimate ConfInt Double)
-> Vector (Estimate ConfInt Double)
forall a b c. (a -> b -> c) -> b -> a -> c
flip (Int -> Double -> Estimate ConfInt Double)
-> Vector Double -> Vector (Estimate ConfInt Double)
forall (v :: * -> *) a b.
(Vector v a, Vector v b) =>
(Int -> a -> b) -> v a -> v b
G.imap (Vector -> Vector Double
forall (v :: * -> *) a (w :: * -> *).
(Vector v a, Vector w a) =>
v a -> w a
G.convert Vector
coeffss) ((Int -> Double -> Estimate ConfInt Double)
 -> Vector (Estimate ConfInt Double))
-> (Int -> Double -> Estimate ConfInt Double)
-> Vector (Estimate ConfInt Double)
forall a b. (a -> b) -> a -> b
$ \Int
i Double
x ->
                Double -> Vector -> Estimate ConfInt Double
est Double
x (Vector -> Estimate ConfInt Double)
-> ((Int -> Double) -> Vector)
-> (Int -> Double)
-> Estimate ConfInt Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> (Int -> Double) -> Vector
forall a. Unbox a => Int -> (Int -> a) -> Vector a
U.generate Int
numResamples ((Int -> Double) -> Estimate ConfInt Double)
-> (Int -> Double) -> Estimate ConfInt Double
forall a b. (a -> b) -> a -> b
$ \Int
k -> (Vector Vector
coeffsv Vector Vector -> Int -> Vector
forall (v :: * -> *) a. Vector v a => v a -> Int -> a
G.! Int
k) Vector -> Int -> Double
forall (v :: * -> *) a. Vector v a => v a -> Int -> a
G.! Int
i
      r2 :: Estimate ConfInt Double
r2      = Double -> Vector -> Estimate ConfInt Double
est Double
r2s (Vector Double -> Vector
forall (v :: * -> *) a (w :: * -> *).
(Vector v a, Vector w a) =>
v a -> w a
G.convert Vector Double
r2v)
      (Vector
coeffss, Double
r2s) = [Vector] -> Vector -> (Vector, Double)
rgrss [Vector]
preds0 Vector
resp0
      est :: Double -> Vector -> Estimate ConfInt Double
est Double
s Vector
v = Double -> (Double, Double) -> CL Double -> Estimate ConfInt Double
forall a. Num a => a -> (a, a) -> CL Double -> Estimate ConfInt a
estimateFromInterval Double
s (Vector
w Vector -> Int -> Double
forall (v :: * -> *) a. Vector v a => v a -> Int -> a
G.! Int
lo, Vector
w Vector -> Int -> Double
forall (v :: * -> *) a. Vector v a => v a -> Int -> a
G.! Int
hi) CL Double
cl
        where w :: Vector
w  = Vector -> Vector
F.sort Vector
v
              bounded :: Int -> Int
bounded Int
i = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min (Vector -> Int
forall a. Unbox a => Vector a -> Int
U.length Vector
w Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
0 Int
i)
              lo :: Int
lo = Int -> Int
bounded (Int -> Int) -> Int -> Int
forall a b. (a -> b) -> a -> b
$ Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round Double
c
              hi :: Int
hi = Int -> Int
bounded (Int -> Int) -> Int -> Int
forall a b. (a -> b) -> a -> b
$ Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
truncate (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
c)
              n :: Double
n  = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
numResamples
              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)
  (Vector (Estimate ConfInt Double), Estimate ConfInt Double)
-> IO (Vector (Estimate ConfInt Double), Estimate ConfInt Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Vector (Estimate ConfInt Double)
coeffs, Estimate ConfInt Double
r2)

-- | Balance units of work across workers.
balance :: Int -> Int -> [Int]
balance :: Int -> Int -> [Int]
balance Int
numSlices Int
numItems = (Int -> Int -> Int) -> [Int] -> [Int] -> [Int]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Int -> Int -> Int
forall a. Num a => a -> a -> a
(+) (Int -> Int -> [Int]
forall a. Int -> a -> [a]
replicate Int
numSlices Int
q)
                                         (Int -> Int -> [Int]
forall a. Int -> a -> [a]
replicate Int
r Int
1 [Int] -> [Int] -> [Int]
forall a. [a] -> [a] -> [a]
++ Int -> [Int]
forall a. a -> [a]
repeat Int
0)
 where (Int
q,Int
r) = Int
numItems Int -> Int -> (Int, Int)
forall a. Integral a => a -> a -> (a, a)
`quotRem` Int
numSlices