diff --git a/CHANGES b/CHANGES --- a/CHANGES +++ b/CHANGES @@ -73,3 +73,6 @@ 0.2.1: Update build configuration : copied from hmatrix + +0.3: + Added vectors of samples 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.2.1 +Version: 0.3 License: BSD3 License-file: LICENSE Copyright: (c) A.V.H. McPhail 2010, 2011, 2013, 2015 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 @@ -1,7 +1,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.GSL.Distribution.Continuous --- Copyright : (c) A. V. H. McPhail 2010 +-- Copyright : (c) A. V. H. McPhail 2010, 2015 -- License : BSD3 -- -- Maintainer : haskell.vivian.mcphail gmail com @@ -20,12 +20,12 @@ , MultiParamDist(..) , BivariateDist(..) , DistFunc(..) - , random_0p, density_0p - , random_1p, density_1p - , random_2p, density_2p - , random_3p, density_3p + , 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, density_biv + , random_biv, random_biv_v, density_biv , spherical_vector ) where @@ -34,6 +34,7 @@ import Data.Packed.Vector --import Data.Packed.Matrix hiding(toLists) import Data.Packed.Development +--import Data.Packed.Development(createVector,vec,app1,app2,app3,app4) --import Numeric.LinearAlgebra.Linear @@ -108,6 +109,16 @@ -> Double -- ^ result random_0p d s = unsafePerformIO $ distribution_random_zero_param (fromIntegral s) (fromei d) +-- | draw samples from a zero parameter distribution +random_0p_v :: ZeroParamDist -- ^ distribution type + -> Int -- ^ random seed + -> Int -- ^ number of samples + -> Vector Double -- ^ result +random_0p_v d s l = unsafePerformIO $ do + r <- createVector l + app1 (distribution_random_zero_param_v (fromIntegral s) (fromei d)) vec r "random_0p_v" + return r + -- | probability of a variate take a value outside the argument density_0p :: ZeroParamDist -- ^ density type -> DistFunc -- ^ distribution function type @@ -121,6 +132,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_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 ----------------------------------------------------------------------------- @@ -136,15 +148,27 @@ r' <- peek r return r' +-- | draw samples from a one parameter distribution +random_1p_v :: OneParamDist -- ^ distribution type + -> Int -- ^ random seed + -> Double -- ^ parameter + -> Int -- ^ number of samples + -> Vector Double -- ^ result +random_1p_v d s p l = unsafePerformIO $ do + r <- createVector l + app1 (distribution_random_one_param_v (fromIntegral s) (fromei d) p) vec r "random_1p_v" + return r + -- | probability of a variate take a value outside the argument density_1p :: OneParamDist -- ^ density type -> DistFunc -- ^ distribution function type -> Double -- ^ parameter -> Double -- ^ value - -> Double -- ^ result + -> Double -- ^ result 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_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 ----------------------------------------------------------------------------- @@ -161,6 +185,18 @@ r' <- peek r return r' +-- | draw samples from a two parameter distribution +random_2p_v :: TwoParamDist -- ^ distribution type + -> Int -- ^ random seed + -> Double -- ^ parameter 1 + -> Double -- ^ parameter 2 + -> Int -- ^ number of samples + -> Vector Double -- ^ result +random_2p_v d s p1 p2 l = unsafePerformIO $ do + r <- createVector l + app1 (distribution_random_two_param_v (fromIntegral s) (fromei d) p1 p2) vec r "random_2p_v" + return r + -- | probability of a variate take a value outside the argument density_2p :: TwoParamDist -- ^ density type -> DistFunc -- ^ distribution function type @@ -183,6 +219,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_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 ----------------------------------------------------------------------------- @@ -200,6 +237,19 @@ r' <- peek r return r' +-- | draw samples from a three parameter distribution +random_3p_v :: ThreeParamDist -- ^ distribution type + -> Int -- ^ random seed + -> Double -- ^ parameter 1 + -> Double -- ^ parameter 2 + -> Double -- ^ parameter 3 + -> Int -- ^ number of samples + -> Vector Double -- ^ result +random_3p_v d s p1 p2 p3 l = unsafePerformIO $ do + r <- createVector l + app1 (distribution_random_three_param_v (fromIntegral s) (fromei d) p1 p2 p3) vec r "random_3p_v" + return r + -- | probability of a variate take a value outside the argument density_3p :: ThreeParamDist -- ^ density type -> DistFunc -- ^ distribution function type @@ -214,11 +264,12 @@ _ -> 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_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 ----------------------------------------------------------------------------- --- | draw a sample from a three parameter distribution +-- | draw a sample from a multi parameter distribution random_mp :: MultiParamDist -- ^ distribution type -> Int -- ^ random seed -> Vector Double -- ^ parameters @@ -264,6 +315,20 @@ 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 + -> Double -- ^ parameter 2 + -> Double -- ^ parameter 3 + -> Int -- ^ number of samples + -> (Vector Double,Vector Double) -- ^ result +random_biv_v d s p1 p2 p3 l = unsafePerformIO $ do + r1 <- createVector l + r2 <- createVector l + app2 (distribution_random_bivariate_v (fromIntegral s) (fromei d) p1 p2 p3) vec r1 vec r2 "random_biv_v" + return (r1,r2) + -- | probability of a variate take a value outside the argument density_biv :: BivariateDist -- ^ density type -> DistFunc -- ^ distribution function type @@ -280,6 +345,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_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 @@ -1,7 +1,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.GSL.Distribution.Discrete --- Copyright : (c) A. V. H. McPhail 2010 +-- Copyright : (c) A. V. H. McPhail 2010, 2015 -- License : BSD3 -- -- Maintainer : haskell.vivian.mcphail gmail com @@ -19,9 +19,9 @@ , TwoParamDist(..), ThreeParamDist(..) , MultiParamDist(..) , DistFunc(..) - , random_1p, density_1p - , random_2p, density_2p - , random_3p, density_3p + , 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 ) where @@ -89,6 +89,17 @@ r' <- peek r return $ fromIntegral r' +-- | draw samples from a one parameter distribution +random_1p_v :: OneParamDist -- ^ distribution type + -> Int -- ^ random seed + -> Double -- ^ parameter + -> Int -- ^ number of samples + -> Vector Int -- ^ result +random_1p_v d s p l = unsafePerformIO $ do + r <- createVector l + app1 (distribution_discrete_one_param_v (fromIntegral s) (fromei d) p) vec r "random_1p_v" + return (mapVector fromIntegral r) + -- | probability of a variate take a value outside the argument density_1p :: OneParamDist -- ^ density type -> DistFunc -- ^ distribution function type @@ -110,6 +121,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_v" distribution_discrete_one_param_v :: CInt -> CInt -> Double -> CInt -> Ptr Int -> IO CInt foreign import ccall "distribution-aux.h discrete1_dist" distribution_dist_one_param :: CInt -> CInt -> CUInt -> Double -> IO Double ----------------------------------------------------------------------------- @@ -126,6 +138,18 @@ r' <- peek r return $ fromIntegral r' +-- | draw samples from a two parameter distribution +random_2p_v :: TwoParamDist -- ^ distribution type + -> Int -- ^ random seed + -> Double -- ^ parameter 1 + -> Int -- ^ parameter 2 + -> Int -- ^ number of samples + -> Vector Int -- ^ result +random_2p_v d s p1 p2 l = unsafePerformIO $ do + r <- createVector l + app1 (distribution_discrete_two_param_v (fromIntegral s) (fromei d) p1 (fromIntegral p2)) vec r "random_2p_v" + return (mapVector fromIntegral r) + -- | probability of a variate take a value outside the argument density_2p :: TwoParamDist -- ^ density type -> DistFunc -- ^ distribution function type @@ -141,6 +165,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_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 ----------------------------------------------------------------------------- @@ -157,6 +182,18 @@ 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 samples from a three parameter distribution +random_3p_v :: ThreeParamDist -- ^ distribution type + -> Int -- ^ random seed + -> Int -- ^ parameter 1 + -> Int -- ^ parameter 2 + -> Int -- ^ parameter 3 + -> Int -- ^ number of samples + -> Vector Int -- ^ result +random_3p_v d s p1 p2 p3 l = unsafePerformIO $ do + r <- createVector l + app1 (distribution_discrete_three_param_v (fromIntegral s) (fromei d) (fromIntegral p1) (fromIntegral p2) (fromIntegral p3)) vec r "random_3p_v" + return (mapVector fromIntegral r) -- | probability of a variate take a value outside the argument density_3p :: ThreeParamDist -- ^ density type @@ -174,6 +211,7 @@ 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_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 ----------------------------------------------------------------------------- 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 @@ -21,6 +21,25 @@ return 0; } +int random0_v(int s, int type, int rs, double* r) +{ + const gsl_rng_type * T; + gsl_rng * rng; + + gsl_rng_env_setup(); + T = gsl_rng_default; + rng = gsl_rng_alloc(T); + gsl_rng_set(rng,s); + + switch(type) { + case 0: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_landau(rng); break; } + } + + gsl_rng_free(rng); + + return 0; +} + double random0_pdf(int type, double x) { switch (type) { @@ -63,6 +82,32 @@ return 0; } +int random1_v(int s, int type, double par, int rs, double* r) +{ + const gsl_rng_type * T; + gsl_rng * rng; + + gsl_rng_env_setup(); + T = gsl_rng_default; + rng = gsl_rng_alloc(T); + gsl_rng_set(rng,s); + + switch(type) { + case 0: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_gaussian(rng,par); break; } + case 1: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_exponential(rng,par); break; } + case 2: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_laplace(rng,par); break; } + case 3: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_cauchy(rng,par); break; } + case 4: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_rayleigh(rng,par); break; } + case 5: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_chisq(rng,par); break; } + case 6: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_tdist(rng,par); break; } + case 7: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_logistic(rng,par); break; } + } + + gsl_rng_free(rng); + + return 0; +} + double random1_pdf(int type, double x, double par) { switch (type) { @@ -177,6 +222,37 @@ return 0; } +int random2_v(int s, int type, double par1, double par2, int rs, double* r) +{ + const gsl_rng_type * T; + gsl_rng * rng; + + gsl_rng_env_setup(); + T = gsl_rng_default; + rng = gsl_rng_alloc(T); + gsl_rng_set(rng,s); + + switch(type) { + case 0: { for (int i = 0; i < rs; i ++) r[i] (*r) = gsl_ran_gaussian_tail(rng,par1,par2); break; } + case 1: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_exppow(rng,par1,par2); break; } + case 2: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_rayleigh_tail(rng,par1,par2); break; } + case 3: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_levy(rng,par1,par2); break; } + case 4: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_gamma(rng,par1,par2); break; } + case 5: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_flat(rng,par1,par2); break; } + case 6: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_lognormal(rng,par1,par2); break; } + case 7: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_fdist(rng,par1,par2); break; } + case 8: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_beta(rng,par1,par2); break; } + case 9: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_pareto(rng,par1,par2); break; } + case 10: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_weibull(rng,par1,par2); break; } + case 11: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_gumbel1(rng,par1,par2); break; } + case 12: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_gumbel2(rng,par1,par2); break; } + } + + gsl_rng_free(rng); + + return 0; +} + double random2_pdf(int type, double x, double par1, double par2) { switch (type) { @@ -289,6 +365,27 @@ return 0; } +int random3_v(int s, int type, double par1, double par2, double par3, int rs, double* r) +{ + const gsl_rng_type * T; + gsl_rng * rng; + + gsl_rng_env_setup(); + T = gsl_rng_default; + rng = gsl_rng_alloc(T); + gsl_rng_set(rng,s); + + switch(type) { + case 0: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_levy_skew(rng,par1,par2,par3); break; } + } + + gsl_rng_free(rng); + + return 0; +} + +for (int i = 0; i < rs; i ++) r[i] + double random3_pdf(int type, double x, double par1, double par2, double par3) { switch (type) { @@ -356,6 +453,25 @@ return 0; } +int random_biv_v(int s, int type, double par1, double par2, double par3, int rs1, double* r1, int rs2, double* r2) +{ + const gsl_rng_type * T; + gsl_rng * rng; + + gsl_rng_env_setup(); + T = gsl_rng_default; + rng = gsl_rng_alloc(T); + gsl_rng_set(rng,s); + + switch(type) { + case 0: { for (int i = 0; i < rs; i ++) gsl_ran_bivariate_gaussian(rng,par1,par2,par3,r1[i],r2[i]); break; } + } + + gsl_rng_free(rng); + + return 0; +} + double random_biv_pdf(int type, double x, double y, double par1, double par2, double par3) { switch (type) { @@ -516,6 +632,30 @@ return 0; } +int discrete1_v(int s, int type, double par, int rs, unsigned int* r) +{ + const gsl_rng_type * T; + gsl_rng * rng; + + gsl_rng_env_setup(); + T = gsl_rng_default; + rng = gsl_rng_alloc(T); + gsl_rng_set(rng,s); + + switch(type) { + case 0: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_poisson(rng,par); break; } + case 1: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_bernoulli(rng,par); break; } + case 2: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_laplace(rng,par); break; } + case 3: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_logarithmic(rng,par); break; } + } + + gsl_rng_free(rng); + + return 0; +} + +for (int i = 0; i < rs; i ++) r[i] + double discrete1_pdf(int type, unsigned int x, double par) { switch (type) { @@ -590,6 +730,27 @@ return 0; } +int discrete2(int s, int type, double par1, unsigned int par2, int rs, unsigned int* r) +{ + const gsl_rng_type * T; + gsl_rng * rng; + + gsl_rng_env_setup(); + T = gsl_rng_default; + rng = gsl_rng_alloc(T); + gsl_rng_set(rng,s); + + switch(type) { + case 0: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_binomial(rng,par1,par2); break; } + case 1: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_negative_binomial(rng,par1,par2*1.0); break; } + case 2: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_pascal(rng,par1,par2); break; } + } + + gsl_rng_free(rng); + + return 0; +} +for (int i = 0; i < rs; i ++) r[i] double discrete2_pdf(int type, unsigned int x, double par1, unsigned int par2) { switch (type) { @@ -663,6 +824,25 @@ return 0; } +int discrete3(int s, int type, unsigned int par1, unsigned int par2, int rs, unsigned int par3, double* r) +{ + const gsl_rng_type * T; + gsl_rng * rng; + + gsl_rng_env_setup(); + T = gsl_rng_default; + rng = gsl_rng_alloc(T); + gsl_rng_set(rng,s); + + switch(type) { + case 0: { for (int i = 0; i < rs; i ++) r[i] = gsl_ran_hypergeometric(rng,par1,par2,par3); break; } + } + + gsl_rng_free(rng); + + return 0; +} +for (int i = 0; i < rs; i ++) r[i] double discrete3_pdf(int type, double x, unsigned int par1, unsigned int par2, unsigned int par3) { switch (type) {