In [1]:
import Graphics.Rendering.Chart
import Data.Colour
import Data.Colour.Names
import Data.Default.Class
import Control.Lens
import qualified Data.Vector.Unboxed as V
import qualified Data.Number.Erf as E
setLinesBlue :: PlotLines a b -> PlotLines a b
setLinesBlue = plot_lines_style . line_color .~ opaque blue
am :: Double -> Double
am x = (sin (x*3.14159/45) + 1) / 2 * (sin (x*3.14159/5))
chart points = toRenderable layout
where
line = plot_lines_values .~ points
$ plot_lines_style . line_color .~ opaque blue
$ def
dots = plot_points_style .~ filledCircles 1 (opaque red)
$ plot_points_values .~ (concat points)
$ def
layout =
layout_plots .~ [toPlot line,
toPlot dots]
$ def
chart [[(x, am x) | x <- [0, 1..400]]]
In [2]:
-- | Faculty, we stay on the save side here and use Integer for big numbers. This is no performance code.
fac :: Integer -> Integer
fac 0 = 1
fac n = n * fac (n-1)
doubleFac :: Integer -> Integer
doubleFac 0 = 1
doubleFac n = product [n - 2*i | i <- [0 .. (n + 1) `div` 2 - 1]]
poissonPdf lambda k
| lambda > 0 = lambda ** fromIntegral k * exp (-lambda) / (fromIntegral $ fac k)
| otherwise = 0
In [3]:
-- Check that the distribution somewhat adds up to 1.
sum [poissonPdf 0.01 k | k <- [0, 1 .. 12]]
In [4]:
-- | Define some crude approximation to the erf function.
-- Taken from http://en.wikipedia.org/wiki/Error_function
erfDropIn1 :: (Floating d, Ord d) => d -> d
erfDropIn1 x
| x < 0 = -erfDropIn1 (-x)
| otherwise = 1 - 1 / (1 + a1*x + a2*x**2 + a3*x**3 + a4*x**4)**4
where
a1 = 0.278393
a2 = 0.230389
a3 = 0.000972
a4 = 0.078108
erfDropIn2 :: (Floating d, Ord d) => d -> d
erfDropIn2 x
| x < 0 = -erfDropIn2 (-x)
| otherwise = 1 - 1 / (1 + a1*x + a2*x**2 + a3*x**3 + a4*x**4 + a5*x**5 + a6*x**6)**16
where
a1 = 0.0705230784
a2 = 0.0422820123
a3 = 0.0092705272
a4 = 0.0001520143
a5 = 0.0002765672
a6 = 0.0000430638
In [5]:
myErf :: (Floating d, Ord d, E.Erf d) => d -> d
myErf x = E.erf x
-- myErf x = erfDropIn2 x
In [6]:
-- | Cummulative density function of a normal Gaussion distribution.
normCdf :: (Floating d, Ord d, E.Erf d) => d -> d
normCdf x = 0.5 * (1 + myErf (x / sqrt 2))
-- | Probability density function of a normal Gaussian distribution.
normPdf :: Floating d => d -> d
normPdf x = exp (-x**2 / 2) / sqrt (2 * pi)
In [7]:
-- | Non-normalized normal.
gauss :: Double -> Double
gauss x = exp (-0.5 * x**2)
normPrimitive :: Double -> Double
normPrimitive x = normCdf x * sqrt (2 * pi)
firstPrimitive :: Double -> Double
firstPrimitive x = -gauss x
secondPrimitive :: Double -> Double
secondPrimitive x = normCdf x * sqrt (2 * pi) - x * gauss x
-- | Primitive for integrals with a monomial of degree 2*i + 1.
oddPrimitive :: Integer -> Double -> Double
oddPrimitive i x =
-gauss x
* sum [fromIntegral (doubleFac (2 * i))
/ fromIntegral (doubleFac (2 * j)) * x ** (2 * fromIntegral j)
| j <- [0 .. i]]
-- | Primitive for integrals with monomials of degree 2*i + 2.
evenPrimitive :: Integer -> Double -> Double
evenPrimitive i x =
-gauss x
* sum [fromIntegral (doubleFac (2 * i + 1))
/ fromIntegral (doubleFac (2 * j + 1)) * x ** (2 * fromIntegral j + 1)
| j <- [0 .. i]]
+ fromIntegral (doubleFac (2 * i + 1)) * normCdf x * sqrt(2*pi)
In [8]:
-- Let's test the above primitives
-- Every odd moment should be 0 for a region large enough.
oddPrimitive 3 10 - oddPrimitive 3 (-10)
-- Every even moment p should be (p - 1)!!
-- Remark the sqrt(2*pi) is not included in the primitives.
(evenPrimitive 0 10 - evenPrimitive 0 (-10)) / sqrt(2*pi) -- p == 2
(evenPrimitive 1 10 - evenPrimitive 1 (-10)) / sqrt(2*pi) -- p == 4
doubleFac (2 - 1)
doubleFac (4 - 1)
In [9]:
-- Let's test it agaist some numerically evaluated integrals.
lower = -2.0
upper = 3.0
i = 1
step = 0.001
sum [x**(2*i + 1) * exp (-0.5 * x**2) | x <- [lower, lower+step .. upper]] * step
oddPrimitive i upper - oddPrimitive i lower
In [10]:
-- Let's test it agaist some numerically evaluated integrals.
lower = -2.0
upper = 3.0
i = 0
step = 0.001
sum [x**(2*i + 2) * exp (-0.5 * x**2) | x <- [lower, lower+step .. upper]] * step
evenPrimitive i upper - evenPrimitive i lower
In [11]:
-- | Convenience functions for integrals of monomials and Gaussians.
monomialGauss :: Integer -> Double -> Double -> Double
monomialGauss 0 lower upper = normPrimitive upper - normPrimitive lower
monomialGauss 1 lower upper = firstPrimitive upper - firstPrimitive lower
monomialGauss 2 lower upper = secondPrimitive upper - secondPrimitive lower
monomialGauss n lower upper
| odd n = oddPrimitive n' upper - oddPrimitive n' lower
| otherwise = evenPrimitive (n' - 1) upper - evenPrimitive (n' - 1) lower
where n' = n `div` 2
In [12]:
-- The sqrt(2*pi) is not included in the Gaussian factor.
monomialGauss 4 (-10) 10 / sqrt(2*pi)
doubleFac 3
In [13]:
-- Again do a numerical test to be sure.
testMonomialGauss n = (analytic, numerical)
where
lower = -2.0
upper = 3.0
analytic = monomialGauss n lower upper
numerical = sum [x**fromIntegral n * exp (-0.5 * x**2)
| x <- [lower, lower+step .. upper]] * step
step = 0.001
[testMonomialGauss n | n <- [0 .. 5]]
In [14]:
prior lambda = normPdf ((lambda - 3) / 0.5)
likelihood k lambda = poissonPdf lambda k
In [15]:
-- Set the plotting interval
lower = -2.0
upper = 1.0
plotMsg fs = chart [[(lambda, f lambda) | lambda <- [lower,lower+0.01..upper]]
| f <- fs]
In [16]:
plotMsg [prior, likelihood 2, \ lambda -> prior lambda * likelihood 2 lambda]
In [17]:
normPdfExp pi_ tau x = exp(tau * x - 0.5 * pi_ * x**2)
-- Plot a normal distribution with mean 2 and var 2
plotMsg [ \ x -> normPdf ((x - 2) / sqrt 2) / sqrt 2
, \ x -> normPdfExp 0.5 (0.5 * 2) x / sqrt(2 * pi)
* sqrt 0.5 * exp (-0.5 * 2 * 0.5 * 2) ]
In [18]:
choose :: Integer -> Integer -> Integer
choose n 0 = 1
choose 0 k = 0
choose n k = choose (n-1) (k-1) * n `div` k
integral :: Integer -> Double -> Double -> Double
integral k pi_ tau =
sum [fromIntegral (choose k i)
* (tau - 1) ** fromIntegral (k - i)
* pi_ ** (0.5 * fromIntegral (-2*k + i - 1))
* monomialGauss i (-(tau - 1) / sqrt pi_) 20
| i <- [0 .. k]]
In [19]:
-- Let's check the integral against a numerically calculated one.
rawMsg k pi_ tau x = x**fromIntegral k * exp (-0.5 * pi_ * (x - (tau - 1)/pi_)**2)
integral 1 37.04419504640045 (-31.76430286131713)
-- sum [rawMsg 3 11.74891663176574 (-5.253978919935352) x | x <- [0,0.00005 .. 4]] * 0.00005
step :: Double
step = 0.0001
step * V.sum (V.map (rawMsg 1 37.04419504640045 (-31.76430286131713)) $ V.iterateN (round (4.0 / step)) (+ step) 0)
In [20]:
approx k pi_ tau = (pi_', tau', mu, sigma2)
where
-- Exclude all but the currently selected factor after the projection.
pi_' = newPi_ - pi_
tau' = newTau - tau
partition = integral k pi_ tau
expectedX = integral (k+1) pi_ tau / partition
expectedX2 = integral (k+2) pi_ tau / partition
mu = expectedX
sigma2 = expectedX2 - expectedX**2
-- Map the expected sufficient statistics back to natural parameters.
newPi_ = 1 / sigma2
newTau = mu / sigma2
approx 3 11.74891663176574 (-5.253978919935352)
In [21]:
trueMsg k pi_ tau x
| x > 0 = x**fromIntegral k * exp (-x) * exp (tau * x - 0.5 * pi_ * x**2)
| otherwise = 0
approxMsg k pi_ tau x = exp ((tau + tau') * x - 0.5 * (pi_ + pi_') * x**2) * 0.05
where
(pi_', tau', _, _) = approx k pi_ tau
In [22]:
meanApproxMsg k pi_ tau = (mean, var)
where
lower = -2.0
upper = 3.0
step = 0.005
moments f = sum [f x * approxMsg k pi_ tau x | x <- [lower, lower+step .. upper]] * step
partition = moments $ const 1
mean = moments id / partition
var = moments (**2)/partition - mean**2
meanApproxMsg 3 11.74891663176574 (-5.253978919935352)
In [23]:
meanTrueMsg k pi_ tau = (mean, var)
where
lower = -2.0
upper = 3.0
step = 0.005
moments f = sum [f x * trueMsg k pi_ tau x | x <- [lower, lower+step .. upper]] * step
partition = moments $ const 1
mean = moments id / partition
var = moments (**2)/partition - mean**2
meanTrueMsg 3 11.74891663176574 (-5.253978919935352)
In [24]:
plotMsg' k pi_ tau= plotMsg [ \ x -> trueMsg k pi_ tau x / scaleTrue
, \ x -> approxMsg k pi_ tau x / scaleApprox]
where
lower = -2.0
upper = 3.0
step = 0.005
scale f = sum [f x | x <- [lower, lower+step .. upper]] * step
scaleTrue = scale $ trueMsg k pi_ tau
scaleApprox = scale $ approxMsg k pi_ tau
plotMsg' 1 37.04419504640045 (-31.76430286131713)
In [ ]:
In [ ]: