{-# LANGUAGE
MultiParamTypeClasses, FlexibleInstances, FlexibleContexts,
UndecidableInstances, ForeignFunctionInterface, BangPatterns,
RankNTypes
#-}moduleData.Random.Distribution.Normal(Normal(..),normal,normalT,stdNormal,stdNormalT,doubleStdNormal,floatStdNormal,realFloatStdNormal,normalTail,normalPair,boxMullerNormalPair,knuthPolarNormalPair)whereimportData.Random.Internal.WordsimportData.BitsimportData.Random.SourceimportData.Random.DistributionimportData.Random.Distribution.UniformimportData.Random.Distribution.ZigguratimportData.Random.RVarimportData.Vector.Generic(Vector)importqualifiedData.VectorasVimportqualifiedData.Vector.UnboxedasUVimportData.Number.Erf-- |A random variable that produces a pair of independent-- normally-distributed values.normalPair::(Floatinga,DistributionStdUniforma)=>RVar(a,a)normalPair=boxMullerNormalPair-- |A random variable that produces a pair of independent-- normally-distributed values, computed using the Box-Muller method.-- This algorithm is slightly slower than Knuth's method but using a -- constant amount of entropy (Knuth's method is a rejection method).-- It is also slightly more general (Knuth's method require an 'Ord'-- instance).{-# INLINE boxMullerNormalPair #-}boxMullerNormalPair::(Floatinga,DistributionStdUniforma)=>RVar(a,a)boxMullerNormalPair=dou<-stdUniformt<-stdUniformletr=sqrt(-2*logu)theta=(2*pi)*tx=r*costhetay=r*sinthetareturn(x,y)-- |A random variable that produces a pair of independent-- normally-distributed values, computed using Knuth's polar method.-- Slightly faster than 'boxMullerNormalPair' when it accepts on the -- first try, but does not always do so.{-# INLINE knuthPolarNormalPair #-}knuthPolarNormalPair::(Floatinga,Orda,DistributionUniforma)=>RVar(a,a)knuthPolarNormalPair=dov1<-uniform(-1)1v2<-uniform(-1)1lets=v1*v1+v2*v2ifs>=1thenknuthPolarNormalPairelsereturn$ifs==0then(0,0)elseletscale=sqrt(-2*logs/s)in(v1*scale,v2*scale)-- |Draw from the tail of a normal distribution (the region beyond the provided value){-# INLINE normalTail #-}normalTail::(DistributionStdUniforma,Floatinga,Orda)=>a->RVarTmanormalTailr=gowherego=do!u<-stdUniformTlet!x=logu/r!v<-stdUniformTlet!y=logvifx*x+y+y>0thengoelsereturn(r-x)-- |Construct a 'Ziggurat' for sampling a normal distribution, given-- @logBase 2 c@ and the 'zGetIU' implementation.normalZ::(RealFloata,Erfa,Vectorva,DistributionUniforma,Integralb)=>b->(forallm.RVarTm(Int,a))->ZigguratvanormalZp=mkZigguratRecTruenormalFnormalFInvnormalFIntnormalFVol(2^p)-- | Ziggurat target function (upper half of a non-normalized gaussian PDF)normalF::(Floatinga,Orda)=>a->anormalFx|x<=0=1|otherwise=exp((-0.5)*x*x)-- | inverse of 'normalF'normalFInv::Floatinga=>a->anormalFInvy=sqrt((-2)*logy)-- | integral of 'normalF'normalFInt::(Floatinga,Erfa,Orda)=>a->anormalFIntx|x<=0=0|otherwise=normalFVol*erf(x*sqrt0.5)-- | volume of 'normalF'normalFVol::Floatinga=>anormalFVol=sqrt(0.5*pi)-- |A random variable sampling from the standard normal distribution-- over any 'RealFloat' type (subject to the rest of the constraints --- it builds and uses a 'Ziggurat' internally, which requires the 'Erf'-- class). -- -- Because it computes a 'Ziggurat', it is very expensive to use for-- just one evaluation, or even for multiple evaluations if not used and-- reused monomorphically (to enable the ziggurat table to be let-floated-- out). If you don't know whether your use case fits this description-- then you're probably better off using a different algorithm, such as-- 'boxMullerNormalPair' or 'knuthPolarNormalPair'. And of course if-- you don't need the full generality of this definition then you're much-- better off using 'doubleStdNormal' or 'floatStdNormal'.---- As far as I know, this should be safe to use in any monomorphic-- @Distribution Normal@ instance declaration.realFloatStdNormal::(RealFloata,Erfa,DistributionUniforma)=>RVarTmarealFloatStdNormal=runZiggurat(normalZpgetIU`asTypeOf`(undefined::ZigguratV.Vectora))wherep::Intp=6getIU::(Numa,DistributionUniforma)=>RVarTm(Int,a)getIU=doi<-getRandomWord8u<-uniformT(-1)1return(fromIntegrali.&.(2^p-1),u)-- |A random variable sampling from the standard normal distribution-- over the 'Double' type.doubleStdNormal::RVarTmDoubledoubleStdNormal=runZigguratdoubleStdNormalZ-- doubleStdNormalC must not be over 2^12 if using wordToDoubleWithExcessdoubleStdNormalC::IntdoubleStdNormalC=512doubleStdNormalR,doubleStdNormalV::DoubledoubleStdNormalR=3.852046150368388doubleStdNormalV=2.4567663515413507e-3{-# NOINLINE doubleStdNormalZ #-}doubleStdNormalZ::ZigguratUV.VectorDoubledoubleStdNormalZ=mkZiggurat_TruenormalFnormalFInvdoubleStdNormalCdoubleStdNormalRdoubleStdNormalVgetIU(normalTaildoubleStdNormalR)wheregetIU::RVarTm(Int,Double)getIU=do!w<-getRandomWord64let(u,i)=wordToDoubleWithExcesswreturn$!(fromIntegrali.&.(doubleStdNormalC-1),u+u-1)-- |A random variable sampling from the standard normal distribution-- over the 'Float' type.floatStdNormal::RVarTmFloatfloatStdNormal=runZigguratfloatStdNormalZ-- floatStdNormalC must not be over 2^9 if using word32ToFloatWithExcessfloatStdNormalC::IntfloatStdNormalC=512floatStdNormalR,floatStdNormalV::FloatfloatStdNormalR=3.852046150368388floatStdNormalV=2.4567663515413507e-3{-# NOINLINE floatStdNormalZ #-}floatStdNormalZ::ZigguratUV.VectorFloatfloatStdNormalZ=mkZiggurat_TruenormalFnormalFInvfloatStdNormalCfloatStdNormalRfloatStdNormalVgetIU(normalTailfloatStdNormalR)wheregetIU::RVarTm(Int,Float)getIU=do!w<-getRandomWord32let(u,i)=word32ToFloatWithExcesswreturn(fromIntegrali.&.(floatStdNormalC-1),u+u-1)normalCdf::(Reala)=>a->a->a->DoublenormalCdfmsx=normcdf((realToFracx-realToFracm)/realToFracs)-- |A specification of a normal distribution over the type 'a'.dataNormala-- |The \"standard\" normal distribution - mean 0, stddev 1=StdNormal-- |@Normal m s@ is a normal distribution with mean @m@ and stddev @sd@.|Normalaa-- mean, sdinstanceDistributionNormalDoublewherervarTStdNormal=doubleStdNormalrvarT(Normalms)=dox<-doubleStdNormalreturn(x*s+m)instanceDistributionNormalFloatwherervarTStdNormal=floatStdNormalrvarT(Normalms)=dox<-floatStdNormalreturn(x*s+m)instance(Reala,DistributionNormala)=>CDFNormalawherecdfStdNormal=normalCdf01cdf(Normalms)=normalCdfms{-# SPECIALIZE stdNormal :: RVar Double #-}{-# SPECIALIZE stdNormal :: RVar Float #-}-- |'stdNormal' is a normal variable with distribution 'StdNormal'.stdNormal::DistributionNormala=>RVarastdNormal=rvarStdNormal-- |'stdNormalT' is a normal process with distribution 'StdNormal'.stdNormalT::DistributionNormala=>RVarTmastdNormalT=rvarTStdNormal-- |@normal m s@ is a random variable with distribution @'Normal' m s@.normal::DistributionNormala=>a->a->RVaranormalms=rvar(Normalms)-- |@normalT m s@ is a random process with distribution @'Normal' m s@.normalT::DistributionNormala=>a->a->RVarTmanormalTms=rvarT(Normalms)