{-# LANGUAGE FlexibleContexts #-}
module Statistics.Gcm
(
gcm
, unsafeGcm
) where
import qualified Data.Vector.Generic as V
import Data.Vector.Generic (Vector)
import Statistics.Pava.Common
pool :: (Ord a, Real a, Show a, Real b, Show b)
=> [a] -> [b] -> [Double] -> a -> b -> ([a], [b], [Double])
pool xs@(x1:_) ys@(y1:_) ss@(s01:_) x2 y2 =
if s01 < s12
then (x2:xs, y2:ys, s12:ss)
else pool (tail xs) (tail ys) (tail ss) x2 y2
where
s12 = slope x1 x2 y1 y2
pool [x] [y] [] x' y' = (x':[x], y':[y], [slope x x' y y'])
pool xs ys ss x y = error $ "pool: xs, ys, ss, x, y: "
++ show xs ++ ", " ++ show ys ++ ", " ++ show ss
++ ", " ++ show x ++ ", " ++ show y ++ "."
{-# SPECIALIZE pool :: [Int] -> [Double] -> [Double] -> Int -> Double
-> ([Int], [Double], [Double]) #-}
{-# SPECIALIZE pool :: [Double] -> [Double] -> [Double] -> Double -> Double
-> ([Double], [Double], [Double]) #-}
gcm :: (Real a, Real b, Show a, Vector v a, Show b, Vector v b, Vector v Bool)
=> v a -> v b -> ([a], [b], [Double])
gcm ps rs | lPs /= lRs = error $ "gcm: Number of predictors is " ++ show lPs
++ ", but number of responses is " ++ show lRs ++ "."
| not (strictlyOrdered ps) = error "gcm: The predictors are not strictly ordered."
| otherwise = unsafeGcm ps rs
where lPs = V.length ps
lRs = V.length rs
unsafeGcm :: (Real a, Real b, Show a, Vector v a, Show b, Vector v b)
=> v a -> v b -> ([a], [b], [Double])
unsafeGcm ps rs | l == 0 = ([], [], [])
| l == 1 = start
| otherwise = reverse3 $ go start (1 :: Int)
where
l = V.length rs
start = ([V.head ps], [V.head rs], [])
go (xs, ys, ss) i | i >= l = (xs, ys, ss)
| otherwise = go (pool xs ys ss x' y') (i+1)
where
x' = ps V.! i
y' = rs V.! i