diff --git a/maxent.cabal b/maxent.cabal --- a/maxent.cabal +++ b/maxent.cabal @@ -10,7 +10,7 @@ -- PVP summary: +-+------- breaking API changes -- | | +----- non-breaking API additions -- | | | +--- code changes with no API change -version: 0.3.1.1 +version: 0.4.0.0 -- A short (one-line) description of the package. synopsis: Compute Maximum Entropy Distributions @@ -33,18 +33,12 @@ . Here is a the example from Probability the Logic of Science . - @ - maxent ([1,2,3], [average 1.5]) - @ - . + >>> 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]) - @ - . + >>> maxent ([1,2,3,4,5,6], [average 4.5]) Right [.05, .07, 0.11, 0.16, 0.23, 0.34] . One can use different constraints besides the average value there. @@ -78,22 +72,44 @@ build-type: Simple -- Constraint on the version of Cabal needed to build this package. -cabal-version: >=1.8 +cabal-version: >=1.12 library -- Modules exported by the library. - exposed-modules: MaxEnt + exposed-modules: Numeric.MaxEnt -- Modules included in this library but not exported. - other-modules: MaxEnt.Internal + other-modules: Numeric.MaxEnt.Internal, Numeric.MaxEnt.Linear, Numeric.MaxEnt.ConjugateGradient -- Other library packages from which modules are imported. build-depends: base ==4.6.*, nonlinear-optimization ==0.3.*, vector ==0.10.*, - ad ==3.2.* + ad ==3.4.* -- Directories containing source files. hs-source-dirs: src + default-language: Haskell2010 + +Test-Suite tests + hs-source-dirs: src, tests + type: exitcode-stdio-1.0 + main-is: Main.hs + build-depends: base ==4.6.*, + nonlinear-optimization ==0.3.*, + vector ==0.10.*, + ad ==3.4.*, + hmatrix ==0.14.*, + QuickCheck == 2.5.*, + test-framework-quickcheck2 ==0.3.*, + test-framework-quickcheck2 ==0.3.*, + test-framework-hunit ==0.3.*, + test-framework == 0.8.* + + default-language: Haskell2010 + + + + diff --git a/src/MaxEnt.hs b/src/MaxEnt.hs deleted file mode 100644 --- a/src/MaxEnt.hs +++ /dev/null @@ -1,55 +0,0 @@ --- | --- The maximum entropy method, or MAXENT, is variational approach for computing probability --- distributions given a list of moment, or expected value, constraints. --- --- Here are some links for background info. --- A good overview of applications: --- --- On the idea of maximum entropy in general: --- --- --- --- Use this package to compute discrete maximum entropy distributions over 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] --- --- One can use different constraints besides the average value there. --- --- As for why you want to maximize the entropy to find the probability constraint, --- I will say this for now. In the case of the average constraint --- it is a kin to choosing a integer partition with the most interger compositions. --- I doubt that makes any sense, but I will try to explain more with a blog post soon. --- -module MaxEnt ( - Constraint, - ExpectationFunction, - constraint, - average, - variance, - maxent, - maxentLinear, - generalMaxent -) where -import MaxEnt.Internal (Constraint, - ExpectationFunction, - constraint, - average, - variance, - maxent, - maxentLinear, - generalMaxent) diff --git a/src/MaxEnt/Internal.hs b/src/MaxEnt/Internal.hs deleted file mode 100644 --- a/src/MaxEnt/Internal.hs +++ /dev/null @@ -1,295 +0,0 @@ -{-# 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) -import Data.Traversable -import Numeric.AD.Types -import Numeric.AD.Internal.Classes ---import Numeric.AD.Lagrangian - -sumMap :: Num b => (a -> b) -> [a] -> b -sumMap f = sum . map f - -sumWith :: Num c => (a -> b -> c) -> [a] -> [b] -> c -sumWith f xs = sum . zipWith f xs - -pOfK :: Floating a => [a] -> [ExpectationFunction a] -> [a] -> Int -> a -pOfK values fs ls k = exp (negate . sumWith (\l f -> l * f (values !! k)) ls $ fs) / - partitionFunc values fs ls - -pOfKLinear :: Floating a => [[a]] -> [a] -> Int -> a -pOfKLinear matrix ls k = - exp (dot (matrix !! k) ls) / - partitionFuncLinear matrix ls - -probs :: (Floating b) - => [b] - -> [ExpectationFunction b] - -> [b] - -> [b] -probs values fs ls = map (pOfK values fs ls) [0..length values - 1] - -partitionFunc :: Floating a - => [a] - -> [ExpectationFunction a] - -> [a] - -> a -partitionFunc values fs ls = sum $ [ exp ((-l) * f x) | - x <- values, - (f, l) <- zip fs ls] - -objectiveFunc :: Floating a - => [a] - -> [ExpectationFunction a] - -> [a] - -> [a] - -> a -objectiveFunc values fs moments ls = log (partitionFunc values fs ls) - + sumWith (*) ls moments - - -dot :: Num a => [a] -> [a] -> a -dot x y = sum . zipWith (*) x $ y - -partitionFuncLinear :: Floating a - => [[a]] - -> [a] - -> a -partitionFuncLinear matrix ws = sum $ [ exp (dot as ws) | - as <- matrix] - --- This is almost the sam as the objectiveFunc -objectiveFuncLinear :: Floating a - => [a] - -> [[a]] - -> [a] - -> [a] - -> a -objectiveFuncLinear values as moments ls = - log (partitionFuncLinear as ls) - dot ls moments - - -linProbs :: (Floating b) - => [[b]] - -> [b] - -> [b] -linProbs matrix ls = - map (pOfKLinear matrix ls) [0..length matrix - 1] - - -entropy :: Floating a => [a] -> a -entropy = negate . sumMap (\p -> p * log p) - - --- This a functions that takes in a list of values and --- a list of probabilities -type GeneralConstraint a = [a] -> [a] -> a - - -lagrangian :: Floating a - => ([a], [GeneralConstraint a], [a]) - -> [a] - -> a -lagrangian (values, fs, constants) lamsAndProbs = result where - result = entropy ps + (sum $ zipWith (*) lams constraints) - constraints = zipWith (-) appliedConstraints constants - appliedConstraints = map (\f -> f values ps) fs - - -- split the args list - ps = take (length values) lamsAndProbs - lams = drop (length values) lamsAndProbs - -squaredGrad :: Num a - => (forall s. Mode s => [AD s a] -> AD s a) -> [a] -> a -squaredGrad f vs = sumMap (\x -> x*x) (grad f vs) - -generalObjectiveFunc :: Floating a => (forall b. Floating b => - ([b], [GeneralConstraint b], [b])) - -> [a] - -> a -generalObjectiveFunc params lamsAndProbs = - squaredGrad lang lamsAndProbs where - - lang :: Floating a => (forall s. Mode s => [AD s a] -> AD s a) - lang = lagrangian params - - -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. A function and the constant it equals. --- --- Think of it as the pair @(f, c)@ in the constraint --- --- @ --- Σ pₐ f(xₐ) = 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 Constraint a = (ExpectationFunction a, a) - --- | A function that takes an index and value and returns a value. --- See 'average' and 'variance' for examples. -type ExpectationFunction a = (a -> a) - --- make a constraint from function and constant -constraint :: Floating a => ExpectationFunction 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 - --- | Most general solver --- This will solve the langrangian by added the constraint that the --- probabilities must add up to zero. --- This is the slowest but most flexible method. - -generalMaxent :: (forall a. Floating a => ([a], [(GeneralConstraint a, 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 -generalMaxent params = result where - obj :: Floating a => [a] -> a - obj = generalObjectiveFunc objInput - - values :: Floating a => [a] - values = fst params - - constraints :: Floating a => [(GeneralConstraint a, a)] - constraints = snd params - - fsmoments :: Floating a => ([GeneralConstraint a], [a]) - fsmoments = unzip constraints - - --The new constraint is always needed for probability problems - - sumToOne :: Floating a => GeneralConstraint a - sumToOne _ xs = sumMap id xs - - gcons' :: Floating a => [GeneralConstraint a] - gcons' = fst fsmoments - - gcons :: Floating a => [GeneralConstraint a] - gcons = sumToOne : gcons' - - moments :: Floating a => [a] - moments = 1 : snd fsmoments - - - objInput :: Floating b => ([b], [GeneralConstraint b], [b]) - objInput = (values, gcons, moments) - - fs :: [[Double] -> [Double] -> Double] - fs = fst fsmoments - - --hmm maybe there is a better way to get rid of the defaulting - guess = U.fromList $ replicate - (length fs) (1.0 :: Double) - - result = case unsafePerformIO (optimize defaultParameters 0.00001 guess - (toFunction obj) - (toGradient obj) - Nothing) of - -- Not sure what is supposed to happen here - -- I guess I can plug the lambdas back in - (vs, ToleranceStatisfied, _) -> Right $ (S.toList vs) - (_, x, y) -> Left (x, y) - - - --- | This is for the linear case Ax = b --- @x@ in this situation is the vector of probablities. --- --- For example. --- --- @ --- maxentLinear ([1,1,1], ([[0.85, 0.1, 0.05], [0.25, 0.5, 0.25], [0.05, 0.1, 0.85]], [0.29, 0.42, 0.29])) --- @ --- --- Right [0.1, 0.8, 0.1] --- --- To be honest I am not sure why I can't use the 'maxent' version to solve --- this type of problem, but it doesn't work. I'm still learning --- -maxentLinear :: (forall a. Floating a => ([a], ([[a]], [a]))) -- ^ The values and a matrix A and column vector b - -> Either (Result, Statistics) [Double] -- ^ Either the a discription of what wrong or the probability distribution -maxentLinear params = result where - obj :: Floating a => [a] -> a - obj = uncurry (objectiveFuncLinear values) fsmoments - - values :: Floating a => [a] - values = fst params - - constraints :: Floating a => ([[a]], [a]) - constraints = snd params - - fsmoments :: Floating a => ([[a]], [a]) - fsmoments = constraints - - fs :: [[Double]] - fs = fst fsmoments - - -- hmm maybe there is a better way to get rid of the defaulting - guess = U.fromList $ replicate - (length fs) (2.0 :: Double) - - result = case unsafePerformIO (optimize defaultParameters 0.00001 guess - (toFunction obj) - (toGradient obj) - Nothing) of - (vs, ToleranceStatisfied, _) -> Right $ linProbs fs (S.toList vs) - (_, x, y) -> Left (x, y) - - --- | The main entry point for computing discrete maximum entropy distributions. --- Where the constraints are all moment constraints. -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 => [(ExpectationFunction a, a)] - constraints = snd params - - fsmoments :: Floating a => ([ExpectationFunction 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 fs) (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) - - - - - - - - diff --git a/src/Numeric/MaxEnt.hs b/src/Numeric/MaxEnt.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/MaxEnt.hs @@ -0,0 +1,60 @@ +-- | +-- The maximum entropy method, or MAXENT, is variational approach for computing probability +-- distributions given a list of moment, or expected value, constraints. +-- +-- Here are some links for background info. +-- A good overview of applications: +-- +-- On the idea of maximum entropy in general: +-- +-- +-- +-- Use this package to compute discrete maximum entropy distributions over 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] +-- +-- One can use different constraints besides the average value there. +-- +-- As for why you want to maximize the entropy to find the probability constraint, +-- I will say this for now. In the case of the average constraint +-- it is a kin to choosing a integer partition with the most interger compositions. +-- I doubt that makes any sense, but I will try to explain more with a blog post soon. +-- +module Numeric.MaxEnt ( + Constraint, + ExpectationFunction, + constraint, + average, + variance, + -- ** Classic moment based + maxent, + -- ** General + GeneralConstraint, + generalMaxent, + -- ** Linear + LinearConstraints(..), + linear, +) where +import Numeric.MaxEnt.Internal (Constraint, + ExpectationFunction, + constraint, + average, + variance, + maxent, + generalMaxent, + GeneralConstraint) +import Numeric.MaxEnt.Linear (linear, LinearConstraints(..)) + + + + + diff --git a/src/Numeric/MaxEnt/ConjugateGradient.hs b/src/Numeric/MaxEnt/ConjugateGradient.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/MaxEnt/ConjugateGradient.hs @@ -0,0 +1,50 @@ +{-# LANGUAGE TupleSections, Rank2Types #-} +module Numeric.MaxEnt.ConjugateGradient 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) + + +dot :: Num a => [a] -> [a] -> a +dot x y = sum . zipWith (*) x $ y + +sumMap :: Num b => (a -> b) -> [a] -> b +sumMap f = sum . map f + +sumWith :: Num c => (a -> b -> c) -> [a] -> [b] -> c +sumWith f xs = sum . zipWith f xs + + +toFunction :: (forall a. RealFloat a => [a] -> a) -> Function Simple +toFunction f = VFunction (f . U.toList) + +toGradient :: (forall a. RealFloat a => [a] -> a) -> Gradient Simple +toGradient f = VGradient (U.fromList . grad f . U.toList) + +toDoubleF :: (forall a. RealFloat a => [a] -> a) -> [Double] -> Double +toDoubleF f x = f x + +squaredGrad :: Num a + => (forall s. Mode s => [AD s a] -> AD s a) -> [a] -> a +squaredGrad f vs = sumMap (\x -> x*x) (grad f vs) + +solve :: Double + -> Int + -> (forall a. RealFloat a => [a] -> a) + -> Either (Result, Statistics) [Double] +solve percision count obj = result where + guess = U.fromList $ replicate + count ((1.0 :: Double) / (fromIntegral count)) + + result = case unsafePerformIO (optimize (defaultParameters { printFinal = False }) percision guess + (toFunction obj) + (toGradient obj) + Nothing) of + (vs, ToleranceStatisfied, _) -> Right . S.toList $ vs + (_, x, y) -> Left (x, y) diff --git a/src/Numeric/MaxEnt/Internal.hs b/src/Numeric/MaxEnt/Internal.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/MaxEnt/Internal.hs @@ -0,0 +1,210 @@ +{-# LANGUAGE TupleSections, Rank2Types #-} +module Numeric.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) +import Data.Traversable +import Numeric.AD.Types +import Numeric.AD.Internal.Classes +import Data.List (transpose) +--import Numeric.AD.Lagrangian +import Numeric.MaxEnt.ConjugateGradient + +pOfK :: RealFloat a => [a] -> [ExpectationFunction 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 :: (RealFloat b) + => [b] + -> [ExpectationFunction b] + -> [b] + -> [b] +probs values fs ls = map (pOfK values fs ls) [0..length values - 1] + +partitionFunc :: RealFloat a + => [a] + -> [ExpectationFunction a] + -> [a] + -> a +partitionFunc values fs ls = sum $ [ exp ((-l) * f x) | + x <- values, + (f, l) <- zip fs ls] + +objectiveFunc :: RealFloat a + => [a] + -> [ExpectationFunction a] + -> [a] + -> [a] + -> a +objectiveFunc values fs moments ls = log (partitionFunc values fs ls) + + sumWith (*) ls moments + +entropy :: RealFloat a => [a] -> a +entropy = negate . sumMap (\p -> p * log p) + + +-- This a functions that takes in a list of values and +-- a list of probabilities +type GeneralConstraint a = [a] -> [a] -> a + + +lagrangian :: RealFloat a + => ([a], [GeneralConstraint a], [a]) + -> [a] + -> a +lagrangian (values, fs, constants) lamsAndProbs = result where + result = entropy ps + (sum $ zipWith (*) lams constraints) + constraints = zipWith (-) appliedConstraints constants + appliedConstraints = map (\f -> f values ps) fs + + -- split the args list + ps = take (length values) lamsAndProbs + lams = drop (length values) lamsAndProbs + +generalObjectiveFunc :: RealFloat a => (forall b. RealFloat b => + ([b], [GeneralConstraint b], [b])) + -> [a] + -> a +generalObjectiveFunc params lamsAndProbs = + squaredGrad lang lamsAndProbs where + + lang :: RealFloat a => (forall s. Mode s => [AD s a] -> AD s a) + lang = lagrangian params + + + + +-- make a constraint from function and constant +constraint :: RealFloat a => ExpectationFunction a -> a -> Constraint a +constraint = (,) + +-- The average constraint +average :: RealFloat a => a -> Constraint a +average m = constraint id m + +-- The variance constraint +variance :: RealFloat a => a -> Constraint a +variance sigma = constraint (^(2 :: Int)) sigma + +-- | Most general solver +-- This is the slowest but most flexible method. Although, I haven't tried using much... + +generalMaxent :: (forall a. RealFloat a => ([a], [(GeneralConstraint a, 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 +generalMaxent params = result where + obj :: RealFloat a => [a] -> a + obj = generalObjectiveFunc objInput + + values :: RealFloat a => [a] + values = fst params + + constraints :: RealFloat a => [(GeneralConstraint a, a)] + constraints = snd params + + fsmoments :: RealFloat a => ([GeneralConstraint a], [a]) + fsmoments = unzip constraints + + --The new constraint is always needed for probability problems + + sumToOne :: RealFloat a => GeneralConstraint a + sumToOne _ xs = sumMap id xs + + gcons' :: RealFloat a => [GeneralConstraint a] + gcons' = fst fsmoments + + gcons :: RealFloat a => [GeneralConstraint a] + gcons = sumToOne : gcons' + + moments :: RealFloat a => [a] + moments = 1 : snd fsmoments + + + objInput :: RealFloat b => ([b], [GeneralConstraint b], [b]) + objInput = (values, gcons, moments) + + fs :: [[Double] -> [Double] -> Double] + fs = fst fsmoments + + --hmm maybe there is a better way to get rid of the defaulting + guess = U.fromList $ replicate + (length fs) (1.0 :: Double) + + result = case unsafePerformIO (optimize defaultParameters 0.00001 guess + (toFunction obj) + (toGradient obj) + Nothing) of + -- Not sure what is supposed to happen here + -- I guess I can plug the lambdas back in + (vs, ToleranceStatisfied, _) -> Right $ (S.toList vs) + (_, x, y) -> Left (x, y) + + +-- | Constraint type. A function and the constant it equals. +-- +-- Think of it as the pair @(f, c)@ in the constraint +-- +-- @ +-- Σ pₐ f(xₐ) = 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 Constraint a = (ExpectationFunction a, a) + +-- | A function that takes an index and value and returns a value. +-- See 'average' and 'variance' for examples. +type ExpectationFunction a = (a -> a) + +-- | Discrete maximum entropy solver where the constraints are all moment constraints. +maxent :: (forall a. RealFloat 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 :: RealFloat a => [a] -> a + obj = uncurry (objectiveFunc values) fsmoments + + values :: RealFloat a => [a] + values = fst params + + constraints :: RealFloat a => [(ExpectationFunction a, a)] + constraints = snd params + + fsmoments :: RealFloat a => ([ExpectationFunction 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 fs) (1.0 :: Double) + + {- + result = case unsafePerformIO (optimize (defaultParameters { printFinal = False }) + 0.001 guess + (toFunction obj) + (toGradient obj) + Nothing) of + (vs, ToleranceStatisfied, _) -> Right $ probs values fs (S.toList vs) + (_, x, y) -> Left (x, y) + -} + + result = case unsafePerformIO (optimize (defaultParameters { printFinal = False }) + 0.001 guess + (toFunction obj) + (toGradient obj) + Nothing) of + (vs, ToleranceStatisfied, _) -> Right $ probs values fs (S.toList vs) + (_, x, y) -> Left (x, y) + + + + + + + + diff --git a/src/Numeric/MaxEnt/Linear.hs b/src/Numeric/MaxEnt/Linear.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/MaxEnt/Linear.hs @@ -0,0 +1,115 @@ +{-# LANGUAGE TupleSections, Rank2Types #-} +module Numeric.MaxEnt.Linear where +import Numeric.MaxEnt.ConjugateGradient +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 + +probs :: (RealFloat b) => [[b]] -> [b] -> [b] +probs matrix ls = map (pOfK matrix ls) [0..length matrix - 1] + +pOfK :: RealFloat a => [[a]] -> [a] -> Int -> a +pOfK matrix ls k = exp (dot (transpose matrix !! k) ls) / partitionFunc matrix ls + +-- have the column and rows backwards +partitionFunc :: RealFloat a => [[a]] -> [a] -> a +partitionFunc matrix ws = sum $ [ exp (dot as ws) | as <- transpose matrix] + +-- This is almost the sam as the objectiveFunc +objectiveFunc :: RealFloat a => [[a]] -> [a]-> [a] -> a +objectiveFunc as moments ls = log (partitionFunc as ls) - dot ls moments + +data LinearConstraints a = LC { + matrix :: [[a]], + output :: [a] + } + deriving (Show, Eq) + +-- | This is for the linear case Ax = b +-- @x@ in this situation is the vector of probablities. +-- +-- Consider the 1 dimensional circular convolution using hmatrix. +-- +-- >>> import Numeric.LinearAlgebra +-- >>> fromLists [[0.68, 0.22, 0.1], [0.1, 0.68, 0.22], [0.22, 0.1, 0.68]] <> fromLists [[0.2], [0.5], [0.3]] +-- (3><1) [0.276, 0.426, 0.298] +-- +-- Now if we were given just the convolution and the output, we can use 'linear' to infer the input. +-- +-- >>> linear 3.0e-17 $ LC [[0.68, 0.22, 0.1], [0.1, 0.68, 0.22], [0.22, 0.1, 0.68]] [0.276, 0.426, 0.298] +-- Right [0.20000000000000004,0.4999999999999999,0.3] +-- +-- +-- +-- +linear :: Double -- ^ Tolerance for the numerical solver + -> (forall a. RealFloat a => LinearConstraints a) -- ^ The matrix A and column vector b + -> Either (Result, Statistics) [Double] -- ^ Either the a discription of what wrong or the probability distribution +linear tolerance constraints = result where + obj :: RealFloat a => [a] -> a + obj = objectiveFunc (matrix constraints) (output constraints) + + as :: [[Double]] + as = matrix constraints + + result = probs as <$> solve tolerance (length as) obj + +--Just for testing +linear1 :: Double -> (forall a. RealFloat a => LinearConstraints a) -- ^ a matrix A and column vector b + -> Either (Result, Statistics) [Double] -- ^ Either the a discription of what wrong or the probability distribution +linear1 precision constraints = result where + obj :: RealFloat a => [a] -> a + obj = objectiveFunc (matrix constraints) (output constraints) + + fs :: [[Double]] + fs = matrix constraints + + guess :: RealFloat a => [a] + guess = replicate + (length fs) ((1.0) / (fromIntegral (length fs))) + + result = Right . probs fs . last $ conjugateGradientDescent obj guess + + + +{- + I need to clean this up and submit it + submit the sharpen + + I need to figure out how to make linear2 + what is the + it might help to clean up the notation for the problem + basically I need to symbolically take the derative of the problem + and get rid of the automatic differentation + which I can do hand + + + + + + +-} + + + + + + + + + + + + + + + + + diff --git a/tests/Main.hs b/tests/Main.hs new file mode 100644 --- /dev/null +++ b/tests/Main.hs @@ -0,0 +1,18 @@ +module Main where +import Test.Framework (defaultMain, testGroup, defaultMainWithArgs) +import Test.Framework.Providers.HUnit +import Test.Framework.Providers.QuickCheck2 (testProperty) +import Numeric.MaxEnt +import LinearTesting + +-- One test is that a list convolved with a +-- guassian and then infered with a gaussian should have less entropy then it +-- started with + +main = defaultMain [ + testGroup "Linear Tests" [ + testProperty "solvableSystemsAreSolvable" solvableSystemsAreSolvable, + testProperty "probsSumToOne" probsSumToOne, + testProperty "solutionFitsConstraints" solutionFitsConstraints + ] + ]