diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2013, Jonathan Fischoff + +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 Jonathan Fischoff 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.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/maxent.cabal b/maxent.cabal new file mode 100644 --- /dev/null +++ b/maxent.cabal @@ -0,0 +1,78 @@ +-- Initial maxent.cabal generated by cabal init. For further +-- documentation, see http://haskell.org/cabal/users-guide/ + +-- The name of the package. +name: maxent + +-- The package version. See the Haskell package versioning policy (PVP) +-- for standards guiding when and how versions should be incremented. +-- http://www.haskell.org/haskellwiki/Package_versioning_policy +-- PVP summary: +-+------- breaking API changes +-- | | +----- non-breaking API additions +-- | | | +--- code changes with no API change +version: 0.1.0.0 + +-- A short (one-line) description of the package. +synopsis: Compute Maximum Entropy Distrubtions + +-- A longer description of the package. +description: Use this package to compute maximum entropy distributions given a list of values and + list of constraints. + . + Here is a the example from Probability the Logic of Science + . + > maxent ([1,2,3], [average 1.5]) + . + Right [0.61, 0.26, 0.11] + . + The classic dice example + . + > maxent ([1,2,3,4,5,6], [average 4.5]) + . + Right [.05, .07, 0.11, 0.16, 0.23, 0.34] + + I will document this more ... soonish + +-- URL for the project homepage or repository. +homepage: https://github.com/jfischoff/maxent + +-- The license under which the package is released. +license: BSD3 + +-- The file containing the license text. +license-file: LICENSE + +-- The package author(s). +author: Jonathan Fischoff + +-- An email address to which users can send suggestions, bug reports, and +-- patches. +maintainer: jonathangfischoff@gmail.com + +-- A copyright notice. +-- copyright: + +category: Math + +build-type: Simple + +-- Constraint on the version of Cabal needed to build this package. +cabal-version: >=1.8 + + +library + -- Modules exported by the library. + exposed-modules: MaxEnt + + -- Modules included in this library but not exported. + other-modules: MaxEnt.Internal + + -- Other library packages from which modules are imported. + build-depends: base ==4.6.*, + nonlinear-optimization ==0.3.*, + vector ==0.9.*, + ad ==3.2.* + + -- Directories containing source files. + hs-source-dirs: src + diff --git a/src/MaxEnt.hs b/src/MaxEnt.hs new file mode 100644 --- /dev/null +++ b/src/MaxEnt.hs @@ -0,0 +1,27 @@ +-- | +-- Use this package to compute maximum entropy distributions given a list of values and +-- list of constraints. +-- +-- Here is a the example from Probability the Logic of Science +-- +-- > maxent ([1,2,3], [average 1.5]) +-- +-- Right [0.61, 0.26, 0.11] +-- +-- The classic dice example +-- +-- > maxent ([1,2,3,4,5,6], [average 4.5]) +-- +-- Right [.05, .07, 0.11, 0.16, 0.23, 0.34] +module MaxEnt ( + Constraint, + constraint, + average, + variance, + maxent +) where +import MaxEnt.Internal (Constraint, + constraint, + average, + variance, + maxent) diff --git a/src/MaxEnt/Internal.hs b/src/MaxEnt/Internal.hs new file mode 100644 --- /dev/null +++ b/src/MaxEnt/Internal.hs @@ -0,0 +1,87 @@ +{-# LANGUAGE TupleSections, Rank2Types #-} +module MaxEnt.Internal 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) + +sumWith :: Num c => (a -> b -> c) -> [a] -> [b] -> c +sumWith f xs = sum . zipWith f xs + +pOfK :: Floating a => [a] -> [a -> a] -> [a] -> Int -> a +pOfK values fs ls k = exp (negate . sumWith (\l f -> l * f (values !! k)) ls $ fs) / + partitionFunc values fs ls + +probs :: Floating b => [b] -> [b -> b] -> [b] -> [b] +probs values fs ls = map (pOfK values fs ls) [0..length values - 1] + +partitionFunc :: Floating a => [a] -> [a -> a] -> [a] -> a +partitionFunc values fs ls = sum $ [ exp ((-l) * f x) | x <- values, (f, l) <- zip fs ls] + +objectiveFunc :: Floating a => [a] -> [a -> a] -> [a] -> [a] -> a +objectiveFunc values fs moments ls = log (partitionFunc values fs ls) + sumWith (*) ls moments + +toFunction :: (forall a. Floating a => [a] -> a) -> Function Simple +toFunction f = VFunction (f . U.toList) + +toGradient :: (forall a. Floating a => [a] -> a) -> Gradient Simple +toGradient f = VGradient (U.fromList . grad f . U.toList) + +toDoubleF :: (forall a. Floating a => [a] -> a) -> [Double] -> Double +toDoubleF f x = f x + +-- | Constraint type. Think of this as f and c in sum pi (f x) = c +type Constraint a = (a -> a, a) + +-- make a constraint from function and constant +constraint :: Floating a => (a -> a) -> a -> Constraint a +constraint = (,) + +-- The average constraint +average :: Floating a => a -> Constraint a +average m = constraint id m + +-- The variance constraint +variance :: Floating a => a -> Constraint a +variance sigma = constraint (^(2 :: Int)) sigma + +-- | The main entry point for computing discrete maximum entropy distributions. +-- +maxent :: (forall a. Floating a => ([a], [Constraint a])) -- ^ A pair of values that the distributions is over and the constraints + -> Either (Result, Statistics) [Double] -- ^ Either the a discription of what wrong or the probability distribution +maxent params = result where + obj :: Floating a => [a] -> a + obj = uncurry (objectiveFunc values) fsmoments + + values :: Floating a => [a] + values = fst params + + constraints :: Floating a => [(a -> a, a)] + constraints = snd params + + fsmoments :: Floating a => ([a -> a], [a]) + fsmoments = unzip constraints + + fs :: [Double -> Double] + fs = fst fsmoments + + -- hmm maybe there is a better way to get rid of the defaulting + guess = U.fromList $ replicate + (length (constraints :: [(Double -> Double, Double)])) (1.0 :: Double) + + result = case unsafePerformIO (optimize defaultParameters 0.00001 guess + (toFunction obj) + (toGradient obj) + Nothing) of + (vs, ToleranceStatisfied, _) -> Right $ probs values fs (S.toList vs) + (_, x, y) -> Left (x, y) + +--test = maxent ([1.0,2.0,3.0], [average 1.5]) + + + + + + +