basic polar method rejection sampler for normal distribution
authorCarter Tazio Schonwald <carter.schonwald@gmail.com>
Sat, 27 May 2017 21:07:30 +0000 (17:07 -0400)
committerCarter Tazio Schonwald <carter.schonwald@gmail.com>
Sat, 27 May 2017 21:07:30 +0000 (17:07 -0400)
random.cabal
src/Data/Distribution/FloatingInterval.hs
src/Data/Distribution/Normal.hs

index 31944b7..597865c 100644 (file)
@@ -70,6 +70,7 @@ library
   build-depends:       base >=4.8 && <4.12
                       ,ghc-prim
                       ,entropy == 0.3.*
+                      ,numeric-extras == 0.1.*
                   -- entropy will later be folded into random, probably
 
   -- Directories containing source files.
index 194e60b..b9e8860 100644 (file)
@@ -74,6 +74,7 @@ is then xor'd against the corrected unit interval number in this specific
 
 
 extracted docs from the original site:
+"""
    #include <stdint.h>
 
     (x >> 11) * (1. / (UINT64_C(1) << 53))
@@ -88,6 +89,9 @@ An alternative, faster multiplication-free operation is
        return u.d - 1.0;
     }
 The code above cooks up by bit manipulation a real number in the interval [1..2), and then subtracts one to obtain a real number in the interval [0..1). If x is chosen uniformly among 64-bit integers, d is chosen uniformly among dyadic rationals of the form k / 2−52.
+"""
+
+
  -}
 sampleUnitIntervalDoubleReallySloppyM ::  forall m . Monad m => (m Word64) -> m Double
 sampleUnitIntervalDoubleReallySloppyM mword = do
index de28105..afdb433 100644 (file)
@@ -1,5 +1,36 @@
+{-# LANGUAGE  ScopedTypeVariables #-}
 module Data.Distribution.Normal where
 
-{-
 
+import Numeric.Extras (hypot)
+{- | For now this will be using the Marsaglia polar method,
+though might be better methods to use in exchange for a teeny bit more complexity.
+
+NB: tail distribution depends on quality of the unit interval generator
+
+This implementation only returns one of the two random variates, and thus
+only needs to generate real samples from [+0,1) rather than (-1,1)
+
+
+if using the x / 2^-53 style uniform interval, you cant get a result r
+such that abs(r) > 15, though thats pretty extreme in the tail distribution
 -}
+unitNormalPolarMethodM :: forall m . Monad m => m Double -> m Bool -> m Double
+unitNormalPolarMethodM unitSampler boolSampler = getSample
+  where
+    getSample :: m Double
+    getSample = do
+      x <- unitSampler
+      y <- unitSampler
+      sqrtSumSq <- return $ hypot x y
+      straightSum <- return $ x*x + y * y
+      if straightSum > 1 || straightSum == 0
+        --- the usual condition is  x^2 + y^2 > 1, but the same bound holds for the sqrt thereof
+        then getSample
+        else
+          do
+            boolS <- boolSampler
+            signed <- return $ if boolS then 1 else -1
+            return $! signed * x * (sqrt ( -2 * log straightSum) /  sqrtSumSq)
+
+