statistics-dirichlet (empty) → 0.6
raw patch · 9 files changed
+1159/−0 lines, 9 filesdep +basedep +deepseqdep +hmatrix-specialsetup-changed
Dependencies added: base, deepseq, hmatrix-special, nonlinear-optimization, vector
Files
- LICENSE +30/−0
- Setup.lhs +3/−0
- src/Math/Statistics/Dirichlet.hs +47/−0
- src/Math/Statistics/Dirichlet/Density.hs +157/−0
- src/Math/Statistics/Dirichlet/Matrix.hs +217/−0
- src/Math/Statistics/Dirichlet/Mixture.hs +533/−0
- src/Math/Statistics/Dirichlet/Options.hs +91/−0
- src/Math/Statistics/Dirichlet/Util.hs +32/−0
- statistics-dirichlet.cabal +49/−0
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2009-2012, Felipe Lessa++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Felipe Lessa nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Setup.lhs view
@@ -0,0 +1,3 @@+#!/usr/bin/runhaskell+> import Distribution.Simple+> main = defaultMain
+ src/Math/Statistics/Dirichlet.hs view
@@ -0,0 +1,47 @@+---------------------------------------------------------------------------+-- | Module : Math.Statistics.Dirichlet+-- Copyright : (c) 2009-2012 Felipe Lessa+-- License : BSD3+--+-- Maintainer : felipe.lessa@gmail.com+-- Stability : experimental+-- Portability : portable+--+-- This module re-exports functions from+-- "Math.Statistics.Dirichlet.Mixture" and+-- "Math.Statistics.Dirichlet.Options" in a more digestable way.+-- Since this library is under-documented, I recommend reading+-- the documentation of the symbols re-exported here.+--+-- This module does not use "Math.Statistics.Dirichlet.Density"+-- in any way. If you don't need mixtures then you should+-- probably use that module directly since it's faster and more+-- reliable (less magic happens there).+--+--------------------------------------------------------------------------++module Math.Statistics.Dirichlet+ ( -- * Data types (re-exported)+ DirichletMixture(..)+ , empty+ , Component+ , fromList+ , toList+ -- * Options (re-exported)+ , TrainingVector+ , TrainingVectors+ , StepSize(..)+ , Delta+ , Predicate(..)+ , Reason(..)+ , Result(..)+ -- * Training data (re-exported)+ , TrainingData+ , prepareTraining+ -- * Functions (re-exported)+ , derive+ , cost+ ) where++import Math.Statistics.Dirichlet.Mixture+import Math.Statistics.Dirichlet.Options
+ src/Math/Statistics/Dirichlet/Density.hs view
@@ -0,0 +1,157 @@+---------------------------------------------------------------------------+-- | Module : Math.Statistics.Dirichlet.Density+-- Copyright : (c) 2009-2012 Felipe Lessa+-- License : BSD3+--+-- Maintainer : felipe.lessa@gmail.com+-- Stability : experimental+-- Portability : portable+--+--------------------------------------------------------------------------++module Math.Statistics.Dirichlet.Density+ ( DirichletDensity(..)+ , empty+ , fromList+ , toList+ , derive+ , cost+ ) where++import qualified Data.Vector as V+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed as U++import Control.DeepSeq (NFData(..))+import Numeric.GSL.Special.Gamma (lngamma)+import Numeric.GSL.Special.Psi (psi)++import Math.Statistics.Dirichlet.Options+import Math.Statistics.Dirichlet.Util++++-- | A Dirichlet density.+newtype DirichletDensity = DD {unDD :: U.Vector Double} deriving (Eq)++instance Show DirichletDensity where+ showsPrec prec (DD v) =+ showParen (prec > 10) $+ showString "fromList " .+ showsPrec 11 (U.toList v)++instance Read DirichletDensity where+ readsPrec p ('(':xs) = let (ys,')':zs) = break (== ')') xs+ in map (\(x,s) -> (x,s++zs)) $+ readsPrec p ys+ readsPrec p xs = let [("fromList",list)] = lex xs+ in map (\(x,s) -> (fromList x,s)) $+ readsPrec p list++instance NFData DirichletDensity where+ rnf DD {} = ()++-- | @empty n x@ is an \"empty\" Dirichlet density with size+-- @n@ and all alphas set to @x@.+empty :: Int -> Double -> DirichletDensity+empty = (DD .) . U.replicate+{-# INLINE empty #-}++-- | @fromList xs@ constructs a Dirichlet density from a list of+-- alpha values.+fromList :: [Double] -> DirichletDensity+fromList = DD . U.fromList+{-# INLINE fromList #-}++-- | @toList d@ deconstructs a Dirichlet density to a list of+-- alpha values.+toList :: DirichletDensity -> [Double]+toList (DD xs) = U.toList xs+{-# INLINE toList #-}++-- | Derive a Dirichlet density using a maximum likelihood method+-- as described by Karplus et al (equation 26). All training+-- vectors should have the same length, however this is not+-- verified.+derive :: DirichletDensity -> Predicate -> StepSize+ -> TrainingVectors -> Result DirichletDensity+derive (DD initial) (Pred maxIter' minDelta_ deltaSteps' _ _)+ (Step step) trainingData+ | V.length trainingData == 0 = err "empty training data"+ | U.length initial < 1 = err "empty initial vector"+ | maxIter' < 1 = err "non-positive maxIter"+ | minDelta_ < 0 = err "negative minDelta"+ | deltaSteps' < 1 = err "non-positive deltaSteps"+ | step <= 0 = err "non-positive step"+ | step >= 1 = err "step greater than one"+ | otherwise = train+ where+ err = error . ("Dirichlet.derive: " ++)++ -- Compensate the different deltaSteps.+ !minDelta' = minDelta_ * fromIntegral deltaSteps'++ -- Number of training sequences.+ !trainingSize = fromIntegral $ V.length trainingData++ -- Sums of each training sequence.+ trainingSums :: U.Vector Double+ !trainingSums = G.unstream $ G.stream $ V.map U.sum trainingData++ -- Functions that work on the alphas only (and not their logs).+ calcSumAs = U.sum . snd . U.unzip+ finish = DD . snd . U.unzip++ -- Start training in the zero-th iteration and with+ -- infinite inital cost.+ train = train' 1 infinity (U.sum initial) $+ U.map (\x -> (log x, x)) initial++ train' !iter !oldCost !sumAs !alphas =+ -- Reestimate alpha's.+ let !alphas' = U.imap calculateAlphas alphas+ !psiSumAs = psi sumAs+ !psiSums = U.sum $ U.map (\sumT -> psi $ sumT + sumAs) trainingSums+ calculateAlphas !i (!w, !a) =+ let !s1 = trainingSize * (psiSumAs - psi a)+ !s2 = V.sum $ V.map (\t -> psi $ t U.! i + a) trainingData+ !w' = w + step * a * (s1 + s2 - psiSums)+ !a' = exp w'+ in (w', a')++ -- Recalculate constants.+ !sumAs' = calcSumAs alphas'+ !calcCost = iter `mod` deltaSteps' == 0+ !cost' = if calcCost then newCost else oldCost+ where newCost = costWorker (snd $ U.unzip alphas') sumAs'+ trainingData trainingSums+ !delta = abs (cost' - oldCost)++ -- Verify convergence. Even with MaxIter we only stop+ -- iterating if the delta was calculated. Otherwise we+ -- wouldn't be able to tell the caller why the delta was+ -- still big when we reached the limit.+ in case (calcCost, delta <= minDelta', iter >= maxIter') of+ (True, True, _) -> Result Delta iter delta cost' $ finish alphas'+ (True, _, True) -> Result MaxIter iter delta cost' $ finish alphas'+ _ -> train' (iter+1) cost' sumAs' alphas'++-- | Cost function for deriving a Dirichlet density (equation+-- 18). This function is minimized by 'derive'.+cost :: TrainingVectors -> DirichletDensity -> Double+cost tv (DD arr) = costWorker arr (U.sum arr) tv $+ G.unstream $ G.stream $ V.map U.sum tv++-- | 'cost' needs to calculate the sum of all training vectors.+-- This functios avoids recalculting this quantity in 'derive'+-- multiple times. This is the used by both 'cost' and 'derive'.+costWorker :: U.Vector Double -> Double -> TrainingVectors -> U.Vector Double -> Double+costWorker !alphas !sumAs !trainingData !trainingSums =+ let !lngammaSumAs = lngamma sumAs+ f t = U.sum $ U.zipWith w t alphas+ where w t_i a_i = lngamma (t_i + a_i) - lngamma (t_i + 1) - lngamma a_i+ g sumT = lngamma (sumT+1) - lngamma (sumT + sumAs)+ in negate $ (V.sum $ V.map f trainingData)+ + (U.sum $ U.map g trainingSums)+ + lngammaSumAs * fromIntegral (U.length trainingSums)+{-# INLINE costWorker #-}
+ src/Math/Statistics/Dirichlet/Matrix.hs view
@@ -0,0 +1,217 @@+---------------------------------------------------------------------------+-- | Module : Math.Statistics.Dirichlet.Matrix+-- Copyright : (c) 2009-2012 Felipe Lessa+-- License : BSD3+--+-- Maintainer : felipe.lessa@gmail.com+-- Stability : experimental+-- Portability : portable+--+-- Implement matrices using plain 'U.Vector's with data stored in+-- row-major order (i.e. the first elements correspond to the+-- first row).+--+--------------------------------------------------------------------------++module Math.Statistics.Dirichlet.Matrix+ ( -- * Basic+ Matrix(..)+ , size+ , (!)+ -- * Constructing+ , replicate+ , replicateRows+ , fromVector+ , fromVectorT+ -- * Rows+ , rows+ , (!!!)+ -- * Columns+ , cols+ , col+ -- * Maps and zips+ , umap+ , map+ , imap+ , rowmap+ , irowmap+ , uzipWith+ , zipWith+ , izipWith+ , rzipWith+ -- * Other+ , transpose+ ) where++import Prelude hiding (replicate, map, zipWith)+import System.IO.Unsafe (unsafePerformIO)+import qualified Data.Vector as V+import qualified Data.Vector.Fusion.Stream as S+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Unboxed.Mutable as MU+++-- | A matrix.+data Matrix = M { mRows :: !Int+ , mCols :: !Int+ , mData :: !(U.Vector Double)}+ deriving (Eq, Ord, Show)++-- | Size of the matrix.+size :: Matrix -> (Int,Int)+size m = (mRows m, mCols m)++-- | Element at position.+(!) :: Matrix -> (Int,Int) -> Double+(!) m (r,c) = mData m U.! (r * mCols m + c)++++-- | A matrix where all elements are of the same value.+replicate :: (Int,Int) -> Double -> Matrix+replicate (r,c) v = M { mRows = r+ , mCols = c+ , mData = U.replicate (r*c) v}++-- | A matrix where all rows are of the same value.+replicateRows :: Int -> U.Vector Double -> Matrix+replicateRows r v =+ let c = U.length v+ in M { mRows = r+ , mCols = c+ , mData = U.generate (r*c) (\i -> v U.! (i `mod` c))}++-- | Creates a matrix from a vector of vectors. It *is not*+-- verified that the vectors have the right length.+fromVector :: (G.Vector v (w Double), G.Vector w Double)+ => v (w Double) -> Matrix+fromVector v =+ M { mRows = G.length v+ , mCols = G.length (G.head v)+ , mData = G.unstream $ S.concatMap G.stream $ G.stream v}++-- | Creates a matrix from a vector of vectors. The vectors are+-- transposed, so @fromVectorT@ is the same as @transpose+-- . fromVector@. It *is* verified that the vectors have the+-- right length.+fromVectorT :: (G.Vector v (w Double), G.Vector w Double)+ => v (w Double) -> Matrix+fromVectorT v =+ M { mRows = c+ , mCols = r+ , mData = unsafePerformIO $ do+ m <- MU.new (r*c)+ fillCol m r+ G.unsafeFreeze m}+ where+ r = G.length v+ c = G.length (G.head v)+ fillCol _ 0 = return ()+ fillCol m j = let j' = j-1+ in fillRow m (v G.! j') j' c >> fillCol m j'+ fillRow _ _ _ 0 = return ()+ fillRow m clm j' i = let i' = i-1+ x = clm G.! i'+ in MU.write m (i' * r + j') x >> fillRow m clm j' i'+++++-- | /O(rows)/ Rows of the matrix. Each element takes /O(1)/ time and+-- storage.+rows :: Matrix -> V.Vector (U.Vector Double)+rows m = G.map (\i -> U.unsafeSlice i (mCols m) (mData m)) $+ G.enumFromStepN 0 (mCols m) (mRows m)++-- | /O(1)/ @m !!! i@ is the @i@-th row of the matrix.+(!!!) :: Matrix -> Int -> U.Vector Double+m !!! i = U.slice (i * mCols m) (mCols m) (mData m)+++++++-- | /O(rows*cols)/ Columns of the matrix. Each element takes+-- /O(rows)/ time and storage.+cols :: Matrix -> V.Vector (U.Vector Double)+cols m = V.generate (mCols m) (m `col`)++-- | /O(rows)/ @m `col` i@ is the @i@-th column of the matrix.+col :: Matrix -> Int -> U.Vector Double+m `col` i = U.backpermute (mData m) $ U.enumFromStepN i (mCols m) (mRows m)+++++++umap :: (U.Vector Double -> U.Vector Double) -> Matrix -> Matrix+umap f m = m {mData = f (mData m)}++map :: (Double -> Double) -> Matrix -> Matrix+map f = umap (U.map f)++imap :: ((Int,Int) -> Double -> Double) -> Matrix -> Matrix+imap f m = umap (U.imap (f . indices m)) m++rowmap :: (U.Vector Double -> Double) -> Matrix -> U.Vector Double+rowmap f m = U.generate (mRows m) (f . s)+ where s i = U.unsafeSlice (i * mCols m) (mCols m) (mData m)++irowmap :: (Int -> U.Vector Double -> Double) -> Matrix -> U.Vector Double+irowmap f m = U.generate (mRows m) (\i -> f i $ s i)+ where s i = U.unsafeSlice (i * mCols m) (mCols m) (mData m)++uzipWith :: (U.Vector Double -> U.Vector Double -> U.Vector Double)+ -> Matrix -> Matrix -> Matrix+uzipWith f m n+ | mRows m /= mRows n = materror "uzipWith" "mRows"+ | mCols m /= mCols n = materror "uzipWith" "mCols"+ | otherwise = m {mData = f (mData m) (mData n)}++zipWith :: (Double -> Double -> Double) -> Matrix -> Matrix -> Matrix+zipWith f = uzipWith (U.zipWith f)++izipWith :: ((Int,Int) -> Double -> Double -> Double)+ -> Matrix -> Matrix -> Matrix+izipWith f m = uzipWith (U.izipWith (f . indices m)) m++-- | @rzipWith f m n@ is a matrix with the same number of rows as+-- @m@. The @i@-th row is obtained by applying @f@ to the @i@-th+-- rows of @m@ and @n@.+rzipWith :: (Int -> U.Vector Double -> U.Vector Double -> U.Vector Double)+ -> Matrix -> Matrix -> Matrix+rzipWith f m n+ | rm /= rn = materror "rzipWithN" $ "mRows " ++ s+ | cm /= cn = materror "rzipWithN" $ "mCols " ++ s+ | otherwise = fromVector $ V.izipWith f (rows m) (rows n)+ where rm = mRows m; cm = mCols m+ rn = mRows n; cn = mCols n+ s = show ((rm,cm),(rn,cn))+++indices :: Matrix -> Int -> (Int, Int)+indices m i = i `divMod` mCols m++++transpose :: Matrix -> Matrix+transpose m =+ let f i = let (r,c) = i `divMod` mRows m+ in m ! (c,r)+ in M { mRows = mCols m+ , mCols = mRows m+ , mData = U.generate (mRows m * mCols m) f}++{-# RULES+ "transpose/transpose" forall m. transpose (transpose m) = m;+ "transpose/fromVector" forall v. transpose (fromVector v) = fromVectorT v;+ "transpose/fromVectorT" forall v. transpose (fromVectorT v) = fromVector v;+ #-}++++materror :: String -> String -> a+materror f e = error $ "Math.Statistics.Dirichlet.Matrix." ++ f ++ ": " ++ e
+ src/Math/Statistics/Dirichlet/Mixture.hs view
@@ -0,0 +1,533 @@+---------------------------------------------------------------------------+-- | Module : Math.Statistics.Dirichlet.Mixture+-- Copyright : (c) 2009-2012 Felipe Lessa+-- License : BSD3+--+-- Maintainer : felipe.lessa@gmail.com+-- Stability : experimental+-- Portability : portable+--+--------------------------------------------------------------------------++module Math.Statistics.Dirichlet.Mixture+ ( -- * Data types+ DirichletMixture(..)+ , dmComponents+ , dmParameters+ , dmDensitiesL+ , (!!!)+ , empty+ , Component+ , fromList+ , toList+ , fromDD+ -- * Training data+ , TrainingData+ , prepareTraining+ -- * Functions+ , derive+ , cost+ , del_cost_w+ ) where++import qualified Data.Vector as V+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed as U++import Control.DeepSeq (NFData(..))+import Control.Monad.ST+import Data.Bits+import Data.Function (fix)+import Numeric.GSL.Special.Gamma (lngamma)+import Numeric.GSL.Special.Psi (psi)++import qualified Numeric.Optimization.Algorithms.HagerZhang05 as CG++import qualified Math.Statistics.Dirichlet.Density as D+import qualified Math.Statistics.Dirichlet.Matrix as M+import Math.Statistics.Dirichlet.Density (DirichletDensity(..))+import Math.Statistics.Dirichlet.Matrix (Matrix (..))+import Math.Statistics.Dirichlet.Options+import Math.Statistics.Dirichlet.Util++++-- | A Dirichlet mixture.+data DirichletMixture =+ DM { dmWeights :: !(U.Vector Double)+ -- ^ Weights of each density.+ , dmDensities :: !M.Matrix+ -- ^ Values of all parameters of all densities. This+ -- matrix has @length dmWeights@ rows.+ } deriving (Eq)++instance Show DirichletMixture where+ showsPrec prec dm =+ showParen (prec > 10) $+ showString "fromList " .+ showsPrec 11 (toList dm)++instance Read DirichletMixture where+ readsPrec p ('(':xs) = let (ys,')':zs) = break (== ')') xs+ in map (\(x,s) -> (x,s++zs)) $+ readsPrec p ys+ readsPrec p xs = let [("fromList",list)] = lex xs+ in map (\(x,s) -> (fromList x,s)) $+ readsPrec p list++instance NFData DirichletMixture where+ rnf DM {} = ()+++-- | Number of components in a dirichlet mixture.+dmComponents :: DirichletMixture -> Int+dmComponents = U.length . dmWeights++-- | Number of parameters each component has.+dmParameters :: DirichletMixture -> Int+dmParameters = mCols . dmDensities++-- | Separated list of densities.+dmDensitiesL :: DirichletMixture -> [DirichletDensity]+dmDensitiesL (DM _ as) = map DD $ V.toList $ M.rows as++-- | @dm !!! i@ is the @i@-th density. No bounding checks are+-- made.+(!!!) :: DirichletMixture -> Int -> U.Vector Double+(DM _ as) !!! i = as M.!!! i+{-# INLINE (!!!) #-}++++++dmap :: (U.Vector Double -> Double) -> DirichletMixture -> U.Vector Double+dmap f = M.rowmap f . dmDensities++++-- | @empty q n x@ is an \"empty\" Dirichlet mixture with @q@+-- components and @n@ parameters. Each component has size @n@,+-- weight inversely proportional to its index and all alphas set+-- to @x@.+empty :: Int -> Int -> Double -> DirichletMixture+empty q n x = let (DD d) = D.empty n x+ f i = fromIntegral (q-i) / sum_+ sum_ = fromIntegral (q*(q+1)`div`2)+ in DM {dmWeights = U.generate q f+ ,dmDensities = M.replicateRows q d}+{-# INLINE empty #-}+++-- | A list representation of a component of a Dirichlet mixture.+-- Used by 'fromList' and 'toList' only.+type Component = (Double, [Double])++-- | @fromList xs@ constructs a Dirichlet mixture from a+-- non-empty list of components. Each component has a weight and+-- a list of alpha values. The weights sum to 1, all lists must+-- have the same number of values and every number must be+-- non-negative. None of these preconditions are verified.+fromList :: [Component] -> DirichletMixture+fromList components =+ let -- Vectors+ qs = U.fromList $ map fst components+ as = M q n $ U.fromList $ concatMap snd components++ -- Properties of the mixture+ q = length components+ n = length (snd $ head components)+ in DM qs as++-- | @toList dm@ is the inverse of @fromList@, constructs a list+-- of components from a Dirichlet mixture. There are no error+-- conditions and @toList . fromList == id@.+toList :: DirichletMixture -> [Component]+toList dm =+ let qs' = U.toList $ dmWeights dm+ as' = map (U.toList . unDD) (dmDensitiesL dm)+ in zip qs' as'++-- | Constructs a Dirichlet mixture of one component from a+-- Dirichlet density.+fromDD :: DirichletDensity -> DirichletMixture+fromDD (DD d) = DM (U.singleton 1) (M.replicateRows 1 d)++++++-- | Prepares training vectors to be used as training data.+-- Anything that depends only on the training vectors is+-- precalculated here.+--+-- We also try to find columns where all training vectors are+-- zero. Those columns are removed from the derivation process+-- and every component will have zero value on that column. Note+-- that at least one column should have non-zero training+-- vectors.+prepareTraining :: TrainingVectors -> TrainingData+prepareTraining ns_0 =+ let zeroes = zeroedCols ns_0+ ns = removeZeroes ns_0 zeroes+ ns_sums = G.unstream $ G.stream $ V.map U.sum ns+ tns = M.fromVectorT ns+ in TD {..}++-- | Pre-processed training vectors (see 'prepareTraining').+data TrainingData = TD { ns :: !TrainingVectors+ , ns_sums :: !(U.Vector Double)+ , tns :: !Matrix+ , zeroes :: ![Int]}+ deriving (Eq, Show)++-- | Return the list of columns that are zeroed, counting from zero.+zeroedCols :: TrainingVectors -> [Int]+zeroedCols =+ -- We set the i-th bit whenever the i-th column was zeroed.+ let fold (acc, mask) 0 = (acc .|. mask, shiftL mask 1)+ fold (acc, mask) _ = (acc :: Integer, shiftL mask 1)+ unBits !_ 0 = []+ unBits !i x = (if testBit x 0 then (i:) else id)+ (unBits (i+1) (shiftR x 1))+ in unBits 0 . V.foldl1' (.&.) . V.map (fst . U.foldl' fold (0,1))++-- | Remove zeroed columns from training vectors.+removeZeroes :: TrainingVectors -> [Int] -> TrainingVectors+removeZeroes ns [] = ns+removeZeroes ns zs =+ let cols_orig = U.length (V.head ns)+ cols_new = U.filter (`notElem` zs) $ U.enumFromN 0 cols_orig+ in V.map (flip U.backpermute cols_new) ns++-- | Remove zeroed columns from a Dirichlet mixture matrix of+-- densities.+removeZeroesM :: [Int] -> Matrix -> Matrix+removeZeroesM [] as = as+removeZeroesM zs as =+ let size = M.mCols as * M.mRows as+ cols_orig = M.mCols as+ cols_new = U.filter ((`notElem` zs) . (`rem` cols_orig)) $+ U.enumFromN 0 size+ in M {mCols = M.mCols as - length zs+ ,mRows = M.mRows as+ ,mData = U.backpermute (M.mData as) cols_new}++-- | Add zeroed columns back to a Dirichlet mixture matrix of+-- densities.+addZeroesM :: [Int] -> Matrix -> Matrix+addZeroesM [] = id+addZeroesM zs' = M.fromVector .+ V.map (U.fromList . add 0 zs' . U.toList) .+ M.rows+ where+ add !_ [] xs = xs+ add _ zs [] = map (const zero) zs+ add i (z:zs) (x:xs) | i == z = zero : add (i+1) zs (x:xs)+ | otherwise = x : add (i+1) (z:zs) xs+ zero = 0.00001++++++-- | /Prob(a_j | n, theta)/ Defined in equation (16), "the+-- posterior probability of the j-th component of the mixture+-- given the vector of counts n". We return the probabilities+-- for all /j/ in each vector.+--+-- The order of the result is inversed for performance. In the+-- outer boxed vector there are /j/ elements. The /i/-th inner+-- unboxed vector contains that probability for each of the+-- training vectors.+--+-- Calculated as per equation (39) using 'logBeta'. If we take+-- the numerator of the right hand side of equation (39) as /Y_j/+-- and the left hand side as /P_j/, then /P_j/ is proportional to+-- /Y_j/ normalized to sum to 1. We may have problems if /P_j/+-- is too large or too small. Using the suggestion from the+-- paper, we may multiply all /P_j/ by a constant before+-- normalizing everything. We calculate /P_j/ using a logarithm,+-- so that means we may freely add or subtract a constant from+-- the logarithm before appling the exponential function. This+-- is really essencial.+prob_a_n_theta :: TrainingVectors -> DirichletMixture -> Matrix+prob_a_n_theta ns dm@(DM qs _) =+ let -- Precalculate logBeta of all components+ !logBetaAlphas = dmap logBeta dm++ -- Calculate the factors for one of the training vectors.+ calc n i lb_a = let !a = dm !!! i+ in logBeta (U.zipWith (+) n a) - lb_a+ factors n = let fs = U.imap (calc n) logBetaAlphas+ !c = U.maximum fs -- see the note above+ fs' = U.zipWith (\q f -> q * exp (f - c)) qs fs+ !total = U.sum fs'+ in U.map (/ total) fs'+ in M.fromVectorT $ V.map factors ns+++-- | Customized version of @prob_a_n_theta@ used when the weights+-- are being estimated. Precomputes everything that doesn't+-- depend on the weight.+prob_a_n_theta_weights :: TrainingVectors -> Matrix+ -> (U.Vector Double -> Matrix)+prob_a_n_theta_weights ns as =+ let -- Precalculate logBeta of all components+ !logBetaAlphas = M.rowmap logBeta as++ -- Precalculate the factors for one of the training vectors.+ precalc n i lb_a = let !a = as M.!!! i+ in logBeta (U.zipWith (+) n a) - lb_a+ norm fs = let !c = U.maximum fs+ in U.map (exp . subtract c) fs+ !prefactors = V.map (norm . flip U.imap logBetaAlphas . precalc) ns++ in \qs ->+ let -- Calculate the final factors.+ calc pfs = let fs = U.zipWith (*) pfs qs+ total = U.sum fs+ in U.map (/ total) fs+ in M.fromVectorT $ V.map calc prefactors+++++++++++++-- | Cost function for deriving a Dirichlet mixture (equation+-- 18). This function is minimized by 'derive'. Calculated+-- using (17) and (54).+cost :: TrainingData -> DirichletMixture -> Double+cost td dm =+ let as_sums = dmap U.sum dm+ in cost_worker td dm as_sums+++-- | Worker of 'cost' function that avoids repeating some+-- computations that are done when reestimating alphas.+cost_worker :: TrainingData -> DirichletMixture+ -> U.Vector Double -> Double+cost_worker TD {ns, ns_sums} dm@(DM !qs _) !as_sums =+ let -- From the equation (54).+ prob_n_a !n !n_sum !a !a_sum !lngamma_a_sum =+ let !s = lngamma (n_sum+1) + lngamma_a_sum - lngamma (n_sum+a_sum)+ f n_i a_i = lngamma (n_i + a_i) - lngamma (n_i + 1) - lngamma a_i+ in exp $ s + U.sum (U.zipWith f n a)++ -- From equation (17).+ prob_n_theta i n =+ let !n_sum = ns_sums U.! i+ in U.sum $ U.zipWith (*) qs $+ U.izipWith (prob_n_a n n_sum . (dm !!!))+ as_sums lngamma_as_sums+ !lngamma_as_sums = U.map lngamma as_sums+ in negate $ V.sum $ V.imap ((log .) . prob_n_theta) ns++-- | Version of 'cost' function that avoids repeating a lot of+-- computations that are done when reestimating weights.+cost_weight :: TrainingData -> Matrix+ -> U.Vector Double -> (U.Vector Double -> Double)+cost_weight TD {ns, ns_sums} !as !as_sums =+ let -- From the equation (54).+ prob_n_a !n !n_sum !a !a_sum !lngamma_a_sum =+ let !s = lngamma (n_sum+1) + lngamma_a_sum - lngamma (n_sum+a_sum)+ f n_i a_i = lngamma (n_i + a_i) - lngamma (n_i + 1) - lngamma a_i+ in exp $ s + U.sum (U.zipWith f n a)++ -- From equation (17).+ prepare_prob_n_theta i n =+ let !n_sum = ns_sums U.! i+ in {- U.sum $ U.zipWith (*) qs $ -}+ U.izipWith (prob_n_a n n_sum . (as M.!!!))+ as_sums lngamma_as_sums+ !lngamma_as_sums = U.map lngamma as_sums+ !prepared = V.imap prepare_prob_n_theta ns++ -- Final worker function.+ final qs = log . U.sum . U.zipWith (*) qs+ in \(!qs) -> negate $ V.sum $ V.map (final qs) prepared++++++++-- | Derivative of the cost function with respect @w_{i,j}@,+-- defined by Equation (22). The result is given in the same+-- size and order as the 'dmDensitites' vector.+del_cost_w :: TrainingData -> DirichletMixture -> Matrix+del_cost_w td dm =+ let as_sums = dmap U.sum dm+ in del_cost_w_worker td dm as_sums+++-- | Worker function of 'del_cost_w'.+del_cost_w_worker :: TrainingData -> DirichletMixture+ -> U.Vector Double -> Matrix+del_cost_w_worker TD {ns, ns_sums, tns} dm !as_sums =+ let -- Calculate Prob(a | n, theta)+ !probs_a_n = prob_a_n_theta ns dm++ -- Calculate all S_j's.+ !sjs = M.rowmap U.sum probs_a_n++ -- @calc j _ i _ _@ calculates the derivative of the+ -- cost function with respect to @w_{i,j}@. The other+ -- arguments come from vector that we @zipWith@ below.+ calc j probs =+ -- Everything that doesn't depend on i, just on j.+ let !a_sum = as_sums U.! j+ !psi_a_sum = psi a_sum+ !sum_prob_psi = U.sum $ U.zipWith (*) probs $+ U.map (psi . (+) a_sum) ns_sums+ -----+ in \i a_i ->+ let !s1 = (sjs U.! j) * (psi_a_sum - psi a_i)+ !s2 = U.sum $ U.zipWith (\p_i n_i -> p_i * psi (n_i + a_i)) probs (tns M.!!! i)+ in - a_i * (s1 + s2 - sum_prob_psi)++ in M.fromVector $ V.imap (\j p_j -> let !f = calc j p_j+ in U.imap f (dm !!! j))+ (M.rows probs_a_n)++++++-- | Derive a Dirichlet mixture using a maximum likelihood method+-- as described by Karplus et al (equation 25) using CG_DESCENT+-- method by Hager and Zhang (see+-- "Numeric.Optimization.Algorithms.HagerZhang05"). All training+-- vectors should have the same length, however this is not+-- verified.+derive :: DirichletMixture -> Predicate -> StepSize+ -> TrainingData -> Result DirichletMixture+derive (DM initial_qs initial_as') (Pred {..}) _ td@(TD {ns,zeroes})+ | V.length ns == 0 = err "empty training data"+ | U.length initial_qs < 1 = err "empty initial weights vector"+ | M.size initial_as < (1,1) = err "empty initial alphas vector"+ | maxIter < 1 = err "non-positive maxIter"+ | minDelta < 0 = err "negative minDelta"+ | jumpDelta < 0 = err "negative jumpDelta"+ | jumpDelta < minDelta = err "minDelta greater than jumpDelta"+ | otherwise = runST train+ where+ err = error . ("Dirichlet.derive: " ++)+ singleDensity = U.length initial_qs == 1++ -- Remove zeroes from initial_as'.+ initial_as = removeZeroesM zeroes initial_as'++ -- Reciprocal of the number of training sequences.+ !recip_m = recip $ fromIntegral $ V.length ns++ -- Calculate the sums of the alphas.+ calc_as_sums = M.rowmap U.sum++ -- Parameters used by CG_DESCENT.+ verbose = False+ parameters = CG.defaultParameters+ { CG.printFinal = verbose+ , CG.printParams = verbose+ , CG.verbose = if verbose then CG.VeryVerbose else CG.Quiet+ , CG.maxItersFac = max 1 $ fromIntegral maxIter / 20+ , CG.estimateError = CG.RelativeEpsilon (1e-6 * s)+ }+ where (w,h) = M.size initial_as+ s = fromIntegral (w * h * V.length ns)++ -- Transform a U.Vector from/to a M.Matrix in the case that+ -- the matrix has the same shape as initial_as (i.e. all+ -- as's and ws's).+ fromMatrix = M.mData+ toMatrix v = initial_as {M.mData = v}++ -- Create specialized functions that are optimized by+ -- CG_DESCENT. They depend only on @qs@, the weights.+ createFunctions !qs =+ let calc f = \ws -> let !as = M.map exp (toMatrix ws)+ !as_sums = calc_as_sums as+ dm = DM qs as+ in f dm as_sums+ grad_worker = ((fromMatrix .) .) . del_cost_w_worker+ func = CG.VFunction $ calc $ cost_worker td+ grad = CG.VGradient $ calc $ grad_worker td+ comb = CG.VCombined $ calc $ \dm as_sums ->+ (cost_worker td dm as_sums+ ,grad_worker td dm as_sums)+ in (func, grad, comb)++ -- Start training in the zero-th iteration and with+ -- infinite inital cost.+ train = trainAlphas 0 infinity initial_qs $ M.map log initial_as++ trainAlphas !iter !oldCost !qs !ws = {-# SCC "trainAlphas" #-} do+ -- Optimize using CG_DESCENT+ let (func, grad, comb) = createFunctions qs+ opt = CG.optimize parameters minDelta (fromMatrix ws)+ func grad (Just comb)++ (!pre_ws', result, stats) <- unsafeIOToST opt+ let !ws' = toMatrix (G.unstream $ G.stream pre_ws')++ -- Recalculate everything.+ let !as' = M.map exp ws'+ as_sums' = calc_as_sums as'+ !iter' = iter + fromIntegral (CG.totalIters stats)+ !cost' = CG.finalValue stats+ !delta = abs (cost' - oldCost)+ dm = DM qs $ addZeroesM zeroes as'++ -- Verify convergence. Even with MaxIter we only stop+ -- iterating if the delta was calculated. Otherwise we+ -- wouldn't be able to tell the caller why the delta was+ -- still big when we reached the limit.+ case (decide result+ ,delta <= minDelta+ ,iter' >= maxIter+ ,singleDensity) of+ (Stop r,_,_,_) -> return $ Result r iter' delta cost' dm+ (_,True,_,_) -> return $ Result Delta iter' delta cost' dm+ (_,_,True,_) -> return $ Result MaxIter iter' delta cost' dm+ (_,_,_,True) -> return $ Result Delta iter' delta cost' dm+ (GoOn,_,_,_) -> trainWeights iter' cost' qs ws' as' as_sums'++ trainWeights !oldIter !veryOldCost !oldQs !ws !as !as_sums =+ {-# SCC "trainWeights" #-}+ -- Prepare invariant parts.+ let !probs_a_n_mk = prob_a_n_theta_weights ns as+ !cost_mk = cost_weight td as as_sums+ in ($ oldQs) . ($ veryOldCost) . ($ maxWeightIter) . fix $+ \again !itersLeft !oldCost !qs ->+ -- Reestimate weight's.+ let !probs_a_n = probs_a_n_mk qs+ qs' = M.rowmap ((*) recip_m . U.sum) probs_a_n++ -- Recalculate constants.+ !cost' = cost_mk qs'+ !delta = abs (cost' - oldCost)++ -- Verify convergence. We never stop the process here.+ in case (delta <= jumpDelta, itersLeft <= 0) of+ (False,False) -> again (itersLeft-1) cost' qs'+ _ -> trainAlphas oldIter cost' qs' ws+++-- | Decide what we should do depending on the result of the+-- CG_DESCENT routine.+decide :: CG.Result -> Decision+decide CG.ToleranceStatisfied = GoOn+decide CG.FunctionChange = GoOn+decide CG.MaxTotalIter = GoOn+decide CG.MaxSecantIter = GoOn+decide other = Stop (CG other)++data Decision = GoOn | Stop Reason
+ src/Math/Statistics/Dirichlet/Options.hs view
@@ -0,0 +1,91 @@+---------------------------------------------------------------------------+-- | Module : Math.Statistics.Dirichlet.Options+-- Copyright : (c) 2009-2012 Felipe Lessa+-- License : BSD3+--+-- Maintainer : felipe.lessa@gmail.com+-- Stability : experimental+-- Portability : portable+--+--------------------------------------------------------------------------++module Math.Statistics.Dirichlet.Options+ ( TrainingVector+ , TrainingVectors+ , StepSize(..)+ , Delta+ , Predicate(..)+ , Reason(..)+ , Result(..)+ ) where++import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as U+import qualified Numeric.Optimization.Algorithms.HagerZhang05 as CG+import Control.DeepSeq (NFData(..))++-- | A vector used for deriving the parameters of a Dirichlet+-- density or mixture.+type TrainingVector = U.Vector Double++-- | A vector of training vectors. This is the only vector that+-- is not unboxed (for obvious reasons).+type TrainingVectors = V.Vector TrainingVector++-- | Usually denoted by lowercase greek letter eta (η), size of+-- each step in the gradient. Should be greater than zero and+-- much less than one.+newtype StepSize = Step Double++-- | Maximum difference between costs to consider that the+-- process converged.+type Delta = Double++-- | Predicate specifying when the training should be over.+data Predicate = Pred+ { maxIter :: !Int -- ^ Maximum number of iterations.+ , minDelta :: !Delta -- ^ Minimum delta to continue iterating.+ -- This is invariant of @deltaSteps@, which+ -- means that if @deltaSteps@ is @2@ then+ -- minDelta will be considered twice bigger+ -- to account for the different @deltaSteps@.+ , deltaSteps :: !Int -- ^ How many estimation steps should be done+ -- before recalculating the delta. If+ -- @deltaSteps@ is @1@ then it will be+ -- recalculated on every step.+ , maxWeightIter :: !Int -- ^ Maximum number of iterations on+ -- each weight step.+ , jumpDelta :: !Delta -- ^ Used only when calculating mixtures.+ -- When the delta drops below this cutoff+ -- the computation changes from estimating+ -- the alphas to estimatating the weights+ -- and vice-versa. Should be greater than+ -- @minDelta@.+ }+ deriving (Eq, Read, Show)++-- | Reason why the derivation was over.+data Reason = Delta -- ^ The difference between+ -- applications of the cost function+ -- dropped below the minimum delta.+ -- In other words, it coverged.+ | MaxIter -- ^ The maximum number of iterations+ -- was reached while the delta was+ -- still greater than the minimum delta.+ | CG CG.Result -- ^ CG_DESCENT returned this result,+ -- which brought the derivation+ -- process to a halt.+ deriving (Eq, Read, Show)++-- | Result of a deriviation.+data Result a =+ Result { reason :: !Reason -- ^ Reason why the derivation was over.+ , iters :: !Int -- ^ Number of iterations spent.+ , lastDelta :: !Delta -- ^ Last difference between costs.+ , lastCost :: !Double -- ^ Last cost (i.e. the cost of the result).+ , result :: !a -- ^ Result obtained.+ }+ deriving (Eq, Read, Show)++instance NFData a => NFData (Result a) where+ rnf = rnf . result
+ src/Math/Statistics/Dirichlet/Util.hs view
@@ -0,0 +1,32 @@+---------------------------------------------------------------------------+-- | Module : Math.Statistics.Dirichlet.Util+-- Copyright : (c) 2009-2012 Felipe Lessa+-- License : BSD3+--+-- Maintainer : felipe.lessa@gmail.com+-- Stability : experimental+-- Portability : portable+--+--------------------------------------------------------------------------++module Math.Statistics.Dirichlet.Util+ ( infinity+ , logBeta+ )+ where++import qualified Data.Vector.Unboxed as U+import Numeric.GSL.Special.Gamma (lngamma, lnbeta)++++-- | Logarithm of the beta function applied to a vector.+logBeta :: U.Vector Double -> Double+logBeta xs | U.length xs == 2 = lnbeta (U.head xs) (U.last xs)+ | otherwise = U.sum (U.map lngamma xs) - lngamma (U.sum xs)++-- | Infinity, currently defined as @1e100@. Used mainly as the+-- initial cost.+infinity :: Double+infinity = 1e100+{-# INLINE infinity #-}
+ statistics-dirichlet.cabal view
@@ -0,0 +1,49 @@+Cabal-Version: >= 1.6+Build-Type: Simple+Tested-With: GHC+Category: Math+Name: statistics-dirichlet+Version: 0.6+Stability: experimental+License: BSD3+License-File: LICENSE+Copyright: (c) 2009-2012 Felipe A. Lessa+Author: Felipe Almeida Lessa+Maintainer: felipe.lessa@gmail.com+Synopsis: Functions for working with Dirichlet densities and mixtures on vectors.++Description:+ Functions for working with Dirichlet densities and mixtures+ on vectors. The focus of this package is on deriving these+ distributions from observed data.+ .+ This package should be treated as experimental code, it has+ not been battle-tested as much as it would be nice to be.+ .+ Note that although this package is BSD3-licensed, it uses the+ @nonlinear-optimization@ package which is GPLed. It should+ be straightforward to use another library in its stead,+ though.++Source-repository head+ type: darcs+ location: http://patch-tag.com/r/felipe/statistics-dirichlet++Library+ Build-Depends:+ base == 4.*+ , deepseq == 1.3.*+ , vector == 0.9.*+ , nonlinear-optimization == 0.3.*+ , hmatrix-special == 0.1.*++ Ghc-Options: -Wall+ Extensions: BangPatterns, NamedFieldPuns, RecordWildCards, FlexibleContexts+ Exposed-Modules:+ Math.Statistics.Dirichlet,+ Math.Statistics.Dirichlet.Density,+ Math.Statistics.Dirichlet.Matrix,+ Math.Statistics.Dirichlet.Mixture,+ Math.Statistics.Dirichlet.Options,+ Math.Statistics.Dirichlet.Util+ hs-Source-Dirs: src/