hmatrix-gsl-stats 0.1.2.1 → 0.1.2.2
raw patch · 12 files changed
+207/−45 lines, 12 files
Files
- CHANGES +4/−1
- hmatrix-gsl-stats.cabal +8/−2
- lib/Numeric/GSL/Distribution/Common.hs +3/−1
- lib/Numeric/GSL/Distribution/Continuous.hs +6/−6
- lib/Numeric/GSL/Distribution/Discrete.hs +7/−5
- lib/Numeric/GSL/Fitting/Linear.hs +3/−1
- lib/Numeric/GSL/Histogram.hs +43/−15
- lib/Numeric/GSL/Histogram2D.hs +74/−11
- lib/Numeric/GSL/Permutation.hs +3/−1
- lib/Numeric/GSL/Sort.hs +3/−1
- lib/Numeric/GSL/Statistics.hs +3/−1
- lib/Numeric/GSL/histogram-aux.c +50/−0
CHANGES view
@@ -26,4 +26,7 @@ fixed Histogram(2D)PDF sampling 0.1.2.2:- + removed warning causes in Distribution.Discrete/Continuous+ added fromVectors, modified fromMatrix in Histogram(2D) + added countInstance, probability to Histogram(2D)+
hmatrix-gsl-stats.cabal view
@@ -1,5 +1,5 @@ Name: hmatrix-gsl-stats-Version: 0.1.2.1+Version: 0.1.2.2 License: GPL License-file: LICENSE Copyright: (c) A.V.H. McPhail 2010@@ -8,7 +8,13 @@ Stability: provisional Homepage: http://code.haskell.org/hmatrix-gsl-stats Synopsis: GSL Statistics interface-Description: Purely functional interface for statistics based on hmatrix and GSL+Description: + Purely functional interface for statistics based on hmatrix and GSL+ .+ Relieves the burden of alloc/free routines and is otherwise essentially 1:1+ .+ Documentation can be found at <http://www.gnu.org/software/gsl/manual/>+ Category: Math, Statistics tested-with: GHC ==6.12.1
lib/Numeric/GSL/Distribution/Common.hs view
@@ -2,7 +2,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.GSL.Distribution.Common--- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- Copyright : (c) A. V. H. McPhail 2010 -- License : GPL-style -- -- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com@@ -10,6 +10,8 @@ -- Portability : uses ffi -- -- GSL common data types for distributions+--+-- <http://www.gnu.org/software/gsl/manual/> -- -----------------------------------------------------------------------------
lib/Numeric/GSL/Distribution/Continuous.hs view
@@ -2,7 +2,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.GSL.Distribution.Continuous--- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- Copyright : (c) A. V. H. McPhail 2010 -- License : GPL-style -- -- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com@@ -11,6 +11,8 @@ -- -- GSL continuous random distribution functions --+-- <http://www.gnu.org/software/gsl/manual/>+-- ----------------------------------------------------------------------------- module Numeric.GSL.Distribution.Continuous (@@ -205,12 +207,10 @@ -> Double -- ^ parameter 3 -> Double -- ^ value -> Double -- ^ result-density_3p d f p1 p2 p3 x = unsafePerformIO $ do+density_3p d _ _ _ _ _ = unsafePerformIO $ do case d of LevySkew -> error "Levy skew has no PDF or CDF"- where density_only f' d' p1' p2' p3' x' = if f' /= Density- then error "distribution has no CDF"- else distribution_dist_three_param (fromei f') (fromei d') x' p1' p2' p3'+ _ -> 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_dist" distribution_dist_three_param :: CInt -> CInt -> Double -> Double -> Double -> Double -> IO Double@@ -239,7 +239,7 @@ where density_only f' d' p' q' = if f' /= Density then error "distribution has no CDF" else alloca $ \r -> do- app2 (distribution_dist_multi_param (fromei f') (fromei d') r) vec p vec q "density_mp"+ app2 (distribution_dist_multi_param (fromei f') (fromei d') r) vec p' vec q' "density_mp" r' <- peek r return r'
lib/Numeric/GSL/Distribution/Discrete.hs view
@@ -2,7 +2,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.GSL.Distribution.Discrete--- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- Copyright : (c) A. V. H. McPhail 2010 -- License : GPL-style -- -- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com@@ -11,6 +11,8 @@ -- -- GSL discrete random distribution functions --+-- <http://www.gnu.org/software/gsl/manual/>+-- ----------------------------------------------------------------------------- module Numeric.GSL.Distribution.Discrete (@@ -101,10 +103,10 @@ _ -> distribution_dist_one_param (fromei f) (fromei d) (fromIntegral x) p where pdf_only f' d' p' x' = if f' /= Density then error "no CDF"- else distribution_dist_one_param (fromei f) (fromei d) (fromIntegral x) p+ else distribution_dist_one_param (fromei f') (fromei d') (fromIntegral x') p' no_inverse f' d' p' x' = if (f' == LowInv || f' == UppInv) then error "No inverse CDF"- else distribution_dist_one_param (fromei f) (fromei d) (fromIntegral x) p+ 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_dist" distribution_dist_one_param :: CInt -> CInt -> CUInt -> Double -> IO Double@@ -168,7 +170,7 @@ HyperGeometric -> no_inverse f d p1 p2 p3 x where no_inverse f' d' p1' p2' p3' x' = if (f' == LowInv || f' == UppInv) then error "No inverse CDF"- else distribution_dist_three_param (fromei f') (fromei d') (fromIntegral x') (fromIntegral p1') (fromIntegral p1') (fromIntegral p1')+ 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_dist" distribution_dist_three_param :: CInt -> CInt -> CUInt -> CUInt -> CUInt -> CUInt -> IO Double@@ -198,7 +200,7 @@ where density_only f' d' p' q' = if f' /= Density then error "distribution has no CDF" else alloca $ \r -> do- app2 (distribution_dist_multi_param (fromei f') (fromei d') r) vec p vec (mapVector (\x -> (fromIntegral x) :: CUInt) q) "density_mp"+ app2 (distribution_dist_multi_param (fromei f') (fromei d') r) vec p' vec (mapVector (\x -> (fromIntegral x) :: CUInt) q') "density_mp" r' <- peek r return r'
lib/Numeric/GSL/Fitting/Linear.hs view
@@ -2,7 +2,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.GSL.Fitting.Linear--- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- Copyright : (c) A. V. H. McPhail 2010 -- License : GPL-style -- -- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com@@ -10,6 +10,8 @@ -- Portability : uses ffi -- -- GSL linear regression functions+--+-- <http://www.gnu.org/software/gsl/manual/> -- -----------------------------------------------------------------------------
lib/Numeric/GSL/Histogram.hs view
@@ -2,7 +2,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.GSL.Histogram--- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- Copyright : (c) A. V. H. McPhail 2010 -- License : GPL-style -- -- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com@@ -11,24 +11,34 @@ -- -- GSL histogram functions --+-- <http://www.gnu.org/software/gsl/manual/>+-- ----------------------------------------------------------------------------- module Numeric.GSL.Histogram (+ -- * Creation Histogram , emptyRanges, emptyLimits , fromRanges, fromLimits+ -- * Loading , addList, addVector, addListWeighted, addVectorWeighted- , toVectors+ -- * Marshalling+ , toVectors, fromVectors+ -- * Information , getBin, getRange , getMax, getMin, getBins+ -- * Querying , find - , count, prob+ , count, prob, countInstance, probability , maxVal, maxBin, minVal, minBin+ -- * Statistics , mean, stddev, sum , equalBins+ -- * Mathematics , add, subtract, multiply, divide, shift, scale+ -- * Files , fwriteHistogram, freadHistogram, fprintfHistogram, fscanfHistogram- --+ -- * PDF , HistogramPDF , fromHistogram , sample@@ -95,14 +105,6 @@ -} ----------------------------------------------------------------------------- -vectorToTuples = toTuples . toList- where toTuples [] = error "need a minimum of two elements"- toTuples [_] = error "need a minimum of two elements"- toTuples [x1,x2] = [(x1,x2)]- toTuples (x1:x2:xs) = (x1,x2) : (toTuples (x2:xs))-------------------------------------------------------------------------------- instance Binary Histogram where put h = do let (b,c) = toVectors h@@ -111,8 +113,7 @@ get = do b <- get c <- get- return $! addVectorWeighted (emptyRanges b) (middle b) c- where middle = fromList . map (\(x1,x2) -> (x1 + x2)/2) . vectorToTuples+ return $! fromVectors b c ----------------------------------------------------------------------------- @@ -174,8 +175,27 @@ return h -----------------------------------------------------------------------------+{-+vectorToTuples = toTuples . toList+ where toTuples [] = error "need a minimum of two elements"+ toTuples [_] = error "need a minimum of two elements"+ toTuples [x1,x2] = [(x1,x2)]+ toTuples (x1:x2:xs) = (x1,x2) : (toTuples (x2:xs)) --- | extract the ranges and bins+middle = fromList . map (\(x1,x2) -> (x1 + x2)/2) . vectorToTuples+-}+-----------------------------------------------------------------------------++-- | create a histogram from the ranges and bin weights+fromVectors :: Vector Double -> Vector Double -> Histogram+fromVectors r b = unsafePerformIO $ do+ h@(H _ h') <- fromRangesIO r+ app2 (\rs r' bs b' -> withForeignPtr h' $ \h'' -> histogram_from_vectors h'' rs r' bs b') vec r vec b "fromVectors"+ return h+ +foreign import ccall "gsl-histogram.h from_vectors" histogram_from_vectors :: HistHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt++-- | extract the ranges and bin weights toVectors :: Histogram -> (Vector Double,Vector Double) -- ^ (ranges,bins) toVectors (H b h) = unsafePerformIO $ do rs <- createVector (b+1)@@ -313,6 +333,14 @@ else return Nothing foreign import ccall "gsl-histogram.h gsl_histogram_find" histogram_find :: HistHandle -> Double -> Ptr CInt -> IO CInt++-- | find the number of occurences for the input+countInstance :: Histogram -> Double -> Double+countInstance h x = let Just x' = find h x in getBin h x'++-- | find the probability of the input+probability :: Histogram -> Double -> Double+probability h x = let Just x' = find h x in (getBin h x') / (sum h) -- | find the number of occurences for each element of the input vector count :: Histogram -> Vector Double -> Vector Double
lib/Numeric/GSL/Histogram2D.hs view
@@ -2,7 +2,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.GSL.Histogram2D--- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- Copyright : (c) A. V. H. McPhail 2010 -- License : GPL-style -- -- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com@@ -11,25 +11,35 @@ -- -- GSL 2D histogram functions --+-- <http://www.gnu.org/software/gsl/manual/>+-- ----------------------------------------------------------------------------- module Numeric.GSL.Histogram2D (+ -- * Creation Histogram2D , emptyRanges, emptyLimits , fromRanges, fromLimits+ -- * Loading , addList, addVector, addListWeighted, addVectorWeighted- , toMatrix+ -- * Marshalling+ , toMatrix, fromMatrix+ -- * Information , getBin, getXRange, getYRange , getXMax, getYMax, getXMin, getYMin, getXBins, getYBins , reset+ -- * Querying , find- , count, prob+ , count, prob, probPaired, countPaired, countInstance, probability , maxVal, maxBin, minVal, minBin+ -- * Statistics , xmean, ymean, xstddev, ystddev, covariance, sum , equalBins+ -- * Mathematics , add, subtract, multiply, divide, shift, scale+ -- * Files , fwriteHistogram2D, freadHistogram2D, fprintfHistogram2D, fscanfHistogram2D- --+ -- * PDF , Histogram2DPDF , fromHistogram2D , sample@@ -50,6 +60,8 @@ import Data.Binary import Foreign hiding(shift)+--import Foreign.Storable+--import Foreign.Ptr --import Foreign.ForeignPtr --import Foreign.Marshal.Alloc(alloca) import Foreign.C.Types(CInt,CChar)@@ -65,6 +77,22 @@ ----------------------------------------------------------------------------- +-- This Storable instance for Complex is binary compatible with C99, C++ and+-- Fortran complex data types.+instance (Storable a) => Storable (a,a) where+ sizeOf z = 2 * sizeOf (fst z)+ alignment z = alignment (fst z)+ peek p = do let q = castPtr p+ r <- peek q+ i <- peekElemOff q 1+ return (r,i)+ poke p (r,i) = do let q = (castPtr p)+ poke q r+ pokeElemOff q 1 i+++-----------------------------------------------------------------------------+ data Hist2D type Hist2DHandle = Ptr Hist2D -- | A histogram structure@@ -190,13 +218,13 @@ foreign import ccall "histogram-aux.h to_matrix" histogram_to_matrix :: Hist2DHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt ------------------------------------------------------------------------------+{- vectorToTuples = toTuples . toList where toTuples [] = error "need a minimum of two elements" toTuples [_] = error "need a minimum of two elements" toTuples [x1,x2] = [(x1,x2)] toTuples (x1:x2:xs) = (x1,x2) : (toTuples (x2:xs))-+-} ----------------------------------------------------------------------------- -- | create from ranges and bins@@ -204,12 +232,13 @@ -> Vector Double -- ^ y ranges -> Matrix Double -- ^ bins -> Histogram2D -- ^result-fromMatrix x y w = let x' = map (\(x1,x2) -> (x1 + x2)/2) $ vectorToTuples x- y' = map (\(x1,x2) -> (x1 + x2)/2) $ vectorToTuples y- w' = toList $ flatten w- xy = concat $ map (\j -> zip x' (replicate (length x') j)) y'- in addListWeighted (emptyRanges x y) $ zipWith (\(a,b) d -> (a,b,d)) xy w'+fromMatrix x y w = unsafePerformIO $ do+ h@(H _ _ h') <- fromRangesIO x y+ app3 (\xs x' ys y' rs cs b -> withForeignPtr h' $ \h'' -> histogram_from_matrix h'' xs x' ys y' rs cs b) vec x vec y mat w "fromMatrix"+ return h +foreign import ccall "histogram-aux.h from_matrix" histogram_from_matrix :: Hist2DHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt+ ----------------------------------------------------------------------------- -- | create a copy of a histogram@@ -372,6 +401,14 @@ foreign import ccall "gsl-histogram2d.h gsl_histogram2d_find" histogram2d_find :: Hist2DHandle -> Double -> Double -> Ptr CInt -> Ptr CInt -> IO CInt +-- | find the number of occurences for the input+countInstance :: Histogram2D -> (Double,Double) -> Double+countInstance h x = let Just x' = find h x in getBin h x'++-- | find the probability of the input+probability :: Histogram2D -> (Double,Double) -> Double+probability h x = let Just x' = find h x in (getBin h x') / (sum h)+ -- | find the number of occurences for each element of the input vector count :: Histogram2D -> (Vector Double, Vector Double) -> Vector Double count (H _ _ h) (x,y) = unsafePerformIO $ do@@ -385,6 +422,21 @@ prob :: Histogram2D -> (Vector Double,Vector Double) -> Vector Double prob h z = (count h z) / (scalar $ sum h) +-- | find the number of occurences for each element of the input vector+countPaired :: Histogram2D -> Vector (Double,Double) -> Vector Double+countPaired (H _ _ h) x = unsafePerformIO $ do+ r <- createVector $ dim x+ app2 (\xs' x' rs' r' -> withForeignPtr h $ \h' -> histogram2d_count_pair h' xs' x' rs' r') vec x vec r "histogram2d_count_pair"+ return r++foreign import ccall "histogram-aux.h hist2d_count_pair" histogram2d_count_pair :: Hist2DHandle -> CInt -> Ptr (Double,Double) -> CInt -> Ptr Double -> IO CInt++-- | find the joint probability of occuring for each element of the input vector pair+probPaired :: Histogram2D -> (Vector (Double,Double)) -> Vector Double+probPaired h z = let l = unzip $ toList z+ z' = (fromList $ fst l,fromList $ snd l)+ in prob h z'+ ----------------------------------------------------------------------------- -- | the maximum value contained in the bins@@ -597,4 +649,15 @@ ----------------------------------------------------------------------------- -----------------------------------------------------------------------------+-----------------------------------------------------------------------------++unzipPair :: Vector (Double,Double) -> (Vector Double,Vector Double)+unzipPair v = unsafePerformIO $ do+ a <- createVector $ dim v+ b <- createVector $ dim v+ app3 histogram2d_unzip_double_pair vec v vec a vec b "unzipPair"+ return (a,b)++foreign import ccall "histogram-aux.h unzip_double_pair" histogram2d_unzip_double_pair :: CInt -> Ptr (Double,Double) -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt+ -----------------------------------------------------------------------------
lib/Numeric/GSL/Permutation.hs view
@@ -2,7 +2,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.GSL.Permutation--- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- Copyright : (c) A. V. H. McPhail 2010 -- License : GPL-style -- -- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com@@ -10,6 +10,8 @@ -- Portability : uses ffi -- -- GSL permutation functions+--+-- <http://www.gnu.org/software/gsl/manual/> -- -----------------------------------------------------------------------------
lib/Numeric/GSL/Sort.hs view
@@ -2,7 +2,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.GSL.Sort--- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- Copyright : (c) A. V. H. McPhail 2010 -- License : GPL-style -- -- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com@@ -10,6 +10,8 @@ -- Portability : uses ffi -- -- GSL sorting functions+--+-- <http://www.gnu.org/software/gsl/manual/> -- -----------------------------------------------------------------------------
lib/Numeric/GSL/Statistics.hs view
@@ -2,7 +2,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.GSL.Statistics--- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- Copyright : (c) A. V. H. McPhail 2010 -- License : GPL-style -- -- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com@@ -10,6 +10,8 @@ -- Portability : uses ffi -- -- GSL statistics functions+--+-- <http://www.gnu.org/software/gsl/manual/> -- -----------------------------------------------------------------------------
lib/Numeric/GSL/histogram-aux.c view
@@ -6,6 +6,14 @@ #include <stdio.h> +int from_vectors(gsl_histogram * H, int rs, const double* r, int bs, const double* b)+{+ if (rs != bs + 1) return 2000; // BAD_SIZE+ int i;+ for (i = 0; i < rs - 1; i++)+ gsl_histogram_accumulate(H,(r[i]+r[i+1])/2,b[i]);+ return 0;+} int to_vectors(gsl_histogram * H, int rs, double* r, int bs, double* b) {@@ -76,6 +84,18 @@ ////////////////////// +int from_matrix(gsl_histogram2d * H, int rxs, const double* rx, int rys, const double* ry, int bx, int by, const double* b)+{+ if (rxs != bx + 1) return 2000; // BAD_SIZE+ if (rys != by + 1) return 2000; // BAD_SIZE++ int i, j;+ for (i = 0; i < rxs - 1; i++)+ for (j = 0; j <rys - 1; j++)+ gsl_histogram2d_accumulate(H,(rx[i]+rx[i+1])/2,(ry[j]+ry[j+1])/2,b[i*by+j]);+ return 0;+}+ int to_matrix(gsl_histogram2d * H, int rxs, double* rx, int rys, double* ry, int bx, int by, double* b) { int bz = (rxs-1)*(rys-1); @@ -126,6 +146,22 @@ return 0; } +int hist2d_count_pair(gsl_histogram2d* H, int xs, const double* x, int rs, double* res)+{+ if (xs != rs) return 2000; // BAD_SIZE+ int i;+ size_t r,c;+ int err;++ for (i = 0; i < xs; i++) {+ err = gsl_histogram2d_find(H,x[i*2],x[i*2+1],&r,&c);+ if (err != 0) return 0;+ else { res[i] = gsl_histogram2d_get(H,r,c); }+ }++ return 0;+}+ int hist2d_fwrite(const char* filename, const gsl_histogram2d* h) { int err;@@ -166,3 +202,17 @@ return err; } +///////////////////////////////////////////////////////////////++int unzip_double_pair(int vs, const double* v, int as, double* a, int bs, double* b)+{+ if (vs != as || vs != bs) return 2000; // BAD_SIZE++ int i;+ for (i = 0; i < vs; i++) {+ a[i] = v[i*2];+ b[i] = v[i*2+1];+ }++ return 0;+}