numeric-prelude-0.2.2: test/Test/MathObj/Gaussian/Variance.hs
{-# LANGUAGE RebindableSyntax #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances #-}
module Test.MathObj.Gaussian.Variance where
import qualified MathObj.Gaussian.Variance as G
import qualified Number.Root as Root
-- import qualified Algebra.Ring as Ring
import qualified Algebra.Laws as Laws
import Test.NumericPrelude.Utility (testUnit)
import Test.QuickCheck (Testable, quickCheck, (==>), Arbitrary, arbitrary, )
import qualified Test.HUnit as HUnit
import Control.Monad (liftM2, liftM3, )
import Data.Function.HT (nest, compose2, )
import NumericPrelude.Base as P
import NumericPrelude.Numeric as NP
newtype PositiveInteger = PositiveInteger Integer
deriving Show
instance Arbitrary PositiveInteger where
arbitrary =
fmap (\p -> PositiveInteger $ 1 + abs p) arbitrary
{- |
For @(HoelderConjugates p q)@ it holds
> 1/p + 1/q = 1
-}
data HoelderConjugates = HoelderConjugates Rational Rational
deriving Show
{-
instance Arbitrary HoelderConjugates where
arbitrary = liftM2
(\(PositiveInteger p) (PositiveInteger q) ->
let s = 1%p + 1%q
in HoelderConjugates (fromInteger p * s) (fromInteger q * s))
arbitrary arbitrary
-}
instance Arbitrary HoelderConjugates where
arbitrary = liftM2
(\(PositiveInteger p) (PositiveInteger q) ->
let s = p + q
in HoelderConjugates (s % p) (s % q))
arbitrary arbitrary
{- |
For @(YoungConjugates p q r)@ it holds
> 1/p + 1/q = 1/r + 1
-}
data YoungConjugates = YoungConjugates Rational Rational Rational
deriving Show
{-
Find positive natural numbers @a, b, c, d@ with
> a + b = c + d
and
> d >= a, d >= b, d >= c
then set
> p=d/a, q=d/b, r=d/c
a+b<=c
b+c<=a
-> 2b <= 0
-}
instance Arbitrary YoungConjugates where
arbitrary = liftM3
(\(PositiveInteger a0) (PositiveInteger b0) (PositiveInteger c0) ->
let guardSwap cond (x,y) =
if cond x y then (x,y) else (y,x)
{-
If a+b<=c, then from b>0 it follows a<c and thus c+b>a.
Swapping a and c is enough and we have not to consider more cases.
-}
(a1,c1) = guardSwap (\a c -> a+b0>c) (a0,c0)
b1 = b0
d1 = a1+b1-c1
((a2,b2),(c2,d2)) =
guardSwap (compose2 (<=) snd)
(guardSwap (<=) (a1,b1),
guardSwap (<=) (c1,d1))
in YoungConjugates (d2%a2) (d2%b2) (d2%c2))
arbitrary arbitrary arbitrary
simple ::
(Testable t) =>
(G.T Rational -> t) -> IO ()
simple f =
quickCheck (\x -> f (x :: G.T Rational))
tests :: HUnit.Test
tests =
HUnit.TestLabel "variance" $
HUnit.TestList $
map testUnit $
testList
testList :: [(String, IO ())]
testList =
{-
("convolution, dirac",
simple $ Laws.identity (+) zero) :
-}
("convolution, commutative",
simple $ Laws.commutative G.convolve) :
("convolution, associative",
simple $ Laws.associative G.convolve) :
("multiplication, one",
simple $ Laws.identity G.multiply G.constant) :
("multiplication, commutative",
simple $ Laws.commutative G.multiply) :
("multiplication, associative",
simple $ Laws.associative G.multiply) :
("convolution via fourier",
simple $ \x y ->
G.fourier (G.convolve x y)
== G.multiply (G.fourier x) (G.fourier y)) :
("fourier identity",
simple $ \x -> nest 4 G.fourier x == x) :
("dilate multiplicative",
simple $ \x a b -> a>0 && b>0 ==>
G.dilate a (G.dilate b x) == G.dilate (a*b) x) :
("dilate by shrink",
simple $ \x a -> a>0 ==>
G.shrink a x == G.dilate (recip a) x) :
("fourier dilate",
simple $ \x a -> a>0 ==>
G.fourier (G.dilate a x) == G.amplify a (G.shrink a (G.fourier x))) :
("fourier, unitary",
simple $ \x y ->
G.scalarProductRoot x y
== G.scalarProductRoot (G.fourier x) (G.fourier y)) :
("norm1 vs. normP 1",
simple $ \x -> G.norm1Root x == G.normPRoot 1 x) :
("norm2 vs. normP 2",
simple $ \x -> G.norm2Root x == G.normPRoot 2 x) :
{-
I would have liked to test for a monotony of norms.
Unfortunately, it does not hold.
Means contain a division by the size of the domain.
Norms do not have this division.
Means are monotonic with respect to the degree.
Norms are not.
We cannot turn the norms into means since the size of the domain
(the complete real axis) is infinitely large.
("norm monotony",
simple $ \x p0 q0 ->
let p = 1 + abs p0
q = 1 + abs q0
in case compare p q of
EQ -> G.normPRoot p x == G.normPRoot q x
LT -> G.normPRoot p x <= G.normPRoot q x
GT -> G.normPRoot p x >= G.normPRoot q x) :
This should also fail,
but QuickCheck does not seem to try counterexamples.
("infinity norm upper bound",
simple $ \x p0 ->
let p = 1 + abs p0
in G.normPRoot p x <= G.normInfRoot x) :
-}
("Cauchy-Schwarz inequality",
simple $ \x y ->
G.scalarProductRoot x y
<= G.norm2Root x `Root.mul` G.norm2Root y) :
("Hoelder conjugates",
quickCheck $ \(HoelderConjugates p q) ->
p>=1 && q>=1 && 1/p + 1/q == 1) :
("Hoelder inequality with infinity norm",
simple $ \x y ->
G.scalarProductRoot x y
<= G.norm1Root x `Root.mul` G.normInfRoot y) :
("Hoelder inequality",
simple $ \x y (HoelderConjugates p q) ->
G.scalarProductRoot x y
<= G.normPRoot p x `Root.mul` G.normPRoot q y) :
("Young inequality with two infinity norms",
simple $ \x y ->
G.normInfRoot (G.convolve x y)
<= G.norm1Root x `Root.mul` G.normInfRoot y) :
("Young inequality with infinity norm",
simple $ \x y (HoelderConjugates p q) ->
G.normInfRoot (G.convolve x y)
<= G.normPRoot p x `Root.mul` G.normPRoot q y) :
("Young conjugates",
quickCheck $ \(YoungConjugates p q r) ->
p>=1 && q>=1 && r>=1 && 1/p + 1/q == 1/r + 1) :
("Young inequality",
simple $ \x y (YoungConjugates p q r) ->
G.normPRoot r (G.convolve x y)
<= G.normPRoot p x `Root.mul` G.normPRoot q y) :
[]