hmatrix-gsl-stats 0.2.1 → 0.3
raw patch · 5 files changed
+300/−13 lines, 5 files
Files
- CHANGES +3/−0
- hmatrix-gsl-stats.cabal +1/−1
- lib/Numeric/GSL/Distribution/Continuous.hs +74/−8
- lib/Numeric/GSL/Distribution/Discrete.hs +42/−4
- lib/Numeric/GSL/distribution-aux.c +180/−0
CHANGES view
@@ -73,3 +73,6 @@ 0.2.1: Update build configuration : copied from hmatrix++0.3:+ Added vectors of samples
hmatrix-gsl-stats.cabal view
@@ -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
lib/Numeric/GSL/Distribution/Continuous.hs view
@@ -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 <at> gmail <dot> 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 -----------------------------------------------------------------------------
lib/Numeric/GSL/Distribution/Discrete.hs view
@@ -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 <at> gmail <dot> 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 -----------------------------------------------------------------------------
lib/Numeric/GSL/distribution-aux.c view
@@ -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) {