diff --git a/NEWS b/NEWS --- a/NEWS +++ b/NEWS @@ -1,3 +1,11 @@ +Changes in 0.4.2 +* Depend on the vector package. +* Felipe Lessa added beta, logistic, Pareto, Weibull, gamma + and Dirichlet distributions. +* Remove -fPIC flag when compiling C code; + fixes https://github.com/patperry/hs-gsl-random/issues/1 + (thanks to Dmitry Tsygankov) + Changes in 0.4.1 * Require cabal >= 1.8 * Fix framework detection bug in configure script on Mac OS @@ -12,7 +20,6 @@ Changes in 0.3.1 * Add exponential distribution -Changes in 0.3: - +Changes in 0.3 * Add Levy alpha-Stable distributions * Fix GHCi bus error on x86-64 diff --git a/gsl-random.cabal b/gsl-random.cabal --- a/gsl-random.cabal +++ b/gsl-random.cabal @@ -1,6 +1,6 @@ name: gsl-random -version: 0.4.1 -homepage: http://github.com/patperry/gsl-random +version: 0.4.2 +homepage: http://github.com/patperry/hs-gsl-random synopsis: Bindings the the GSL random number generation facilities. description: Bindings to the GNU Scientific Library random and quasi-random number @@ -19,7 +19,7 @@ tested-with: GHC == 6.12.3 extra-source-files: NEWS - + library hs-source-dirs: lib exposed-modules: GSL.Random.Gen @@ -30,7 +30,6 @@ c-sources: cbits/qrng_extras.c ghc-options: -Wall - if !os(darwin) - cc-options: -fPIC extensions: ForeignFunctionInterface - build-depends: base >= 4 && < 5 + build-depends: base >= 4 && < 5 + , vector >= 0.7.0.1 && < 0.8 diff --git a/lib/GSL/Random/Dist.hs b/lib/GSL/Random/Dist.hs --- a/lib/GSL/Random/Dist.hs +++ b/lib/GSL/Random/Dist.hs @@ -35,7 +35,7 @@ getUGaussian, getUGaussianRatioMethod, - + -- * The Flat (Uniform) Distribution flatPdf, @@ -43,7 +43,7 @@ flatQ, flatPInv, flatQInv, - + getFlat, -- * The Exponential Distribution @@ -53,13 +53,13 @@ exponentialQ, exponentialPInv, exponentialQInv, - + getExponential, -- * The Levy alpha-Stable Distributions getLevy, getLevySkew, - + -- * The Poisson Distribution poissonPdf, @@ -67,7 +67,7 @@ poissonQ, getPoisson, - + -- * The Cauchy Distribution getCauchy, @@ -75,18 +75,73 @@ cauchyP, cauchyQ, cauchyPInv, - cauchyQInv + cauchyQInv, + -- * The Beta Distribution + getBeta, + + betaPdf, + betaP, + betaQ, + betaPInv, + betaQInv, + + -- * The Logistic Distribution + getLogistic, + + logisticPdf, + logisticP, + logisticQ, + logisticPInv, + logisticQInv, + + -- * The Pareto Distribution + getPareto, + + paretoPdf, + paretoP, + paretoQ, + paretoPInv, + paretoQInv, + + -- * The Weibull Distribution + getWeibull, + + weibullPdf, + weibullP, + weibullQ, + weibullPInv, + weibullQInv, + + -- * The Gamma Distribution + getGamma, + getGammaKnuth, + + gammaPdf, + gammaP, + gammaQ, + gammaPInv, + gammaQInv, + + -- * The Dirichlet Distribution + getDirichlet, + + dirichletPdf, + dirichletLnPdf, + ) where -import Control.Monad -import Foreign.C.Types ( CUInt, CDouble ) -import Foreign.ForeignPtr ( withForeignPtr ) +import Control.Applicative ( (<$>) ) +import Foreign.C.Types ( CUInt, CDouble, CSize ) +import Foreign.ForeignPtr ( withForeignPtr, mallocForeignPtrArray ) import Foreign.Ptr ( Ptr ) +import System.IO.Unsafe ( unsafePerformIO ) +import qualified Data.Vector.Storable as VS + import GSL.Random.Gen.Internal ( RNG(..) ) --- | @gaussianPdf x sigma@ computes the probabililty density p(x) for +-- | @gaussianPdf x sigma@ computes the probabililty density p(x) for -- a Gaussian distribution with mean @0@ and standard deviation @sigma@. gaussianPdf :: Double -> Double -> Double gaussianPdf = liftDouble2 gsl_ran_gaussian_pdf @@ -94,7 +149,7 @@ foreign import ccall unsafe "gsl/gsl_randist.h" gsl_ran_gaussian_pdf :: CDouble -> CDouble -> CDouble --- | @gaussianP x sigma@ computes the cumulative distribution function P(x) for +-- | @gaussianP x sigma@ computes the cumulative distribution function P(x) for -- a Gaussian distribution with mean @0@ and standard deviation @sigma@. gaussianP :: Double -> Double -> Double gaussianP = liftDouble2 gsl_cdf_gaussian_P @@ -102,7 +157,7 @@ foreign import ccall unsafe "gsl/gsl_randist.h" gsl_cdf_gaussian_P :: CDouble -> CDouble -> CDouble --- | @gaussianQ x sigma@ computes the cumulative distribution function Q(x) for +-- | @gaussianQ x sigma@ computes the cumulative distribution function Q(x) for -- a Gaussian distribution with mean @0@ and standard deviation @sigma@. gaussianQ :: Double -> Double -> Double gaussianQ = liftDouble2 gsl_cdf_gaussian_Q @@ -110,7 +165,7 @@ foreign import ccall unsafe "gsl/gsl_randist.h" gsl_cdf_gaussian_Q :: CDouble -> CDouble -> CDouble --- | @gaussianPInv p sigma@ computes the inverse of the cumulative distribution +-- | @gaussianPInv p sigma@ computes the inverse of the cumulative distribution -- function of a Gaussian distribution with mean @0@ and standard deviation -- @sigma@. It returns @x@ such that @P(x) = p@. gaussianPInv :: Double -> Double -> Double @@ -119,7 +174,7 @@ foreign import ccall unsafe "gsl/gsl_randist.h" gsl_cdf_gaussian_Pinv :: CDouble -> CDouble -> CDouble --- | @gaussianPInv q sigma@ computes the inverse of the cumulative distribution +-- | @gaussianPInv q sigma@ computes the inverse of the cumulative distribution -- function of a Gaussian distribution with mean @0@ and standard deviation -- @sigma@. It returns @x@ such that @Q(x) = q@. gaussianQInv :: Double -> Double -> Double @@ -129,39 +184,36 @@ gsl_cdf_gaussian_Qinv :: CDouble -> CDouble -> CDouble -- | @getGaussian r sigma@ gets a normal random variable with mean --- @0@ and standard deviation @sigma@. +-- @0@ and standard deviation @sigma@. -- This uses the Box-Mueller algorithm. getGaussian :: RNG -> Double -> IO Double -getGaussian = getGaussianHelp gsl_ran_gaussian +getGaussian = liftRan1 gsl_ran_gaussian +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_ran_gaussian :: Ptr () -> CDouble -> IO CDouble + -- | @getGaussianZiggurat r sigma@ gets a normal random variable with mean --- @0@ and standard deviation @sigma@. +-- @0@ and standard deviation @sigma@. -- This uses the Marsaglia-Tsang ziggurat algorithm. getGaussianZiggurat :: RNG -> Double -> IO Double -getGaussianZiggurat = getGaussianHelp gsl_ran_gaussian_ziggurat +getGaussianZiggurat = liftRan1 gsl_ran_gaussian_ziggurat +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_ran_gaussian_ziggurat :: Ptr () -> CDouble -> IO CDouble + -- | @getGaussianRatioMethod r sigma@ gets a normal random variable with mean --- @0@ and standard deviation @sigma@. +-- @0@ and standard deviation @sigma@. -- This uses the Kinderman-Monahan-Leva ratio method. getGaussianRatioMethod:: RNG -> Double -> IO Double -getGaussianRatioMethod = getGaussianHelp gsl_ran_gaussian_ratio_method - -getGaussianHelp :: (Ptr () -> CDouble -> IO CDouble) - -> RNG -> Double -> IO Double -getGaussianHelp ran_gaussian (MkRNG fptr) sigma = - let sigma' = realToFrac sigma - in withForeignPtr fptr $ \ptr -> do - x <- ran_gaussian ptr sigma' - return $ realToFrac x +getGaussianRatioMethod = liftRan1 gsl_ran_gaussian_ratio_method foreign import ccall unsafe "gsl/gsl_randist.h" - gsl_ran_gaussian :: Ptr () -> CDouble -> IO CDouble -foreign import ccall unsafe "gsl/gsl_randist.h" - gsl_ran_gaussian_ziggurat :: Ptr () -> CDouble -> IO CDouble -foreign import ccall unsafe "gsl/gsl_randist.h" gsl_ran_gaussian_ratio_method :: Ptr () -> CDouble -> IO CDouble --- | @ugaussianPdf x@ computes the probabililty density p(x) for + + + +-- | @ugaussianPdf x@ computes the probabililty density p(x) for -- a Gaussian distribution with mean @0@ and standard deviation @1@. ugaussianPdf :: Double -> Double ugaussianPdf = liftDouble gsl_ran_ugaussian_pdf @@ -169,7 +221,7 @@ foreign import ccall unsafe "gsl/gsl_randist.h" gsl_ran_ugaussian_pdf :: CDouble -> CDouble --- | @ugaussianP x@ computes the cumulative distribution function P(x) for +-- | @ugaussianP x@ computes the cumulative distribution function P(x) for -- a Gaussian distribution with mean @0@ and standard deviation @1@. ugaussianP :: Double -> Double ugaussianP = liftDouble gsl_cdf_ugaussian_P @@ -177,7 +229,7 @@ foreign import ccall unsafe "gsl/gsl_randist.h" gsl_cdf_ugaussian_P :: CDouble -> CDouble --- | @ugaussianQ x@ computes the cumulative distribution function Q(x) for +-- | @ugaussianQ x@ computes the cumulative distribution function Q(x) for -- a Gaussian distribution with mean @0@ and standard deviation @1@. ugaussianQ :: Double -> Double ugaussianQ = liftDouble gsl_cdf_ugaussian_Q @@ -185,7 +237,7 @@ foreign import ccall unsafe "gsl/gsl_randist.h" gsl_cdf_ugaussian_Q :: CDouble -> CDouble --- | @ugaussianPInv p@ computes the inverse of the cumulative distribution +-- | @ugaussianPInv p@ computes the inverse of the cumulative distribution -- function of a Gaussian distribution with mean @0@ and standard deviation -- @1@. It returns @x@ such that @P(x) = p@. ugaussianPInv :: Double -> Double @@ -194,7 +246,7 @@ foreign import ccall unsafe "gsl/gsl_randist.h" gsl_cdf_ugaussian_Pinv :: CDouble -> CDouble --- | @ugaussianPInv q@ computes the inverse of the cumulative distribution +-- | @ugaussianPInv q@ computes the inverse of the cumulative distribution -- function of a Gaussian distribution with mean @0@ and standard deviation -- @1@. It returns @x@ such that @Q(x) = q@. ugaussianQInv :: Double -> Double @@ -204,69 +256,27 @@ gsl_cdf_ugaussian_Qinv :: CDouble -> CDouble -- | @getUGaussian r@ gets a normal random variable with mean --- @0@ and standard deviation @1@. +-- @0@ and standard deviation @1@. -- This uses the Box-Mueller algorithm. getUGaussian :: RNG -> IO Double -getUGaussian = getUGaussianHelp gsl_ran_ugaussian +getUGaussian = liftRan0 gsl_ran_ugaussian +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_ran_ugaussian :: Ptr () -> IO CDouble + -- | @getUGaussianRatioMethod r@ gets a normal random variable with mean --- @0@ and standard deviation @1@. +-- @0@ and standard deviation @1@. -- This uses the Kinderman-Monahan-Leva ratio method. getUGaussianRatioMethod:: RNG -> IO Double -getUGaussianRatioMethod = getUGaussianHelp gsl_ran_ugaussian_ratio_method - -getUGaussianHelp :: (Ptr () -> IO CDouble) - -> RNG -> IO Double -getUGaussianHelp ran_ugaussian (MkRNG fptr) = - withForeignPtr fptr $ \ptr -> do - x <- ran_ugaussian ptr - return $ realToFrac x +getUGaussianRatioMethod = liftRan0 gsl_ran_ugaussian_ratio_method foreign import ccall unsafe "gsl/gsl_randist.h" - gsl_ran_ugaussian :: Ptr () -> IO CDouble -foreign import ccall unsafe "gsl/gsl_randist.h" gsl_ran_ugaussian_ratio_method :: Ptr () -> IO CDouble --- | @getExponential r mu@ gets a random exponential with mean @mu@. -getExponential :: RNG -> Double -> IO Double -getExponential (MkRNG f) mu = withForeignPtr f $ \p -> - liftM realToFrac $ gsl_ran_exponential p (realToFrac mu) - -foreign import ccall unsafe "gsl/gsl_randist.h" - gsl_ran_exponential :: Ptr () -> CDouble -> IO CDouble --- | @exponentialPdf x mu@ computes the density at @x@ of an exponential --- with mean @mu@. -exponentialPdf :: Double -> Double -> Double -exponentialPdf = liftDouble2 gsl_ran_exponential_pdf -foreign import ccall unsafe "gsl/gsl_randist.h" - gsl_ran_exponential_pdf :: CDouble -> CDouble -> CDouble -exponentialP :: Double -> Double -> Double -exponentialP = liftDouble2 gsl_cdf_exponential_P -foreign import ccall unsafe "gsl/gsl_randist.h" - gsl_cdf_exponential_P :: CDouble -> CDouble -> CDouble - -exponentialQ :: Double -> Double -> Double -exponentialQ = liftDouble2 gsl_cdf_exponential_Q - -foreign import ccall unsafe "gsl/gsl_randist.h" - gsl_cdf_exponential_Q :: CDouble -> CDouble -> CDouble - -exponentialPInv :: Double -> Double -> Double -exponentialPInv = liftDouble2 gsl_cdf_exponential_Pinv - -foreign import ccall unsafe "gsl/gsl_randist.h" - gsl_cdf_exponential_Pinv :: CDouble -> CDouble -> CDouble - -exponentialQInv :: Double -> Double -> Double -exponentialQInv = liftDouble2 gsl_cdf_exponential_Qinv - -foreign import ccall unsafe "gsl/gsl_randist.h" - gsl_cdf_exponential_Qinv :: CDouble -> CDouble -> CDouble - -- | @flatPdf x a b@ computes the probability density @p(x)@ at @x@ for -- a uniform distribution from @a@ to @b@. flatPdf :: Double -> Double -> Double -> Double @@ -307,24 +317,62 @@ -- | @getFlat r a b@ gets a value uniformly chosen in @[a,b)@. getFlat :: RNG -> Double -> Double -> IO (Double) -getFlat (MkRNG fptr) a b = - let a' = realToFrac a - b' = realToFrac b - in withForeignPtr fptr $ \ptr -> do - x <- gsl_ran_flat ptr a' b' - return $ realToFrac x - +getFlat = liftRan2 gsl_ran_flat + foreign import ccall unsafe "gsl/gsl_randist.h" gsl_ran_flat :: Ptr () -> CDouble -> CDouble -> IO CDouble + + + +-- | @getExponential r mu@ gets a random exponential with mean @mu@. +getExponential :: RNG -> Double -> IO Double +getExponential = liftRan1 gsl_ran_exponential + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_ran_exponential :: Ptr () -> CDouble -> IO CDouble + +-- | @exponentialPdf x mu@ computes the density at @x@ of an exponential +-- with mean @mu@. +exponentialPdf :: Double -> Double -> Double +exponentialPdf = liftDouble2 gsl_ran_exponential_pdf + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_ran_exponential_pdf :: CDouble -> CDouble -> CDouble + +exponentialP :: Double -> Double -> Double +exponentialP = liftDouble2 gsl_cdf_exponential_P + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_exponential_P :: CDouble -> CDouble -> CDouble + +exponentialQ :: Double -> Double -> Double +exponentialQ = liftDouble2 gsl_cdf_exponential_Q + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_exponential_Q :: CDouble -> CDouble -> CDouble + +exponentialPInv :: Double -> Double -> Double +exponentialPInv = liftDouble2 gsl_cdf_exponential_Pinv + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_exponential_Pinv :: CDouble -> CDouble -> CDouble + +exponentialQInv :: Double -> Double -> Double +exponentialQInv = liftDouble2 gsl_cdf_exponential_Qinv + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_exponential_Qinv :: CDouble -> CDouble -> CDouble + + + + -- | @getLevy r c alpha@ gets a variate from the Levy symmetric stable -- distribution with scale @c@ and exponent @alpha@. The algorithm only -- works for @0 <= alpha <= 2@. getLevy :: RNG -> Double -> Double -> IO (Double) -getLevy (MkRNG f) c alpha = - withForeignPtr f $ \p -> - realToFrac `fmap` gsl_ran_levy p (realToFrac c) (realToFrac alpha) - +getLevy = liftRan2 gsl_ran_levy + foreign import ccall unsafe "gsl/gsl_randist.h" gsl_ran_levy :: Ptr () -> CDouble -> CDouble -> IO CDouble @@ -333,14 +381,16 @@ -- @beta@. The skewness parameter must lie in the range @[-1,1]@. The -- algorithm only works for @0 <= alpha <= 2@. getLevySkew :: RNG -> Double -> Double -> Double -> IO (Double) -getLevySkew (MkRNG f) c alpha beta = - withForeignPtr f $ \p -> - realToFrac `fmap` gsl_ran_levy_skew p (realToFrac c) (realToFrac alpha) (realToFrac beta) - +getLevySkew = liftRan3 gsl_ran_levy_skew + foreign import ccall unsafe "gsl/gsl_randist.h" gsl_ran_levy_skew :: Ptr () -> CDouble -> CDouble -> CDouble -> IO CDouble --- | @poissonPdf k mu@ evaluates the probability density @p(k)@ at @k@ for + + + + +-- | @poissonPdf k mu@ evaluates the probability density @p(k)@ at @k@ for -- a Poisson distribution with mean @mu@. poissonPdf :: Int -> Double -> Double poissonPdf k = liftDouble $ gsl_ran_poisson_pdf (fromIntegral k) @@ -348,7 +398,7 @@ foreign import ccall unsafe "gsl/gsl_randist.h" gsl_ran_poisson_pdf :: CUInt -> CDouble -> CDouble --- | @poissonP k mu@ evaluates the cumulative distribution function @P(k)@ +-- | @poissonP k mu@ evaluates the cumulative distribution function @P(k)@ -- at @k@ for a Poisson distribution with mean @mu@. poissonP :: Int -> Double -> Double poissonP k = liftDouble $ gsl_cdf_poisson_P (fromIntegral k) @@ -356,7 +406,7 @@ foreign import ccall unsafe "gsl/gsl_randist.h" gsl_cdf_poisson_P :: CUInt -> CDouble -> CDouble --- | @poissonQ k mu@ evaluates the cumulative distribution function @Q(k)@ +-- | @poissonQ k mu@ evaluates the cumulative distribution function @Q(k)@ -- at @k@ for a Poisson distribution with mean @mu@. poissonQ :: Int -> Double -> Double poissonQ k = liftDouble $ gsl_cdf_poisson_Q (fromIntegral k) @@ -364,20 +414,21 @@ foreign import ccall unsafe "gsl/gsl_randist.h" gsl_cdf_poisson_Q :: CUInt -> CDouble -> CDouble - + -- | @getPoisson r mu@ gets a poisson random variable with mean @mu@. getPoisson :: RNG -> Double -> IO Int getPoisson (MkRNG fptr) mu = let mu' = realToFrac mu in withForeignPtr fptr $ \ptr -> do - x <- gsl_ran_poisson ptr mu' + x <- gsl_ran_poisson ptr mu' return $ fromIntegral x - + foreign import ccall unsafe "gsl/gsl_randist.h" gsl_ran_poisson :: Ptr () -> CDouble -> IO CUInt + -- | @cauchyPdf x a@ evaluates the probability density @p(x)@ at @x@ -- for a Cauchy distribution with scale parameter @a@. The density -- is given by @p(x) dx = { 1 \over a\pi (1 + (x/a^2)) } dx@. @@ -389,8 +440,7 @@ -- | @getCauchy r a@ gets a random cauchy with scale @a@. getCauchy :: RNG -> Double -> IO Double -getCauchy (MkRNG f) a = withForeignPtr f $ \p -> - liftM realToFrac $ gsl_ran_cauchy p (realToFrac a) +getCauchy = liftRan1 gsl_ran_cauchy foreign import ccall unsafe "gsl/gsl_randist.h" gsl_ran_cauchy :: Ptr () -> CDouble -> IO CDouble @@ -421,18 +471,326 @@ -liftDouble :: (CDouble -> CDouble) + +-- | @betaPdf x a b@ evaluates the probability density @p(x)@ at @x@ +-- for a Beta distribution with parameters @a@ and @b@. The density +-- is given by @p(x) dx = {\Gamma(a+b) \over \Gamma(a) \Gamma(b)} x^{a-1} (1-x)^{b-1} dx@ +-- for @0 <= x <= 1@. +betaPdf :: Double -> Double -> Double -> Double +betaPdf = liftDouble3 gsl_ran_beta_pdf + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_ran_beta_pdf :: CDouble -> CDouble -> CDouble -> CDouble + +-- | @getBeta r a b@ gets a random beta with parameters @a@ and @b@. +getBeta :: RNG -> Double -> Double -> IO Double +getBeta = liftRan2 gsl_ran_beta + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_ran_beta :: Ptr () -> CDouble -> CDouble -> IO CDouble + +betaP :: Double -> Double -> Double -> Double +betaP = liftDouble3 gsl_cdf_beta_P + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_beta_P :: CDouble -> CDouble -> CDouble -> CDouble + +betaQ :: Double -> Double -> Double -> Double +betaQ = liftDouble3 gsl_cdf_beta_Q + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_beta_Q :: CDouble -> CDouble -> CDouble -> CDouble + +betaPInv :: Double -> Double -> Double -> Double +betaPInv = liftDouble3 gsl_cdf_beta_Pinv + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_beta_Pinv :: CDouble -> CDouble -> CDouble -> CDouble + +betaQInv :: Double -> Double -> Double -> Double +betaQInv = liftDouble3 gsl_cdf_beta_Qinv + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_beta_Qinv :: CDouble -> CDouble -> CDouble -> CDouble + + + + +-- | @logisticPdf x a@ evaluates the probability density @p(x)@ at @x@ +-- for a logistic distribution with scale parameter @a@. The density +-- is given by @p(x) dx = { \exp(-x/a) \over a (1 + \exp(-x/a))^2 } dx@. +logisticPdf :: Double -> Double -> Double +logisticPdf = liftDouble2 gsl_ran_logistic_pdf + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_ran_logistic_pdf :: CDouble -> CDouble -> CDouble + +-- | @getLogistic r a@ gets a random logistic with scale @a@. +getLogistic :: RNG -> Double -> IO Double +getLogistic = liftRan1 gsl_ran_logistic + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_ran_logistic :: Ptr () -> CDouble -> IO CDouble + +logisticP :: Double -> Double -> Double +logisticP = liftDouble2 gsl_cdf_logistic_P + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_logistic_P :: CDouble -> CDouble -> CDouble + +logisticQ :: Double -> Double -> Double +logisticQ = liftDouble2 gsl_cdf_logistic_Q + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_logistic_Q :: CDouble -> CDouble -> CDouble + +logisticPInv :: Double -> Double -> Double +logisticPInv = liftDouble2 gsl_cdf_logistic_Pinv + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_logistic_Pinv :: CDouble -> CDouble -> CDouble + +logisticQInv :: Double -> Double -> Double +logisticQInv = liftDouble2 gsl_cdf_logistic_Qinv + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_logistic_Qinv :: CDouble -> CDouble -> CDouble + + + + +-- | @paretoPdf x a b@ evaluates the probability density @p(x)@ at @x@ +-- for a Pareto distribution with exponent @a@ and scale @b@. The density +-- is given by @p(x) dx = (a/b) / (x/b)^{a+1} dx@ for @x >= b@. +paretoPdf :: Double -> Double -> Double -> Double +paretoPdf = liftDouble3 gsl_ran_pareto_pdf + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_ran_pareto_pdf :: CDouble -> CDouble -> CDouble -> CDouble + +-- | @getPareto r a b@ gets a random Pareto with exponent @a@ and scale @b@. +getPareto :: RNG -> Double -> Double -> IO Double +getPareto = liftRan2 gsl_ran_pareto + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_ran_pareto :: Ptr () -> CDouble -> CDouble -> IO CDouble + +paretoP :: Double -> Double -> Double -> Double +paretoP = liftDouble3 gsl_cdf_pareto_P + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_pareto_P :: CDouble -> CDouble -> CDouble -> CDouble + +paretoQ :: Double -> Double -> Double -> Double +paretoQ = liftDouble3 gsl_cdf_pareto_Q + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_pareto_Q :: CDouble -> CDouble -> CDouble -> CDouble + +paretoPInv :: Double -> Double -> Double -> Double +paretoPInv = liftDouble3 gsl_cdf_pareto_Pinv + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_pareto_Pinv :: CDouble -> CDouble -> CDouble -> CDouble + +paretoQInv :: Double -> Double -> Double -> Double +paretoQInv = liftDouble3 gsl_cdf_pareto_Qinv + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_pareto_Qinv :: CDouble -> CDouble -> CDouble -> CDouble + + + + +-- | @weibullPdf x a b@ evaluates the probability density @p(x)@ at @x@ +-- for a Weibull distribution with scale @a@ and exponent @b@. The density +-- is given by @p(x) dx = {b \over a^b} x^{b-1} \exp(-(x/a)^b) dx@ for @x >= 0@. +weibullPdf :: Double -> Double -> Double -> Double +weibullPdf = liftDouble3 gsl_ran_weibull_pdf + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_ran_weibull_pdf :: CDouble -> CDouble -> CDouble -> CDouble + +-- | @getWeibull r a b@ gets a random Weibull with scale @a@ and exponent @b@. +getWeibull :: RNG -> Double -> Double -> IO Double +getWeibull = liftRan2 gsl_ran_weibull + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_ran_weibull :: Ptr () -> CDouble -> CDouble -> IO CDouble + +weibullP :: Double -> Double -> Double -> Double +weibullP = liftDouble3 gsl_cdf_weibull_P + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_weibull_P :: CDouble -> CDouble -> CDouble -> CDouble + +weibullQ :: Double -> Double -> Double -> Double +weibullQ = liftDouble3 gsl_cdf_weibull_Q + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_weibull_Q :: CDouble -> CDouble -> CDouble -> CDouble + +weibullPInv :: Double -> Double -> Double -> Double +weibullPInv = liftDouble3 gsl_cdf_weibull_Pinv + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_weibull_Pinv :: CDouble -> CDouble -> CDouble -> CDouble + +weibullQInv :: Double -> Double -> Double -> Double +weibullQInv = liftDouble3 gsl_cdf_weibull_Qinv + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_weibull_Qinv :: CDouble -> CDouble -> CDouble -> CDouble + + + + +-- | @gammaPdf x a b@ evaluates the probability density @p(x)@ at @x@ +-- for a gamma distribution with parameters @a@ and @b@. The density +-- is given by @p(x) dx = p(x) dx = {1 \over \Gamma(a) b^a} x^{a-1} e^{-x/b} dx@ +-- for @x > 0@. +gammaPdf :: Double -> Double -> Double -> Double +gammaPdf = liftDouble3 gsl_ran_gamma_pdf + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_ran_gamma_pdf :: CDouble -> CDouble -> CDouble -> CDouble + +-- | @getGamma r a b@ gets a random gamma with parameters @a@ and @b@. +-- Uses the Marsagli-Tsang fast gamma method. +getGamma :: RNG -> Double -> Double -> IO Double +getGamma = liftRan2 gsl_ran_gamma + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_ran_gamma :: Ptr () -> CDouble -> CDouble -> IO CDouble + +-- | @getGammaKnuth r a b@ gets a random gamma with parameters @a@ and @b@. +-- Uses the algorithms from Knuth (vol 2). +getGammaKnuth :: RNG -> Double -> Double -> IO Double +getGammaKnuth = liftRan2 gsl_ran_gamma_knuth + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_ran_gamma_knuth :: Ptr () -> CDouble -> CDouble -> IO CDouble + +gammaP :: Double -> Double -> Double -> Double +gammaP = liftDouble3 gsl_cdf_gamma_P + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_gamma_P :: CDouble -> CDouble -> CDouble -> CDouble + +gammaQ :: Double -> Double -> Double -> Double +gammaQ = liftDouble3 gsl_cdf_gamma_Q + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_gamma_Q :: CDouble -> CDouble -> CDouble -> CDouble + +gammaPInv :: Double -> Double -> Double -> Double +gammaPInv = liftDouble3 gsl_cdf_gamma_Pinv + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_gamma_Pinv :: CDouble -> CDouble -> CDouble -> CDouble + +gammaQInv :: Double -> Double -> Double -> Double +gammaQInv = liftDouble3 gsl_cdf_gamma_Qinv + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_cdf_gamma_Qinv :: CDouble -> CDouble -> CDouble -> CDouble + + + +-- | @dirichletPdf xs alphas@ evaluates the probability density +-- @p(xs)@ at @xs@ for a Dirichlet distribution with parameters +-- @alphas@, where all @alphas@ are positive (strictly greater +-- than zero). Note that @xs@ and @alphas@ should have the same +-- length. +dirichletPdf :: VS.Vector Double -- ^ @xs@ + -> VS.Vector Double -- ^ @alphas@ + -> Double +dirichletPdf = dirichletPdfHelper gsl_ran_dirichlet_pdf + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_ran_dirichlet_pdf :: CSize -> Ptr Double -> Ptr Double -> CDouble + +-- | @dirichletLnPdf xs alphas == log (dirichletPdf xs alphas)@, +-- but more efficient. +dirichletLnPdf :: VS.Vector Double -- ^ @xs@ + -> VS.Vector Double -- ^ @alphas@ + -> Double +dirichletLnPdf = dirichletPdfHelper gsl_ran_dirichlet_lnpdf + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_ran_dirichlet_lnpdf :: CSize -> Ptr Double -> Ptr Double -> CDouble + +dirichletPdfHelper :: (CSize -> Ptr Double -> Ptr Double -> CDouble) + -> VS.Vector Double -> VS.Vector Double -> Double +dirichletPdfHelper f xs alphas + | len /= len2 = error "dirichletPdf*: different lengths" + | otherwise = unsafePerformIO $ + VS.unsafeWith xs $ \xs_ptr -> + VS.unsafeWith alphas $ \alphas_ptr -> + return $ realToFrac $ f lenS alphas_ptr xs_ptr + where + len = VS.length xs + len2 = VS.length alphas + lenS = fromIntegral len + +-- | @getDirichlet r alphas@ gets a random sample from a +-- Dirichlet distribution with parameters @alphas@, where all +-- @alphas@ are positive. +getDirichlet :: RNG -> VS.Vector Double -> IO (VS.Vector Double) +getDirichlet (MkRNG rng_fptr) alphas = + withForeignPtr rng_fptr $ \rng_ptr -> + VS.unsafeWith alphas $ \alphas_ptr -> do + let len = VS.length alphas + ret_fptr <- mallocForeignPtrArray len + withForeignPtr ret_fptr $ \ret_ptr -> + gsl_ran_dirichlet rng_ptr (fromIntegral len) alphas_ptr ret_ptr + return (VS.unsafeFromForeignPtr ret_fptr 0 len) + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_ran_dirichlet :: Ptr () -> CSize -> Ptr Double -> Ptr Double -> IO () + + + + + +-- Helper functions + +liftDouble :: (CDouble -> CDouble) -> Double -> Double liftDouble f x = realToFrac $ f (realToFrac x) -liftDouble2 :: (CDouble -> CDouble -> CDouble) +liftDouble2 :: (CDouble -> CDouble -> CDouble) -> Double -> Double -> Double liftDouble2 f x y = realToFrac $ f (realToFrac x) (realToFrac y) -liftDouble3 :: (CDouble -> CDouble -> CDouble -> CDouble) +liftDouble3 :: (CDouble -> CDouble -> CDouble -> CDouble) -> Double -> Double -> Double -> Double liftDouble3 f x y z = realToFrac $ f (realToFrac x) (realToFrac y) (realToFrac z) + + +liftRan0 :: (Ptr () -> IO CDouble) -> RNG -> IO Double +liftRan0 ran_fn (MkRNG fptr) = + withForeignPtr fptr $ \ptr -> + realToFrac <$> ran_fn ptr + +liftRan1 :: (Ptr () -> CDouble -> IO CDouble) + -> RNG -> Double -> IO Double +liftRan1 ran_fn (MkRNG fptr) p = + withForeignPtr fptr $ \ptr -> + realToFrac <$> ran_fn ptr (realToFrac p) + +liftRan2 :: (Ptr () -> CDouble -> CDouble -> IO CDouble) + -> RNG -> Double -> Double -> IO Double +liftRan2 ran_fn (MkRNG fptr) p q = + withForeignPtr fptr $ \ptr -> + realToFrac <$> ran_fn ptr (realToFrac p) (realToFrac q) + +liftRan3 :: (Ptr () -> CDouble -> CDouble -> CDouble -> IO CDouble) + -> RNG -> Double -> Double -> Double -> IO Double +liftRan3 ran_fn (MkRNG fptr) p q r = + withForeignPtr fptr $ \ptr -> + realToFrac <$> ran_fn ptr (realToFrac p) (realToFrac q) (realToFrac r)