mwc-probability

commit 5b96ae405e27bdee1cd2ab10209e3861fcebc212
parent 782945d540278ed52d5e3cf49a64a8b5bd448e4f
Author: Marco Zocca <marco.zocca@recordunion.com>
Date:   Wed,  9 May 2018 13:01:29 +0200

add inverse Gaussian

Diffstat:
M
CHANGELOG
 | 
3
+++
M
mwc-probability.cabal
 | 
2
+-
M
src/System/Random/MWC/Probability.hs
 | 
17
+++++++++++++++++

3 files changed, 21 insertions(+), 1 deletion(-)
diff --git a/CHANGELOG b/CHANGELOG
@@ -1,5 +1,8 @@
 	# Changelog
 
+	- 2.0.3 (2018-05-09)
+	* Add inverse Gaussian (Wald) distribution
+	
 	- 2.0.2 (2018-01-30)
 	* Add negative binomial distribution
 
diff --git a/mwc-probability.cabal b/mwc-probability.cabal
@@ -1,5 +1,5 @@
 name:                mwc-probability
-version:             2.0.2
+version:             2.0.3
 homepage:            http://github.com/jtobin/mwc-probability
 license:             MIT
 license-file:        LICENSE
diff --git a/src/System/Random/MWC/Probability.hs b/src/System/Random/MWC/Probability.hs
@@ -73,6 +73,7 @@ module System.Random.MWC.Probability (
   , isoNormal    
   , logNormal
   , exponential
+  , inverseGaussian
   , laplace
   , gamma
   , inverseGamma
@@ -334,6 +335,22 @@ isoNormal
 isoNormal ms sd = traverse (`normal` sd) ms
 {-# INLINABLE isoNormal #-}
 
+-- | The inverse Gaussian (also known as Wald) distribution.
+--
+-- Both the mean parameter 'my' and the shape parameter 'lambda' must be positive.
+inverseGaussian :: PrimMonad m => Double -> Double -> Prob m Double
+inverseGaussian lambda my = do
+  nu <- standardNormal
+  let y = nu ** 2
+      x = my + 1 / (2 * lambda) * (my ** 2 * y + my * sqrt (4 * my * lambda * y) + my ** 2  * y ** 2)
+      thresh = my / (my + x)
+  z <- uniform
+  if z <= thresh
+    then return x
+    else return (my ** 2 / x)
+{-# INLINABLE inverseGaussian #-}    
+  
+
 -- | The Poisson distribution.
 poisson :: PrimMonad m => Double -> Prob m Int
 poisson l = Prob $ genFromTable table where