author Carter Tazio Schonwald Sat, 27 May 2017 21:07:30 +0000 (17:07 -0400) committer Carter Tazio Schonwald Sat, 27 May 2017 21:07:30 +0000 (17:07 -0400)

index 31944b7..597865c 100644 (file)
@@ -70,6 +70,7 @@ library
build-depends:       base >=4.8 && <4.12
,ghc-prim
,entropy == 0.3.*
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.
-- 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:

extracted docs from the original site:
+"""
#include <stdint.h>

(x >> 11) * (1. / (UINT64_C(1) << 53))
#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.
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
-}
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

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)
+
+