numeric-prelude-0.1.1: src/MathObj/Gaussian/Example.hs
{-# LANGUAGE NoImplicitPrelude #-}
{-
Reciprocal of variance of a Gaussian bell curve.
We describe the curve only in terms of its variance
thus we represent a bell curve at the coordinate origin
neglecting its amplitude.
We could also define the amplitude as @root 4 c@,
thus preserving L2 norm being one,
but then @dilate@ and @shrink@ also include an amplification.
We could do some projective geometry in the exponent
in order to also have zero variance,
which corresponds to the dirac impulse.
-}
module MathObj.Gaussian.Example where
import qualified MathObj.Gaussian.Polynomial as PolyBell
import qualified MathObj.Gaussian.Bell as Bell
import qualified MathObj.Gaussian.Variance as Var
import qualified MathObj.Polynomial as Poly
import qualified Algebra.Transcendental as Trans
import qualified Algebra.Algebraic as Algebraic
import qualified Algebra.Field as Field
-- import qualified Algebra.Real as Real
import qualified Algebra.Ring as Ring
-- import qualified Algebra.Additive as Additive
import qualified Number.Complex as Complex
import Algebra.Transcendental (pi, )
import Algebra.Algebraic (root, )
import Algebra.Ring ((*), (^), )
import Number.Complex ((+:), )
import qualified Numerics.Function as Func
import qualified Numerics.Fourier as Fourier
import qualified Numerics.Integration as Integ
import qualified Numerics.Differentiation as Diff
import qualified Graphics.Gnuplot.Simple as GP
import Control.Applicative (liftA2, )
-- import System.Exit (ExitCode, )
-- import Prelude (($))
import NumericPrelude
import PreludeBase
import qualified Prelude as P
curve0 :: Var.T Double
curve0 = curve0a
curve0a :: Var.T Double
curve0a = Var.Cons 1.4 3.3
curve0b :: Var.T Double
curve0b = Var.Cons 2.2 1.7
variance0 :: (Double, Double)
variance0 =
(Var.variance curve0,
(Integ.rectangular 1000 (-2,2) $ liftA2 (*) (^2) (Var.evaluate curve0)) /
(Integ.rectangular 1000 (-2,2) $ Var.evaluate curve0))
norm10 :: (Double, Double)
norm10 =
(Integ.rectangular 1000 (-2,2) $ Var.evaluate curve0,
Var.norm1 curve0)
norm20 :: (Double, Double)
norm20 =
(sqrt $ Integ.rectangular 1000 (-2,2) $ (^2) . Var.evaluate curve0,
Var.norm2 curve0)
norm30 :: (Double, Double)
norm30 =
(root 3 $ Integ.rectangular 1000 (-2,2) $ (^3) . Var.evaluate curve0,
Var.normP 3 curve0)
fourier0 :: IO ()
fourier0 =
GP.plotFuncs []
(GP.linearScale 100 (-2,2))
[Var.evaluate $ Var.fourier curve0,
Fourier.analysisTransformOneReal 100 (-2,2) $ Var.evaluate curve0]
multiply0 :: IO ()
multiply0 =
GP.plotFuncs []
(GP.linearScale 100 (-1,1))
[Var.evaluate $ Var.multiply curve0a curve0b,
liftA2 (*) (Var.evaluate curve0a) (Var.evaluate curve0b)]
convolve0 :: IO ()
convolve0 =
GP.plotFuncs []
(GP.linearScale 100 (-2,2))
[Var.evaluate $ Var.convolve curve0a curve0b,
Integ.convolve 1000 (-3,3) (Var.evaluate curve0a) (Var.evaluate curve0b)]
curve1 :: Bell.T Double
curve1 = curve1a
curve1a :: Bell.T Double
curve1a = Bell.Cons 1.4 (0.1+:0.3) ((-0.2)+:1.4) 2.3
curve1b :: Bell.T Double
curve1b = Bell.Cons 2.2 ((-0.3)+:2.1) (0.2+:(-0.4)) 1.7
variance1 :: (Double, Double)
variance1 =
(Bell.variance curve1,
(Integ.rectangular 1000 (-2,2) $
liftA2 (*) (^2)
(Complex.magnitudeSqr .
Func.translateRight
(Complex.real (Bell.c1 curve1) / (2 * Bell.c2 curve1))
(Bell.evaluate curve1))) /
(Integ.rectangular 1000 (-2,2) $ Complex.magnitude . Bell.evaluate curve1))
{- the norm depends on too much things
norm0vs1 :: (Double, Double)
norm0vs1 =
((Integ.rectangular 1000 (-5,5) $ Var.evaluate curve0)
* exp (- Complex.real (Bell.c0 curve1)),
Integ.rectangular 1000 (-5,5) $ Complex.magnitude . Bell.evaluate curve1)
-}
fourier1 :: IO ()
fourier1 =
GP.plotFuncs []
(GP.linearScale 100 (-5,5))
[Complex.real . (Bell.evaluate $ Bell.fourier curve1),
fourierAnalysisReal 100 (-2,2) $ Bell.evaluate curve1]
curve2 :: PolyBell.T Double
curve2 =
PolyBell.Cons
-- Bell.unit
-- (Bell.Cons 1.4 (0.1+:0.3) 0 1.2)
-- (Bell.Cons 1.4 (0.1+:0.3) ((-0.2)+:1.4) 1)
curve1
-- (Poly.fromCoeffs [one])
-- (Poly.fromCoeffs [zero,one])
-- (Poly.fromCoeffs [zero,zero,one])
-- (Poly.fromCoeffs [0,Complex.imaginaryUnit])
(Poly.fromCoeffs [1.4+:(-0.1),0.8+:(0.1),(-1.1)+:0.3])
differentiate2 :: IO ()
differentiate2 =
GP.plotFuncs []
(GP.linearScale 100 (-2,2))
[Complex.real . (PolyBell.evaluateSqRt $ PolyBell.differentiate curve2),
((/ sqrt pi) . ) $ Diff.diff (1e-5) $ Complex.real . PolyBell.evaluateSqRt curve2]
fourier2 :: IO ()
fourier2 =
GP.plotFuncs []
(GP.linearScale 100 (-5,5))
[Complex.real . (PolyBell.evaluateSqRt $ PolyBell.fourier curve2),
fourierAnalysisReal 100 (-2,2) $ PolyBell.evaluateSqRt curve2]
fourierAnalysisReal ::
(P.Floating a) =>
Integer -> (a, a) -> (a -> Complex.T a) -> a -> a
fourierAnalysisReal n rng f =
liftA2 (P.-)
(Fourier.analysisTransformOneReal n rng (Complex.real . f))
(Fourier.analysisTransformOneImag n rng (Complex.imag . f))
{- |
Try to approximate @\x -> exp (-x^2) * x@
by a difference of translated Gaussian bells.
exp(-x^2) * x
== exp(-(a+b*x+c*x^2)) - exp(-(a-b*x+c*x^2))
== exp(-(a+c*x^2)) * (exp(-b*x) - exp(b*x))
== exp(-(a+c*x^2)) * 2*sinh (b*x)
It holds
lim (\b x -> sinh (b*x) / b) = id
-}
diffApprox :: IO ()
diffApprox =
let amp = (2*b)^- (-2)
a = 0
{-
amp = 1
a = log (2 * abs b)
-}
b = -0.1
c = 1
ac = Complex.fromReal a
bc = Complex.fromReal b
in GP.plotFuncs []
(GP.linearScale 100 (-2,2::Double))
[Complex.real .
(PolyBell.evaluateSqRt $
PolyBell.Cons Bell.unit (Poly.fromCoeffs [zero,one])),
Complex.real .
liftA2 (-)
(PolyBell.evaluateSqRt $
PolyBell.Cons (Bell.Cons amp ac bc c) (Poly.fromCoeffs [one]))
(PolyBell.evaluateSqRt $
PolyBell.Cons (Bell.Cons amp ac (-bc) c) (Poly.fromCoeffs [one]))]
polyApprox :: IO ()
polyApprox =
GP.plotFuncs []
(GP.linearScale 100 (-2,2::Double))
[Complex.real .
PolyBell.evaluateSqRt curve2,
Complex.real . sum .
mapM (\(amp,b) -> \x -> amp * Bell.evaluateSqRt b x)
(PolyBell.approximateByBells 0.1 curve2)]