diff --git a/CHANGES b/CHANGES --- a/CHANGES +++ b/CHANGES @@ -7,3 +7,6 @@ 0.1.1.2: fixed bug - not dividing by sum - prob for Histogram2D + +0.1.1.3: + minor improvements 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.1.1.2 +Version: 0.1.1.3 License: GPL License-file: LICENSE Copyright: (c) A.V.H. McPhail 2010 diff --git a/lib/Numeric/GSL/Histogram.hs b/lib/Numeric/GSL/Histogram.hs --- a/lib/Numeric/GSL/Histogram.hs +++ b/lib/Numeric/GSL/Histogram.hs @@ -21,7 +21,7 @@ , getBin, getRange , getMax, getMin, getBins , find - , prob + , count, prob , maxVal, maxBin, minVal, minBin , mean, stddev, sum , equalBins @@ -40,6 +40,8 @@ import Data.Packed.Development --import Numeric.LinearAlgebra.Linear +import Numeric.LinearAlgebra.Algorithms hiding(add,multiply,divide,scale) +import Numeric.LinearAlgebra.Interface() --import Control.Monad @@ -278,12 +280,14 @@ foreign import ccall "gsl-histogram.h gsl_histogram_find" histogram_find :: HistHandle -> Double -> Ptr CInt -> IO CInt +-- | find the number of occurences for each element of the input vector +count :: Histogram -> Vector Double -> Vector Double +count h = mapVector (\x -> let Just x' = find h x + in getBin h x') + -- | find the probability of occurring for each element of the input vector prob :: Histogram -> Vector Double -> Vector Double -prob h = mapVector (\x -> let s = sum h - Just x' = find h x - b = getBin h x' - in b/s) +prob h x = (count h x) / (scalar $ sum h) ----------------------------------------------------------------------------- @@ -294,8 +298,8 @@ -- | the index of the bin containing the maximum value maxBin :: Histogram -> Int maxBin (H _ h) = unsafePerformIO $ do - i <- withForeignPtr h histogram_max_bin - return $ fromIntegral i + j <- withForeignPtr h histogram_max_bin + return $ fromIntegral j -- | the minimum value contained in the bins minVal :: Histogram -> Double @@ -304,8 +308,8 @@ -- | the index of the bin containing the minimum value minBin :: Histogram -> Int minBin (H _ h) = unsafePerformIO $ do - i <- withForeignPtr h histogram_min_bin - return $ fromIntegral i + j <- withForeignPtr h histogram_min_bin + return $ fromIntegral j -- | the mean of the values, accuracy limited by bin width mean :: Histogram -> Double @@ -332,9 +336,9 @@ -- | returns True of all the individual bin ranges of the two histograms are identical equalBins :: Histogram -> Histogram -> Bool equalBins (H _ h1) (H _ h2) = unsafePerformIO $ do - i <- withForeignPtr h1 $ \p1 -> do + j <- withForeignPtr h1 $ \p1 -> do withForeignPtr h2 $ \p2 -> histogram_equal_bins p1 p2 - if (fromIntegral i) == (1 :: Int) + if (fromIntegral j) == (1 :: Int) then return True else return False diff --git a/lib/Numeric/GSL/Histogram2D.hs b/lib/Numeric/GSL/Histogram2D.hs --- a/lib/Numeric/GSL/Histogram2D.hs +++ b/lib/Numeric/GSL/Histogram2D.hs @@ -22,7 +22,7 @@ , getXMax, getYMax, getXMin, getYMin, getXBins, getYBins , reset , find - , prob + , count, prob , maxVal, maxBin, minVal, minBin , xmean, ymean, xstddev, ystddev, covariance, sum , equalBins @@ -40,6 +40,7 @@ import Data.Packed.Matrix import Data.Packed.Development +import Numeric.LinearAlgebra.Algorithms hiding(add,multiply,divide,scale) import Numeric.LinearAlgebra.Interface() --import Control.Monad @@ -320,14 +321,18 @@ foreign import ccall "gsl-histogram2d.h gsl_histogram2d_find" histogram2d_find :: Hist2DHandle -> Double -> Double -> Ptr CInt -> Ptr CInt -> IO CInt --- | find the probability of occuring for each element of the input vector -prob :: Histogram2D -> (Vector Double,Vector Double) -> Vector Double -prob hi@(H _ _ h) (x,y) = unsafePerformIO $ do +-- | 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 r <- createVector $ dim x - app3 (\xs' x' ys' y' rs' r' -> withForeignPtr h $ \h' -> histogram2d_prob h' xs' x' ys' y' rs' r') vec x vec y vec r "histogram2d_prob" - return $ r/(constant (sum hi) (dim x)) + app3 (\xs' x' ys' y' rs' r' -> withForeignPtr h $ \h' -> histogram2d_count h' xs' x' ys' y' rs' r') vec x vec y vec r "histogram2d_count" + return r -foreign import ccall "histogram-aux.h hist2d_prob" histogram2d_prob :: Hist2DHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt +foreign import ccall "histogram-aux.h hist2d_count" histogram2d_count :: Hist2DHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt + +-- | find the joint probability of occuring for each element of the input vector pair +prob :: Histogram2D -> (Vector Double,Vector Double) -> Vector Double +prob h z = (count h z) / (scalar $ sum h) ----------------------------------------------------------------------------- diff --git a/lib/Numeric/GSL/histogram-aux.c b/lib/Numeric/GSL/histogram-aux.c --- a/lib/Numeric/GSL/histogram-aux.c +++ b/lib/Numeric/GSL/histogram-aux.c @@ -109,18 +109,20 @@ return 0; } -int hist2d_prob(gsl_histogram2d* H, int xs, const double* x, int ys, const double* y, int rs, double* res) +int hist2d_count(gsl_histogram2d* H, int xs, const double* x, int ys, const double* y, int rs, double* res) { if (xs != ys) return 2000; // BAD_SIZE if (xs != rs) return 2000; // BAD_SIZE int i; - int r,c; + size_t r,c; int err; + for (i = 0; i < xs; i++) { err = gsl_histogram2d_find(H,x[i],y[i],&r,&c); if (err != 0) return 0; else { res[i] = gsl_histogram2d_get(H,r,c); } } + return 0; }