maxent 0.6.0.0 → 0.6.0.1

raw patch · 3 files changed

+132/−2 lines, 3 files

Files

• maxent.cabal +4/−2
• src/Numeric/MaxEnt/General.hs +35/−0
• src/Numeric/MaxEnt/Moment.hs +93/−0

@@ -10,7 +10,7 @@ -- PVP summary:      +-+------- breaking API changes --                   | | +----- non-breaking API additions --                   | | | +--- code changes with no API change-version:             0.6.0.0+version:             0.6.0.1  -- A short (one-line) description of the package. synopsis:            Compute Maximum Entropy Distributions@@ -77,7 +77,9 @@   -- Modules included in this library but not exported.   other-modules:     Numeric.MaxEnt.Internal,                       Numeric.MaxEnt.Linear, -                     Numeric.MaxEnt.ConjugateGradient +                     Numeric.MaxEnt.ConjugateGradient,+                     Numeric.MaxEnt.Moment,+                     Numeric.MaxEnt.General        -- Other library packages from which modules are imported.   build-depends: base ==4.6.*,

@@ -0,0 +1,35 @@+{-# LANGUAGE TupleSections, Rank2Types #-}+module Numeric.MaxEnt.General (+    Constraint,+    general+ ) where+import Numeric.Optimization.Algorithms.HagerZhang05+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Storable as S+import Numeric.AD+import GHC.IO                   (unsafePerformIO)+import Data.Traversable+import Numeric.AD.Types+import Numeric.AD.Internal.Tower+import Numeric.AD.Internal.Classes+import Data.List (transpose)+import Control.Applicative+import Numeric.AD.Lagrangian++entropy :: Floating a => [a] -> a+entropy xs = negate . sum . map (\x -> x * log x) $ xs++-- | A more general solver. This directly solves the lagrangian of the constraints and the+--  the additional constraint that the probabilities must sum to one.+general :: Double +        -- ^ Tolerance for the numerical solver+        -> Int+        -- ^ the count of probabilities+        -> [Constraint Double]+        -- ^  constraints+        -> Either (Result, Statistics) (S.Vector Double) +        -- ^ Either the a discription of what wrong or the probability distribution+general tolerance count constraints = +    fst <$> maximize tolerance entropy ((sum <=> 1.0) : constraints) count+ +   

@@ -0,0 +1,93 @@+{-# LANGUAGE TupleSections, Rank2Types, NoMonomorphismRestriction #-}+module Numeric.MaxEnt.Moment (+        ExpectationConstraint,+        (.=.),+        ExpectationFunction,+        average,+        variance,+        maxent,+        UU(..)+    ) where+import Numeric.Optimization.Algorithms.HagerZhang05+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Storable as S+import Numeric.AD+import GHC.IO                   (unsafePerformIO)+import Data.Traversable+import Numeric.AD.Types+import Numeric.AD.Internal.Classes+import Data.List (transpose)+import Control.Applicative+import Numeric.MaxEnt.ConjugateGradient+import Data.List (foldl')+--import Data.Vector++-- | Constraint type. A function and the constant it equals.+-- +--   Think of it as the pair @(f, c)@ in the constraint +--+-- @+--     &#931; p&#8336; f(x&#8336;) = c+-- @+--+--  such that we are summing over all values .+--+--  For example, for a variance constraint the @f@ would be @(\\x -> x*x)@ and @c@ would be the variance.+type ExpectationConstraint a = (UU a, a)++--+infixr 1 .=.+(.=.) :: (forall s. Mode s => AD s a -> AD s a) -> a -> ExpectationConstraint a+f .=. c = (UU f, c)++-- | A function that takes an index and value and returns a value.+--   See 'average' and 'variance' for examples.+type ExpectationFunction a = (a -> a)++newtype UU a = UU {unUU :: forall s. Mode s => ExpectationFunction (AD s a) }++-- The average constraint+average :: Num a => a -> ExpectationConstraint a+average m = id .=. m++-- The variance constraint+variance :: Num a => a -> ExpectationConstraint a+variance sigma = (^(2 :: Int)) .=. sigma++pOfK :: [Double] -> [ExpectationFunction Double] -> S.Vector Double -> Int -> Double+pOfK values fs ls k = +    exp (negate . sum . zipWith (\l f -> l * f (values !! k)) lsList $ fs) / +        (partitionFunc values fs lsList) where+            lsList = S.toList ls++probs values fs ls = S.map (pOfK values fs ls) . S.enumFromN 0 $ length values ++partitionFunc values fs ls = sum $ [ exp ((-l) * f x) | +                                x <- values, +                                (f, l) <- zip fs ls]++objectiveFunc fs moments values ls = +    log (partitionFunc values fs ls) + (sum $ zipWith (\x y -> x * y) ls moments)++-- My thoughts are that I should maybe split this up+-- | Discrete maximum entropy solver where the constraints are all moment constraints. +maxent :: Double +       -- ^ Tolerance for the numerical solver+       -> [Double]+       -- ^ values that the distributions is over+       -> [ExpectationConstraint Double]+       -- ^ The constraints+       -> Either (Result, Statistics) (S.Vector Double) +       -- ^ Either the a discription of what wrong or the probability distribution +maxent tolerance values constraints = result where+    obj = objectiveFunc (map unUU fs') (map auto moments) (map auto values)+    +    count = length fs+        +    (fs', moments) = unzip constraints +    +    fs = map (\x -> lowerUU $ unUU x) fs'+    +    guess = U.fromList $ replicate count (1.0 / fromIntegral count :: Double) +    +    result =  probs values fs <$> minimize tolerance count obj