diff --git a/CHANGES b/CHANGES --- a/CHANGES +++ b/CHANGES @@ -88,3 +88,6 @@ 0.3.0.3: removed inline modifier from statistics-aux.c fixed CInt -> CUInt in discrete_1p_v FFI call + +0.4: + added RNG type to pass between calls to RNG diff --git a/hmatrix-gsl-stats.cabal b/hmatrix-gsl-stats.cabal --- a/hmatrix-gsl-stats.cabal +++ b/hmatrix-gsl-stats.cabal @@ -1,5 +1,5 @@ Name: hmatrix-gsl-stats -Version: 0.3.0.3 +Version: 0.4 License: BSD3 License-file: LICENSE Copyright: (c) A.V.H. McPhail 2010, 2011, 2013, 2015 @@ -50,7 +50,7 @@ Numeric.GSL.Distribution.Discrete Numeric.GSL.Distribution.Common Numeric.GSL.Fitting.Linear - other-modules: + other-modules: Numeric.GSL.Distribution.Internal C-sources: lib/Numeric/GSL/statistics-aux.c lib/Numeric/GSL/sort-aux.c lib/Numeric/GSL/histogram-aux.c diff --git a/lib/Numeric/GSL/Distribution/Common.hs b/lib/Numeric/GSL/Distribution/Common.hs --- a/lib/Numeric/GSL/Distribution/Common.hs +++ b/lib/Numeric/GSL/Distribution/Common.hs @@ -1,7 +1,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.GSL.Distribution.Common --- Copyright : (c) A. V. H. McPhail 2010 +-- Copyright : (c) A. V. H. McPhail 2010, 2015 -- License : BSD3 -- -- Maintainer : haskell.vivian.mcphail gmail com @@ -16,10 +16,16 @@ module Numeric.GSL.Distribution.Common ( DistFunc(..) + ,RNG() + ,newRNG,seedRNG ) where ----------------------------------------------------------------------------- +import Numeric.GSL.Distribution.Internal + +----------------------------------------------------------------------------- + data DistFunc = Density -- ^ pdf | Lower -- ^ lower cdf | Upper -- ^ upper cdf @@ -28,3 +34,4 @@ deriving(Enum,Eq) ----------------------------------------------------------------------------- + diff --git a/lib/Numeric/GSL/Distribution/Continuous.hs b/lib/Numeric/GSL/Distribution/Continuous.hs --- a/lib/Numeric/GSL/Distribution/Continuous.hs +++ b/lib/Numeric/GSL/Distribution/Continuous.hs @@ -20,12 +20,13 @@ , MultiParamDist(..) , BivariateDist(..) , DistFunc(..) - , random_0p, random_0p_v, density_0p - , random_1p, random_1p_v, density_1p - , random_2p, random_2p_v, density_2p - , random_3p, random_3p_v, density_3p - , random_mp, density_mp - , random_biv, random_biv_v, density_biv + , random_0p,random_0p_s,random_0p_v,density_0p + , random_1p,random_1p_s,random_1p_v,density_1p + , random_2p,random_2p_s,random_2p_v,density_2p + , random_3p,random_3p_s,random_3p_v,density_3p + , random_mp,random_mp_s,density_mp + , random_biv,random_biv_s + , random_biv_v,density_biv , spherical_vector ) where @@ -55,6 +56,7 @@ --import Prelude hiding(reverse) import Numeric.GSL.Distribution.Common +import Numeric.GSL.Distribution.Internal import System.IO.Unsafe(unsafePerformIO) @@ -109,6 +111,12 @@ -> Double -- ^ result random_0p d s = unsafePerformIO $ distribution_random_zero_param (fromIntegral s) (fromei d) +-- | draw a sample from a zero parameter distribution +random_0p_s :: RNG -- ^ the random number generator + -> ZeroParamDist -- ^ distribution type + -> IO Double -- ^ result +random_0p_s (RNG rng) d = withForeignPtr rng $ \r -> distribution_random_zero_param_s r (fromei d) + -- | draw samples from a zero parameter distribution random_0p_v :: ZeroParamDist -- ^ distribution type -> Int -- ^ random seed @@ -132,6 +140,7 @@ else distribution_dist_zero_param (fromei f') (fromei d') x' foreign import ccall "distribution-aux.h random0" distribution_random_zero_param :: CInt -> CInt -> IO Double +foreign import ccall "distribution-aux.h random0_s" distribution_random_zero_param_s :: RNGHandle -> CInt -> IO Double foreign import ccall "distribution-aux.h random0_v" distribution_random_zero_param_v :: CInt -> CInt -> CInt -> Ptr Double -> IO CInt foreign import ccall "distribution-aux.h random0_dist" distribution_dist_zero_param :: CInt -> CInt -> Double -> IO Double @@ -148,6 +157,16 @@ r' <- peek r return r' +-- | draw a sample from a one parameter distribution +random_1p_s :: RNG -- ^ the random number generator + -> OneParamDist -- ^ distribution type + -> Double -- ^ parameter + -> IO Double -- ^ result +random_1p_s (RNG rng) d p = alloca $ \r -> do + check "random_1p_s" $ withForeignPtr rng $ \rg -> distribution_random_one_param_s rg (fromei d) p r + r' <- peek r + return r' + -- | draw samples from a one parameter distribution random_1p_v :: OneParamDist -- ^ distribution type -> Int -- ^ random seed @@ -168,6 +187,7 @@ density_1p d f p x = unsafePerformIO $ distribution_dist_one_param (fromei f) (fromei d) x p foreign import ccall "distribution-aux.h random1" distribution_random_one_param :: CInt -> CInt -> Double -> Ptr Double -> IO CInt +foreign import ccall "distribution-aux.h random1_s" distribution_random_one_param_s :: RNGHandle -> CInt -> Double -> Ptr Double -> IO CInt foreign import ccall "distribution-aux.h random1_v" distribution_random_one_param_v :: CInt -> CInt -> Double -> CInt -> Ptr Double -> IO CInt foreign import ccall "distribution-aux.h random1_dist" distribution_dist_one_param :: CInt -> CInt -> Double -> Double -> IO Double @@ -185,6 +205,17 @@ r' <- peek r return r' +-- | draw a sample from a two parameter distribution +random_2p_s :: RNG -- ^ the random number generator + -> TwoParamDist -- ^ distribution type + -> Double -- ^ parameter 1 + -> Double -- ^ parameter 2 + -> IO Double -- ^ result +random_2p_s (RNG rng) d p1 p2 = alloca $ \r -> do + check "random_2p_s" $ withForeignPtr rng $ \rg -> distribution_random_two_param_s rg (fromei d) p1 p2 r + r' <- peek r + return r' + -- | draw samples from a two parameter distribution random_2p_v :: TwoParamDist -- ^ distribution type -> Int -- ^ random seed @@ -219,6 +250,7 @@ else distribution_dist_two_param (fromei f') (fromei d') x' p1' p2' foreign import ccall "distribution-aux.h random2" distribution_random_two_param :: CInt -> CInt -> Double -> Double -> Ptr Double -> IO CInt +foreign import ccall "distribution-aux.h random2_s" distribution_random_two_param_s :: RNGHandle -> CInt -> Double -> Double -> Ptr Double -> IO CInt foreign import ccall "distribution-aux.h random2_v" distribution_random_two_param_v :: CInt -> CInt -> Double -> Double -> CInt -> Ptr Double -> IO CInt foreign import ccall "distribution-aux.h random2_dist" distribution_dist_two_param :: CInt -> CInt -> Double -> Double -> Double -> IO Double @@ -237,6 +269,18 @@ r' <- peek r return r' +-- | draw a sample from a three parameter distribution +random_3p_s :: RNG -- ^ the random number generator + -> ThreeParamDist -- ^ distribution type + -> Double -- ^ parameter 1 + -> Double -- ^ parameter 2 + -> Double -- ^ parameter 3 + -> IO Double -- ^ result +random_3p_s (RNG rng) d p1 p2 p3 = alloca $ \r -> do + check "random_3p_s" $ withForeignPtr rng $ \rg -> distribution_random_three_param_s rg (fromei d) p1 p2 p3 r + r' <- peek r + return r' + -- | draw samples from a three parameter distribution random_3p_v :: ThreeParamDist -- ^ distribution type -> Int -- ^ random seed @@ -264,6 +308,7 @@ _ -> error "unknown 3 parameter distribution" foreign import ccall "distribution-aux.h random3" distribution_random_three_param :: CInt -> CInt -> Double -> Double -> Double -> Ptr Double -> IO CInt +foreign import ccall "distribution-aux.h random3_s" distribution_random_three_param_s :: RNGHandle -> CInt -> Double -> Double -> Double -> Ptr Double -> IO CInt foreign import ccall "distribution-aux.h random3_v" distribution_random_three_param_v :: CInt -> CInt -> Double -> Double -> Double -> CInt -> Ptr Double -> IO CInt foreign import ccall "distribution-aux.h random3_dist" distribution_dist_three_param :: CInt -> CInt -> Double -> Double -> Double -> Double -> IO Double @@ -279,6 +324,17 @@ app2 (distribution_random_multi_param (fromIntegral s) (fromei d)) vec p vec r "random_mp" return r +-- | draw a sample from a multi parameter distribution +random_mp_s :: RNG -- ^ the random number generator + -> MultiParamDist -- ^ distribution type + -> Vector Double -- ^ parameters + -> IO (Vector Double) -- ^ result +random_mp_s (RNG rng) d p = do + let l = dim p + r <- createVector l + withForeignPtr rng $ \rg -> app2 (distribution_random_multi_param_s rg (fromei d)) vec p vec r "random_mp_s" + return r + -- | probability of a variate take a value outside the argument density_mp :: MultiParamDist -- ^ density type -> DistFunc -- ^ distribution function type @@ -296,6 +352,7 @@ return r' foreign import ccall "distribution-aux.h random_mp" distribution_random_multi_param :: CInt -> CInt -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt +foreign import ccall "distribution-aux.h random_mp_s" distribution_random_multi_param_s :: RNGHandle -> CInt -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt foreign import ccall "distribution-aux.h random_mp_dist" distribution_dist_multi_param :: CInt -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt ----------------------------------------------------------------------------- @@ -316,6 +373,21 @@ return (r1',r2') -- | draw a sample from a bivariate distribution +random_biv_s :: RNG -- ^ the random number generator + -> BivariateDist -- ^ distribution type + -> Double -- ^ parameter 1 + -> Double -- ^ parameter 2 + -> Double -- ^ parameter 3 + -> IO (Double,Double) -- ^ result +random_biv_s (RNG rng) d p1 p2 p3 = do + alloca $ \r1 -> + alloca $ \r2 -> do + withForeignPtr rng $ \r -> distribution_random_bivariate_s r (fromei d) p1 p2 p3 r1 r2 + r1' <- peek r1 + r2' <- peek r2 + return (r1',r2') + +-- | draw a sample from a bivariate distribution random_biv_v :: BivariateDist -- ^ distribution type -> Int -- ^ random seed -> Double -- ^ parameter 1 @@ -345,6 +417,7 @@ else distribution_dist_bivariate (fromei f') (fromei d') x' y' p1' p2' p3' foreign import ccall "distribution-aux.h random_biv" distribution_random_bivariate :: CInt -> CInt -> Double -> Double -> Double -> Ptr Double -> Ptr Double -> IO () +foreign import ccall "distribution-aux.h random_biv_s" distribution_random_bivariate_s :: RNGHandle -> CInt -> Double -> Double -> Double -> Ptr Double -> Ptr Double -> IO () foreign import ccall "distribution-aux.h random_biv_v" distribution_random_bivariate_v :: CInt -> CInt -> Double -> Double -> Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt foreign import ccall "distribution-aux.h random_biv_dist" distribution_dist_bivariate :: CInt -> CInt -> Double -> Double -> Double -> Double -> Double -> IO Double diff --git a/lib/Numeric/GSL/Distribution/Discrete.hs b/lib/Numeric/GSL/Distribution/Discrete.hs --- a/lib/Numeric/GSL/Distribution/Discrete.hs +++ b/lib/Numeric/GSL/Distribution/Discrete.hs @@ -50,6 +50,7 @@ --import Prelude hiding(reverse) import Numeric.GSL.Distribution.Common +import Numeric.GSL.Distribution.Internal import System.IO.Unsafe(unsafePerformIO) @@ -89,6 +90,18 @@ r' <- peek r return $ fromIntegral r' +-- | draw a sample from a one parameter distribution +random_1p_s :: RNG -- ^ RNG + -> OneParamDist -- ^ distribution type + -> Double -- ^ parameter + -> IO Int -- ^ result +random_1p_s (RNG rng) d p = do + alloca $ \r -> + withForeignPtr rng $ \rg -> do + check "random1p_s" $ distribution_discrete_one_param_s rg (fromei d) p r + r' <- peek r + return $ fromIntegral r' + -- | draw samples from a one parameter distribution random_1p_v :: OneParamDist -- ^ distribution type -> Int -- ^ random seed @@ -121,6 +134,7 @@ else distribution_dist_one_param (fromei f') (fromei d') (fromIntegral x') p' foreign import ccall "distribution-aux.h discrete1" distribution_discrete_one_param :: CInt -> CInt -> Double -> Ptr CUInt -> IO CInt +foreign import ccall "distribution-aux.h discrete1_s" distribution_discrete_one_param_s :: RNGHandle -> CInt -> Double -> Ptr CUInt -> IO CInt foreign import ccall "distribution-aux.h discrete1_v" distribution_discrete_one_param_v :: CInt -> CInt -> Double -> CInt -> Ptr CUInt -> IO CInt foreign import ccall "distribution-aux.h discrete1_dist" distribution_dist_one_param :: CInt -> CInt -> CUInt -> Double -> IO Double @@ -138,6 +152,18 @@ r' <- peek r return $ fromIntegral r' +-- | draw a sample from a two parameter distribution +random_2p_s :: RNG -- ^ RNG + -> TwoParamDist -- ^ distribution type + -> Double -- ^ parameter 1 + -> Int -- ^ parameter 2 + -> IO Int -- ^ result +random_2p_s (RNG rng) d p1 p2 = alloca $ \r -> + withForeignPtr rng $ \rg -> do + check "random2p_s" $ distribution_discrete_two_param_s rg (fromei d) p1 (fromIntegral p2) r + r' <- peek r + return $ fromIntegral r' + -- | draw samples from a two parameter distribution random_2p_v :: TwoParamDist -- ^ distribution type -> Int -- ^ random seed @@ -165,6 +191,7 @@ else distribution_dist_two_param (fromei f') (fromei d') (fromIntegral x') p1' (fromIntegral p2') foreign import ccall "distribution-aux.h discrete2" distribution_discrete_two_param :: CInt -> CInt -> Double -> CUInt -> Ptr CUInt -> IO CInt +foreign import ccall "distribution-aux.h discrete2_s" distribution_discrete_two_param_s :: RNGHandle -> CInt -> Double -> CUInt -> Ptr CUInt -> IO CInt foreign import ccall "distribution-aux.h discrete2_v" distribution_discrete_two_param_v :: CInt -> CInt -> Double -> CUInt -> CInt -> Ptr CUInt -> IO CInt foreign import ccall "distribution-aux.h discrete2_dist" distribution_dist_two_param :: CInt -> CInt -> CUInt -> Double -> CUInt -> IO Double @@ -182,6 +209,20 @@ check "random_3p" $ distribution_discrete_three_param (fromIntegral s) (fromei d) (fromIntegral p1) (fromIntegral p2) (fromIntegral p3) r r' <- peek r return $ fromIntegral r' + +-- | draw a sample from a three parameter distribution +random_3p_s :: RNG -- ^ RNG + -> ThreeParamDist -- ^ distribution type + -> Int -- ^ parameter 1 + -> Int -- ^ parameter 2 + -> Int -- ^ parameter 3 + -> IO Int -- ^ result +random_3p_s (RNG rng) d p1 p2 p3 = alloca $ \r -> + withForeignPtr rng $ \rg -> do + check "random_3p_s" $ distribution_discrete_three_param_s rg (fromei d) (fromIntegral p1) (fromIntegral p2) (fromIntegral p3) r + r' <- peek r + return $ fromIntegral r' + -- | draw samples from a three parameter distribution random_3p_v :: ThreeParamDist -- ^ distribution type -> Int -- ^ random seed @@ -211,12 +252,13 @@ else distribution_dist_three_param (fromei f') (fromei d') (fromIntegral x') (fromIntegral p1') (fromIntegral p2') (fromIntegral p3') foreign import ccall "distribution-aux.h discrete3" distribution_discrete_three_param :: CInt -> CInt -> CUInt -> CUInt -> CUInt -> Ptr CUInt -> IO CInt +foreign import ccall "distribution-aux.h discrete3_s" distribution_discrete_three_param_s :: RNGHandle -> CInt -> CUInt -> CUInt -> CUInt -> Ptr CUInt -> IO CInt foreign import ccall "distribution-aux.h discrete3_v" distribution_discrete_three_param_v :: CInt -> CInt -> CUInt -> CUInt -> CUInt -> CInt -> Ptr CUInt -> IO CInt foreign import ccall "distribution-aux.h discrete3_dist" distribution_dist_three_param :: CInt -> CInt -> CUInt -> CUInt -> CUInt -> CUInt -> IO Double ----------------------------------------------------------------------------- --- | draw a sample from a three parameter distribution +-- | draw a sample from a multi-parameter distribution random_mp :: MultiParamDist -- ^ distribution type -> Int -- ^ random seed -> Int -- ^ trials @@ -227,6 +269,17 @@ app2 (distribution_discrete_multi_param (fromIntegral s) (fromei d) (fromIntegral t)) vec p vec r "random_mp" return $ mapVector (\x -> (fromIntegral x) :: Int) r +-- | draw a sample from a multi-parameter distribution +random_mp_s :: RNG -- ^ RNG + -> MultiParamDist -- ^ distribution type + -> Int -- ^ trials + -> Vector Double -- ^ parameters + -> IO (Vector Int) -- ^ result +random_mp_s (RNG rng) d t p = withForeignPtr rng $ \rg -> do + r <- createVector $ dim p + app2 (distribution_discrete_multi_param_s rg (fromei d) (fromIntegral t)) vec p vec r "random_mp_s" + return $ mapVector (\x -> (fromIntegral x) :: Int) r + -- | probability of a variate take a value outside the argument density_mp :: MultiParamDist -- ^ density type -> DistFunc -- ^ distribution function type @@ -244,6 +297,7 @@ return r' foreign import ccall "distribution-aux.h discrete_mp" distribution_discrete_multi_param :: CInt -> CInt -> CUInt -> CInt -> Ptr Double -> CInt -> Ptr CUInt -> IO CInt +foreign import ccall "distribution-aux.h discrete_mp_s" distribution_discrete_multi_param_s :: RNGHandle -> CInt -> CUInt -> CInt -> Ptr Double -> CInt -> Ptr CUInt -> IO CInt foreign import ccall "distribution-aux.h discrete_mp_dist" distribution_dist_multi_param :: CInt -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr CUInt -> IO CInt ----------------------------------------------------------------------------- diff --git a/lib/Numeric/GSL/Distribution/Internal.hs b/lib/Numeric/GSL/Distribution/Internal.hs new file mode 100644 --- /dev/null +++ b/lib/Numeric/GSL/Distribution/Internal.hs @@ -0,0 +1,54 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.GSL.Distribution.Internal +-- Copyright : (c) A. V. H. McPhail 2015 +-- License : BSD3 +-- +-- Maintainer : haskell.vivian.mcphail gmail com +-- Stability : provisional +-- Portability : uses ffi +-- +-- GSL common data types for distributions +-- +-- +-- +----------------------------------------------------------------------------- + +module Numeric.GSL.Distribution.Internal ( + RNGHandle + ,RNG(..) + ,newRNG,seedRNG + ) where + +----------------------------------------------------------------------------- + +import Foreign.Ptr +import Foreign.ForeignPtr +import Foreign.C.Types(CULong(..)) + +----------------------------------------------------------------------------- + +data GSL_RNG +type RNGHandle = Ptr GSL_RNG +data RNG = RNG { gsl_rng :: {-# UNPACK #-} !(ForeignPtr GSL_RNG) } + +----------------------------------------------------------------------------- + +-- | create a random number generator (RNG) object +newRNG :: IO RNG +newRNG = do + rng <- rng_new + rng' <- newForeignPtr rng_free rng + return $ RNG rng' + +-- | assign a new seed to an RNG +seedRNG :: RNG -> Int -> IO () +seedRNG (RNG rng) s = withForeignPtr rng $ \r -> rng_seed r (fromIntegral s) + +----------------------------------------------------------------------------- + +foreign import ccall "distribution-aux.h new_rng" rng_new :: IO RNGHandle +foreign import ccall "gsl_rng.h gsl_rng_set" rng_seed :: RNGHandle -> CULong -> IO () +foreign import ccall "gsl_rng.h &gsl_rng_free" rng_free :: FunPtr (RNGHandle -> IO ()) + +----------------------------------------------------------------------------- diff --git a/lib/Numeric/GSL/distribution-aux.c b/lib/Numeric/GSL/distribution-aux.c --- a/lib/Numeric/GSL/distribution-aux.c +++ b/lib/Numeric/GSL/distribution-aux.c @@ -2,6 +2,18 @@ #include #include +gsl_rng* new_rng() +{ + const gsl_rng_type * T; + gsl_rng * rng; + + gsl_rng_env_setup(); + T = gsl_rng_default; + rng = gsl_rng_alloc(T); + + return rng; +} + int random0(int s, int type, double* r) { const gsl_rng_type * T; @@ -21,6 +33,14 @@ return 0; } +int random0_s(gsl_rng* rng, int type, double* r) +{ + switch(type) { + case 0: { (*r) = gsl_ran_landau(rng); break; } + } + return 0; +} + int random0_v(int s, int type, int rs, double* r) { const gsl_rng_type * T; @@ -84,6 +104,21 @@ return 0; } +int random1_s(gsl_rng* rng, int type, double par, double* r) +{ + switch(type) { + case 0: { (*r) = gsl_ran_gaussian(rng,par); break; } + case 1: { (*r) = gsl_ran_exponential(rng,par); break; } + case 2: { (*r) = gsl_ran_laplace(rng,par); break; } + case 3: { (*r) = gsl_ran_cauchy(rng,par); break; } + case 4: { (*r) = gsl_ran_rayleigh(rng,par); break; } + case 5: { (*r) = gsl_ran_chisq(rng,par); break; } + case 6: { (*r) = gsl_ran_tdist(rng,par); break; } + case 7: { (*r) = gsl_ran_logistic(rng,par); break; } + } + return 0; +} + int random1_v(int s, int type, double par, int rs, double* r) { const gsl_rng_type * T; @@ -226,6 +261,26 @@ return 0; } +int random2_s(gsl_rng* rng, int type, double par1, double par2, double* r) +{ + switch(type) { + case 0: { (*r) = gsl_ran_gaussian_tail(rng,par1,par2); break; } + case 1: { (*r) = gsl_ran_exppow(rng,par1,par2); break; } + case 2: { (*r) = gsl_ran_rayleigh_tail(rng,par1,par2); break; } + case 3: { (*r) = gsl_ran_levy(rng,par1,par2); break; } + case 4: { (*r) = gsl_ran_gamma(rng,par1,par2); break; } + case 5: { (*r) = gsl_ran_flat(rng,par1,par2); break; } + case 6: { (*r) = gsl_ran_lognormal(rng,par1,par2); break; } + case 7: { (*r) = gsl_ran_fdist(rng,par1,par2); break; } + case 8: { (*r) = gsl_ran_beta(rng,par1,par2); break; } + case 9: { (*r) = gsl_ran_pareto(rng,par1,par2); break; } + case 10: { (*r) = gsl_ran_weibull(rng,par1,par2); break; } + case 11: { (*r) = gsl_ran_gumbel1(rng,par1,par2); break; } + case 12: { (*r) = gsl_ran_gumbel2(rng,par1,par2); break; } + } + return 0; +} + int random2_v(int s, int type, double par1, double par2, int rs, double* r) { const gsl_rng_type * T; @@ -371,6 +426,14 @@ return 0; } +int random3_s(gsl_rng* rng, int type, double par1, double par2, double par3, double* r) +{ + switch(type) { + case 0: { (*r) = gsl_ran_levy_skew(rng,par1,par2,par3); break; } + } + return 0; +} + int random3_v(int s, int type, double par1, double par2, double par3, int rs, double* r) { const gsl_rng_type * T; @@ -459,6 +522,14 @@ return 0; } +int random_biv_s(gsl_rng* rng, int type, double par1, double par2, double par3, double* r1, double* r2) +{ + switch(type) { + case 0: { gsl_ran_bivariate_gaussian(rng,par1,par2,par3,r1,r2); break; } + } + return 0; +} + int random_biv_v(int s, int type, double par1, double par2, double par3, int rs1, double* r1, int rs2, double* r2) { if (rs1 != rs2) return 2000; // BAD_SIZE @@ -551,6 +622,14 @@ return 0; } +int random_mp_s(gsl_rng* rng, int type, int ps, const double* p, int rs, double* r) +{ + switch(type) { + case 0: { gsl_ran_dirichlet(rng,ps,p,r); break; } + } + return 0; +} + double random_mp_pdf(int type, int ps, const double* p, int qs, const double* q) { switch (type) { @@ -642,6 +721,17 @@ return 0; } +int discrete1_s(gsl_rng* rng, int type, double par, unsigned int* r) +{ + switch(type) { + case 0: { (*r) = gsl_ran_poisson(rng,par); break; } + case 1: { (*r) = gsl_ran_bernoulli(rng,par); break; } + case 2: { (*r) = gsl_ran_laplace(rng,par); break; } + case 3: { (*r) = gsl_ran_logarithmic(rng,par); break; } + } + return 0; +} + int discrete1_v(int s, int type, double par, int rs, unsigned int* r) { const gsl_rng_type * T; @@ -740,6 +830,16 @@ return 0; } +int discrete2_s(gsl_rng* rng, int type, double par1, unsigned int par2, unsigned int* r) +{ + switch(type) { + case 0: { (*r) = gsl_ran_binomial(rng,par1,par2); break; } + case 1: { (*r) = gsl_ran_negative_binomial(rng,par1,par2*1.0); break; } + case 2: { (*r) = gsl_ran_pascal(rng,par1,par2); break; } + } + return 0; +} + int discrete2_v(int s, int type, double par1, unsigned int par2, int rs, unsigned int* r) { const gsl_rng_type * T; @@ -836,6 +936,15 @@ return 0; } +int discrete3_s(gsl_rng* rng, int type, unsigned int par1, unsigned int par2, unsigned int par3, double* r) +{ + switch(type) { + case 0: { (*r) = gsl_ran_hypergeometric(rng,par1,par2,par3); break; } + } + + return 0; +} + int discrete3_v(int s, int type, unsigned int par1, unsigned int par2, unsigned int par3, int rs, double* r) { const gsl_rng_type * T; @@ -923,6 +1032,14 @@ gsl_rng_free(rng); + return 0; +} + +int discrete_mp_s(gsl_rng* rng, int type, unsigned int n, int ps, const double* p, int rs, unsigned int* r) +{ + switch(type) { + case 0: { gsl_ran_multinomial(rng,ps,n,p,r); break; } + } return 0; }