{-# LANGUAGE CPP #-}
{-
	Copyright (C) 2011-2017 Dr. Alistair Ward

	This program is free software: you can redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation, either version 3 of the License, or
	(at your option) any later version.

	This program is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with this program.  If not, see <http://www.gnu.org/licenses/>.
-}
{- |
 [@AUTHOR@]	Dr. Alistair Ward

 [@DESCRIPTION@]	Miscellaneous statistics functions.
-}

module Factory.Math.Statistics(
-- * Functions
	getMean,
	getRootMeanSquare,
	getWeightedMean,
--	getDispersionFromMean,
	getVariance,
	getStandardDeviation,
	getAverageAbsoluteDeviation,
	getCoefficientOfVariance,
	nCr,
	nPr
) where

import			Control.Arrow((***))
import qualified	Control.Exception
import			Control.Parallel(par, pseq)
import qualified	Data.Foldable
import qualified	Data.List
import qualified	Factory.Math.Factorial			as Math.Factorial
import qualified	Factory.Math.Implementations.Factorial	as Math.Implementations.Factorial
import qualified	Factory.Math.Power			as Math.Power

{- |
	* Determines the /mean/ of the specified numbers; <https://en.wikipedia.org/wiki/Mean>.

	* Should the caller define the result-type as 'Rational', then it will be free from rounding-errors.
-}
getMean :: (
	Data.Foldable.Foldable	foldable,
	Fractional		result,
	Real			value
 )
	=> foldable value
	-> result
{-# SPECIALISE getMean :: Data.Foldable.Foldable foldable => foldable Double -> Double #-}
#if MIN_TOOL_VERSION_ghc(7,10,0)
{-# SPECIALISE getMean :: Data.Foldable.Foldable foldable => foldable Rational -> Rational #-}
#endif
getMean :: foldable value -> result
getMean foldable value
foldable	= Bool -> result -> result
forall a. (?callStack::CallStack) => Bool -> a -> a
Control.Exception.assert (Int
denominator Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
0) (result -> result) -> result -> result
forall a b. (a -> b) -> a -> b
$ value -> result
forall a b. (Real a, Fractional b) => a -> b
realToFrac value
numerator result -> result -> result
forall a. Fractional a => a -> a -> a
/ Int -> result
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
denominator	where
	denominator :: Int
	(value
numerator, Int
denominator)	= ((value, Int) -> value -> (value, Int))
-> (value, Int) -> foldable value -> (value, Int)
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
Data.Foldable.foldl' (
		\(value, Int)
acc value
x	-> let
			acc' :: (value, Int)
acc'@(value
n, Int
d)	= (value -> value -> value
forall a. Num a => a -> a -> a
+ value
x) (value -> value) -> (Int -> Int) -> (value, Int) -> (value, Int)
forall (a :: * -> * -> *) b c b' c'.
Arrow a =>
a b c -> a b' c' -> a (b, b') (c, c')
*** Int -> Int
forall a. Enum a => a -> a
succ ((value, Int) -> (value, Int)) -> (value, Int) -> (value, Int)
forall a b. (a -> b) -> a -> b
$ (value, Int)
acc
		in value
n value -> (value, Int) -> (value, Int)
`seq` Int
d Int -> (value, Int) -> (value, Int)
`seq` (value, Int)
acc'
	 ) (value
0, Int
0) foldable value
foldable

-- | Determines the /root mean square/ of the specified numbers; <https://en.wikipedia.org/wiki/Root_mean_square>.
getRootMeanSquare :: (
	Data.Foldable.Foldable	foldable,
	Floating		result,
	Real			value
 )
	=> foldable value
	-> result
{-# SPECIALISE getRootMeanSquare :: Data.Foldable.Foldable foldable => foldable Double -> Double #-}
getRootMeanSquare :: foldable value -> result
getRootMeanSquare foldable value
foldable	= Bool -> result -> result
forall a. (?callStack::CallStack) => Bool -> a -> a
Control.Exception.assert (Int
denominator Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
0) (result -> result) -> (result -> result) -> result -> result
forall b c a. (b -> c) -> (a -> b) -> a -> c
. result -> result
forall a. Floating a => a -> a
sqrt (result -> result) -> result -> result
forall a b. (a -> b) -> a -> b
$ value -> result
forall a b. (Real a, Fractional b) => a -> b
realToFrac value
numerator result -> result -> result
forall a. Fractional a => a -> a -> a
/ Int -> result
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
denominator	where
	denominator :: Int
	(value
numerator, Int
denominator)	= ((value, Int) -> value -> (value, Int))
-> (value, Int) -> foldable value -> (value, Int)
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
Data.Foldable.foldl' (
		\(value, Int)
acc value
x -> let
			acc' :: (value, Int)
acc'@(value
n, Int
d)	= (value -> value -> value
forall a. Num a => a -> a -> a
+ value -> value
forall n. Num n => n -> n
Math.Power.square value
x) (value -> value) -> (Int -> Int) -> (value, Int) -> (value, Int)
forall (a :: * -> * -> *) b c b' c'.
Arrow a =>
a b c -> a b' c' -> a (b, b') (c, c')
*** Int -> Int
forall a. Enum a => a -> a
succ ((value, Int) -> (value, Int)) -> (value, Int) -> (value, Int)
forall a b. (a -> b) -> a -> b
$ (value, Int)
acc
		 in value
n value -> (value, Int) -> (value, Int)
`seq` Int
d Int -> (value, Int) -> (value, Int)
`seq` (value, Int)
acc'
	 ) (value
0, Int
0) foldable value
foldable

{- |
	* Determines the /weighted mean/ of the specified numbers; <https://en.wikipedia.org/wiki/Weighted_arithmetic_mean>.

	* The specified value is only evaluated if the corresponding weight is non-zero.

	* CAVEAT: because the operand is more general than a list, no optimisation is performed when supplied a singleton.
-}
getWeightedMean :: (
	Data.Foldable.Foldable	foldable,
	Fractional		result,
	Real			value,
	Real			weight
 )
	=> foldable (value, weight)	-- ^ Each pair consists of a value & the corresponding weight.
	-> result
{-# SPECIALISE getWeightedMean :: Data.Foldable.Foldable foldable => foldable (Double, Double) -> Double #-}
#if MIN_TOOL_VERSION_ghc(7,10,0)
{-# SPECIALISE getWeightedMean :: Data.Foldable.Foldable foldable => foldable (Rational, Rational) -> Rational #-}
#endif
getWeightedMean :: foldable (value, weight) -> result
getWeightedMean foldable (value, weight)
foldable	= Bool -> result -> result
forall a. (?callStack::CallStack) => Bool -> a -> a
Control.Exception.assert (weight
denominator weight -> weight -> Bool
forall a. Eq a => a -> a -> Bool
/= weight
0) (result -> result) -> result -> result
forall a b. (a -> b) -> a -> b
$ result
numerator result -> result -> result
forall a. Fractional a => a -> a -> a
/ weight -> result
forall a b. (Real a, Fractional b) => a -> b
realToFrac weight
denominator	where
	(result
numerator, weight
denominator)	= ((result, weight) -> (value, weight) -> (result, weight))
-> (result, weight) -> foldable (value, weight) -> (result, weight)
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
Data.Foldable.foldl' (
		\(result, weight)
acc (value
value, weight
weight)	-> if weight
weight weight -> weight -> Bool
forall a. Eq a => a -> a -> Bool
== weight
0
			then (result, weight)
acc	-- Avoid unnecessary evaluation.
			else let
				acc' :: (result, weight)
acc'@(result
n, weight
d)	= (result -> result -> result
forall a. Num a => a -> a -> a
+ weight -> result
forall a b. (Real a, Fractional b) => a -> b
realToFrac weight
weight result -> result -> result
forall a. Num a => a -> a -> a
* value -> result
forall a b. (Real a, Fractional b) => a -> b
realToFrac value
value) (result -> result)
-> (weight -> weight) -> (result, weight) -> (result, weight)
forall (a :: * -> * -> *) b c b' c'.
Arrow a =>
a b c -> a b' c' -> a (b, b') (c, c')
*** (weight -> weight -> weight
forall a. Num a => a -> a -> a
+ weight
weight) ((result, weight) -> (result, weight))
-> (result, weight) -> (result, weight)
forall a b. (a -> b) -> a -> b
$ (result, weight)
acc
			in result
n result -> (result, weight) -> (result, weight)
`seq` weight
d weight -> (result, weight) -> (result, weight)
`seq` (result, weight)
acc'
	 ) (result
0, weight
0) foldable (value, weight)
foldable

{- |
	* Measures the /dispersion/ of a /population/ of results from the /mean/ value; <https://en.wikipedia.org/wiki/Statistical_dispersion>.

	* Should the caller define the result-type as 'Rational', then it will be free from rounding-errors.
-}
getDispersionFromMean :: (
	Data.Foldable.Foldable	foldable,
	Fractional		result,
	Functor			foldable,
	Real			value
 ) => (Rational -> Rational) -> foldable value -> result
{-# SPECIALISE getDispersionFromMean :: (Data.Foldable.Foldable foldable, Functor foldable) => (Rational -> Rational) -> foldable Double -> Double #-}
#if MIN_TOOL_VERSION_ghc(7,10,0)
{-# SPECIALISE getDispersionFromMean :: (Data.Foldable.Foldable foldable, Functor foldable) => (Rational -> Rational) -> foldable Rational -> Rational #-}
#endif
getDispersionFromMean :: (Rational -> Rational) -> foldable value -> result
getDispersionFromMean Rational -> Rational
weight foldable value
foldable	= foldable Rational -> result
forall (foldable :: * -> *) result value.
(Foldable foldable, Fractional result, Real value) =>
foldable value -> result
getMean (foldable Rational -> result) -> foldable Rational -> result
forall a b. (a -> b) -> a -> b
$ (value -> Rational) -> foldable value -> foldable Rational
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (Rational -> Rational
weight (Rational -> Rational) -> (value -> Rational) -> value -> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
subtract Rational
mean (Rational -> Rational) -> (value -> Rational) -> value -> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. value -> Rational
forall a. Real a => a -> Rational
toRational) foldable value
foldable	where
	mean :: Rational
mean	= foldable value -> Rational
forall (foldable :: * -> *) result value.
(Foldable foldable, Fractional result, Real value) =>
foldable value -> result
getMean foldable value
foldable

{- |
	* Determines the exact /variance/ of the specified numbers; <https://en.wikipedia.org/wiki/Variance>.

	* Should the caller define the result-type as 'Rational', then it will be free from rounding-errors.
-}
getVariance :: (
	Data.Foldable.Foldable	foldable,
	Fractional		variance,
	Functor			foldable,
	Real			value
 ) => foldable value -> variance
{-# SPECIALISE getVariance :: (Data.Foldable.Foldable foldable, Functor foldable) => foldable Double -> Double #-}
#if MIN_TOOL_VERSION_ghc(7,10,0)
{-# SPECIALISE getVariance :: (Data.Foldable.Foldable foldable, Functor foldable) => foldable Rational -> Rational #-}
#endif
getVariance :: foldable value -> variance
getVariance	= (Rational -> Rational) -> foldable value -> variance
forall (foldable :: * -> *) result value.
(Foldable foldable, Fractional result, Functor foldable,
 Real value) =>
(Rational -> Rational) -> foldable value -> result
getDispersionFromMean Rational -> Rational
forall n. Num n => n -> n
Math.Power.square

-- | Determines the /standard-deviation/ of the specified numbers; <https://en.wikipedia.org/wiki/Standard_deviation>.
getStandardDeviation :: (
	Data.Foldable.Foldable	foldable,
	Floating		result,
	Functor			foldable,
	Real			value
 ) => foldable value -> result
{-# SPECIALISE getStandardDeviation :: (Data.Foldable.Foldable foldable, Functor foldable) => foldable Double -> Double #-}
getStandardDeviation :: foldable value -> result
getStandardDeviation	= result -> result
forall a. Floating a => a -> a
sqrt (result -> result)
-> (foldable value -> result) -> foldable value -> result
forall b c a. (b -> c) -> (a -> b) -> a -> c
. foldable value -> result
forall (foldable :: * -> *) variance value.
(Foldable foldable, Fractional variance, Functor foldable,
 Real value) =>
foldable value -> variance
getVariance

{- |
	* Determines the /average absolute deviation/ of the specified numbers; <https://en.wikipedia.org/wiki/Absolute_deviation#Average_absolute_deviation>.

	* Should the caller define the result-type as 'Rational', then it will be free from rounding-errors.
-}
getAverageAbsoluteDeviation :: (
	Data.Foldable.Foldable	foldable,
	Fractional		result,
	Functor			foldable,
	Real			value
 ) => foldable value -> result
{-# SPECIALISE getAverageAbsoluteDeviation :: (Data.Foldable.Foldable foldable, Functor foldable) => foldable Double -> Double #-}
#if MIN_TOOL_VERSION_ghc(7,10,0)
{-# SPECIALISE getAverageAbsoluteDeviation :: (Data.Foldable.Foldable foldable, Functor foldable) => foldable Rational -> Rational #-}
#endif
getAverageAbsoluteDeviation :: foldable value -> result
getAverageAbsoluteDeviation	= (Rational -> Rational) -> foldable value -> result
forall (foldable :: * -> *) result value.
(Foldable foldable, Fractional result, Functor foldable,
 Real value) =>
(Rational -> Rational) -> foldable value -> result
getDispersionFromMean Rational -> Rational
forall n. Num n => n -> n
abs

-- | Determines the /coefficient-of-variance/ of the specified numbers; <https://en.wikipedia.org/wiki/Coefficient_of_variation>.
getCoefficientOfVariance :: (
	Data.Foldable.Foldable	foldable,
	Eq			result,
	Floating		result,
	Functor			foldable,
	Real			value
 ) => foldable value -> result
{-# SPECIALISE getCoefficientOfVariance :: (Data.Foldable.Foldable foldable, Functor foldable) => foldable Double -> Double #-}
getCoefficientOfVariance :: foldable value -> result
getCoefficientOfVariance foldable value
l	= Bool -> result -> result
forall a. (?callStack::CallStack) => Bool -> a -> a
Control.Exception.assert (result
mean result -> result -> Bool
forall a. Eq a => a -> a -> Bool
/= result
0) (result -> result) -> result -> result
forall a b. (a -> b) -> a -> b
$ foldable value -> result
forall (foldable :: * -> *) result value.
(Foldable foldable, Floating result, Functor foldable,
 Real value) =>
foldable value -> result
getStandardDeviation foldable value
l result -> result -> result
forall a. Fractional a => a -> a -> a
/ result -> result
forall n. Num n => n -> n
abs result
mean	where
	mean :: result
mean	= foldable value -> result
forall (foldable :: * -> *) result value.
(Foldable foldable, Fractional result, Real value) =>
foldable value -> result
getMean foldable value
l

-- | The number of unordered /combinations/ of /r/ objects taken from /n/; <https://en.wikipedia.org/wiki/Combination>.
nCr :: (Math.Factorial.Algorithmic factorialAlgorithm, Integral i, Show i)
	=> factorialAlgorithm
	-> i	-- ^ The total number of items from which to select.
	-> i	-- ^ The number of items in a sample.
	-> i	-- ^ The number of combinations.
nCr :: factorialAlgorithm -> i -> i -> i
nCr factorialAlgorithm
_ i
0 i
_	= i
1
nCr factorialAlgorithm
_ i
_ i
0	= i
1
nCr factorialAlgorithm
factorialAlgorithm i
n i
r
	| i
n i -> i -> Bool
forall a. Ord a => a -> a -> Bool
< i
r		= i
0
	| Bool
otherwise	= Bool -> i -> i
forall a. (?callStack::CallStack) => Bool -> a -> a
Control.Exception.assert (i
n i -> i -> Bool
forall a. Ord a => a -> a -> Bool
>= i
0 Bool -> Bool -> Bool
&& i
r i -> i -> Bool
forall a. Ord a => a -> a -> Bool
>= i
0) (i -> i) -> i -> i
forall a b. (a -> b) -> a -> b
$ i
numerator i -> i -> i
forall a b. a -> b -> b
`par` (i
denominator i -> i -> i
forall a b. a -> b -> b
`pseq` i
numerator i -> i -> i
forall a. Integral a => a -> a -> a
`div` i
denominator)
	where
		[i
smaller, i
bigger]	= [i] -> [i]
forall a. Ord a => [a] -> [a]
Data.List.sort [i
r, i
n i -> i -> i
forall a. Num a => a -> a -> a
- i
r]
		numerator :: i
numerator		= i -> i -> i
forall i. (Integral i, Show i) => i -> i -> i
Math.Implementations.Factorial.risingFactorial (i -> i
forall a. Enum a => a -> a
succ i
bigger) (i
n i -> i -> i
forall a. Num a => a -> a -> a
- i
bigger)
		denominator :: i
denominator		= factorialAlgorithm -> i -> i
forall algorithm i.
(Algorithmic algorithm, Integral i, Show i) =>
algorithm -> i -> i
Math.Factorial.factorial factorialAlgorithm
factorialAlgorithm i
smaller

-- | The number of /permutations/ of /r/ objects taken from /n/; <https://en.wikipedia.org/wiki/Permutations>.
nPr :: (Integral i, Show i)
	=> i	-- ^ The total number of items from which to select.
	-> i	-- ^ The number of items in a sample.
	-> i	-- ^ The number of permutations.
nPr :: i -> i -> i
nPr i
0 i
_	= i
1
nPr i
_ i
0	= i
1
nPr i
n i
r
	| i
n i -> i -> Bool
forall a. Ord a => a -> a -> Bool
< i
r		= i
0
	| Bool
otherwise	= Bool -> i -> i
forall a. (?callStack::CallStack) => Bool -> a -> a
Control.Exception.assert (i
n i -> i -> Bool
forall a. Ord a => a -> a -> Bool
>= i
0 Bool -> Bool -> Bool
&& i
r i -> i -> Bool
forall a. Ord a => a -> a -> Bool
>= i
0) (i -> i) -> i -> i
forall a b. (a -> b) -> a -> b
$ i -> i -> i
forall i. (Integral i, Show i) => i -> i -> i
Math.Implementations.Factorial.fallingFactorial i
n i
r