diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,3 @@ +#!/usr/bin/runhaskell +> import Distribution.Simple +> main = defaultMain diff --git a/src/Math/Statistics/Dirichlet.hs b/src/Math/Statistics/Dirichlet.hs new file mode 100644 --- /dev/null +++ b/src/Math/Statistics/Dirichlet.hs @@ -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 diff --git a/src/Math/Statistics/Dirichlet/Density.hs b/src/Math/Statistics/Dirichlet/Density.hs new file mode 100644 --- /dev/null +++ b/src/Math/Statistics/Dirichlet/Density.hs @@ -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 #-} diff --git a/src/Math/Statistics/Dirichlet/Matrix.hs b/src/Math/Statistics/Dirichlet/Matrix.hs new file mode 100644 --- /dev/null +++ b/src/Math/Statistics/Dirichlet/Matrix.hs @@ -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 diff --git a/src/Math/Statistics/Dirichlet/Mixture.hs b/src/Math/Statistics/Dirichlet/Mixture.hs new file mode 100644 --- /dev/null +++ b/src/Math/Statistics/Dirichlet/Mixture.hs @@ -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 diff --git a/src/Math/Statistics/Dirichlet/Options.hs b/src/Math/Statistics/Dirichlet/Options.hs new file mode 100644 --- /dev/null +++ b/src/Math/Statistics/Dirichlet/Options.hs @@ -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 diff --git a/src/Math/Statistics/Dirichlet/Util.hs b/src/Math/Statistics/Dirichlet/Util.hs new file mode 100644 --- /dev/null +++ b/src/Math/Statistics/Dirichlet/Util.hs @@ -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 #-} diff --git a/statistics-dirichlet.cabal b/statistics-dirichlet.cabal new file mode 100644 --- /dev/null +++ b/statistics-dirichlet.cabal @@ -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/