diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,15 @@ Changelog ========= +Version 0.3.0.0 +--------------- + + + +* **(Breaking change)** Moved the top-level modules from `Data` to `Numeric`, + to better reflect the nature of the library and to align with the + convention of other similar libraries. + Version 0.2.0.0 --------------- @@ -9,5 +18,5 @@ * Initial release, re-written from the unreleased `0.1.0.0` by re-implementing error propagation with the [ad][] library. - [ad]: https://hackage.haskell.org/package/ad +[ad]: https://hackage.haskell.org/package/ad diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -12,7 +12,7 @@ ## Usage ```haskell -import Data.Uncertain +import Numeric.Uncertain ``` ### Create numbers @@ -67,7 +67,7 @@ samples. ```haskell -ghci> import Data.Uncertain.Correlated +ghci> import Numeric.Uncertain.Correlated ghci> evalCorr $ do x <- sampleUncert $ 12.5 +/- 0.8 y <- sampleUncert $ 15.9 +/- 0.5 @@ -99,7 +99,7 @@ [Monte Carlo simulations]: https://en.wikipedia.org/wiki/Monte_Carlo_method ```haskell -ghci> import qualified Data.Uncertain.MonteCarlo as MC +ghci> import qualified Numeric.Uncertain.MonteCarlo as MC ghci> import System.Random.MWC ghci> let x = 1.52 +/- 0.07 ghci> let y = 781.4 +/- 0.3 @@ -128,7 +128,7 @@ functions to the data with the uncertainties and properly propagating the errors with sound statistical principles. -[data-interval]: https://hackage.haskell.org/package/intervals +[intervals]: https://hackage.haskell.org/package/intervals [rounding]: https://hackage.haskell.org/package/rounding For a clear example, take diff --git a/src/Data/Hople.hs b/src/Data/Hople.hs --- a/src/Data/Hople.hs +++ b/src/Data/Hople.hs @@ -1,7 +1,7 @@ -{-# LANGUAGE CPP #-} {-# LANGUAGE DeriveFoldable #-} {-# LANGUAGE DeriveFunctor #-} {-# LANGUAGE DeriveTraversable #-} +{-# LANGUAGE NoImplicitPrelude #-} {-# OPTIONS_HADDOCK hide #-} {-# OPTIONS_HADDOCK prune #-} @@ -27,10 +27,7 @@ ) where -#if __GLASGOW_HASKELL__ < 710 -import Data.Foldable -import Data.Traversable -#endif +import Prelude.Compat data H1 a = H1 !a diff --git a/src/Data/Uncertain.hs b/src/Data/Uncertain.hs deleted file mode 100644 --- a/src/Data/Uncertain.hs +++ /dev/null @@ -1,603 +0,0 @@ -{-# LANGUAGE CPP #-} -{-# LANGUAGE DeriveDataTypeable #-} -{-# LANGUAGE DeriveGeneric #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE LambdaCase #-} -{-# LANGUAGE RankNTypes #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE ViewPatterns #-} - -#if __GLASGOW_HASKELL__ >= 708 -{-# LANGUAGE PatternSynonyms #-} -#endif - --- | --- Module : Data.Uncertain --- Copyright : (c) Justin Le 2016 --- License : BSD3 --- --- Maintainer : justin@jle.im --- Stability : experimental --- Portability : non-portable --- - -module Data.Uncertain - ( -- * 'Uncert' - Uncert -#if __GLASGOW_HASKELL__ >= 708 - , pattern (:+/-) -#endif - -- ** Creating 'Uncert' values - , (+/-), exact, withPrecision, withPrecisionAtBase, withVar, fromSamples - -- ** Inspecting properties - , uMean, uVar, uStd, uMeanVar, uMeanStd, uRange - -- * Applying arbitrary functions - , liftU - , liftU2, liftU3, liftU4, liftU5, liftUF - -- * Utility functions - , uNormalize, uNormalizeAtBase - , uShow, uShowsPrec - ) - where - -import Data.Data -import Data.Foldable (toList, foldl') -import Data.Function -import Data.Hople -import Data.Ord -import GHC.Generics -import Numeric.AD.Mode.Sparse -import qualified Numeric.AD.Mode.Tower as T - -#if __GLASGOW_HASKELL__ < 710 -import Data.Functor ((<$>)) -import Data.Traversable (Traversable) -#endif - - --- | Represents an independent experimental value centered around a mean --- value with "inherent" and independent uncertainty. --- --- Mostly useful due to its instances of numeric typeclasses like `Num`, --- `Fractional`, etc., which allows you to add and multiply and apply --- arbitrary numerical functions to them and have the uncertainty --- propagate appropriately. You can also lift arbitrary (sufficiently --- polymorphic) functions with 'liftU', 'liftUF', 'liftU2' and family. --- --- @ --- ghci> let x = 1.52 '+/-' 0.07 --- ghci> let y = 781.4 +/- 0.3 --- ghci> let z = 1.53e-1 `'withPrecision'` 3 --- ghci> cosh x --- 2.4 +/- 0.2 --- ghci> exp x / z * sin (y ** z) --- 10.9 +/- 0.9 --- ghci> pi + 3 * logBase x y --- 52 +/- 5 --- @ --- --- Uncertaintly is properly propagated according to the second-degree --- taylor series approximations of the applied functions. However, if the --- higher-degree terms are large with respect to to the means and --- variances of the uncertain values, these approximations may be --- inaccurate. --- --- Can be created with 'exact' to represent an "exact" measurement with no --- uncertainty, '+/-' and ':+/-' to specify a standard deviation as --- a range, 'withPrecision' to specify through decimal precision, and --- 'withVar' to specify with a variance. Can also be inferred from a list --- of samples with 'fromSamples' --- --- @ --- 7.13 '+/-' 0.05 --- 91800 +/- 100 --- 12.5 `'withVar'` 0.36 --- 'exact' 7.9512 --- 81.42 `'withPrecision'` 4 --- 7 :: Uncertain Double --- 9.18 :: Uncertain Double --- 'fromSamples' [12.5, 12.7, 12.6, 12.6, 12.5] --- @ --- --- Can be deconstructed with ':+/-', the pattern synonym/pseudo-constructor --- which matches on the mean and a standard deviation (supported on GHC --- 7.8+, with bidirectional constructor functionality supported on GHC --- 7.10+). You can also access properties with 'uMean', 'uStd', 'uVar', --- 'uMeanStd', 'uMeanVar', 'uRange', etc. --- --- It's important to remember that each "occurrence" represents a unique --- independent sample, so: --- --- @ --- ghci> let x = 15 '+/-' 2 in x + x --- 30 +/- 3 --- --- ghci> let x = 15 +/- 2 in x*2 --- 30 +/- 4 --- @ --- --- @x + x@ does not represent adding the same sample to itself twice, it --- represents /independently/ sampling two values within the range @15 +/- 2@ --- and adding them together. In general, errors and deviations will cancel --- each-other out, leading to a smaller uncertainty. --- --- However, @x*2@ represents taking /one/ sample and multiplying it by two. --- This yields a greater uncertainty, because errors and deviations are --- amplified. --- --- Also be aware that the 'Show' instance "normalizes" the result, and --- won't show any mean/central point to a decimal precision smaller than --- the uncertainty, rounding off the excess. --- -data Uncert a = Un { _uMean :: !a - , _uVar :: !a -- ^ maintained to be positive! - } - deriving (Data, Typeable, Generic, Generic1) - --- | Get the mean/central value/expected value of an 'Uncert'. -uMean :: Uncert a -> a -uMean = _uMean -{-# INLINE uMean #-} - --- | Get the /variance/ of the uncertainty of an 'Uncert', proportional to --- the square of "how uncertain" a value is. Is the square of 'uStd'. -uVar :: Uncert a -> a -uVar = _uVar -{-# INLINE uVar #-} - --- | Get the /standard deviation/ of the uncertainty of an 'Uncert', --- proportional to "how uncertain" a value is. --- --- Very informally, it can be thought of as the interval above and below --- the mean that about 68% of sampled values will fall under after repeated --- sampling, or as the range that one is 68% sure the true value is within. --- --- Is the square root of 'uVar'. -uStd :: Floating a => Uncert a -> a -uStd = sqrt . uVar -{-# INLINE uStd #-} - --- | Create an 'Uncert' with an exact value and 0 uncertainty. -exact - :: Num a - => a -- ^ The exact value - -> Uncert a -exact x = Un x 0 -{-# INLINE exact #-} - -infixl 6 +/- -#if __GLASGOW_HASKELL__ >= 708 -infixl 6 :+/- -#endif - --- | Create an 'Uncert' around a central value and a given "range" of --- uncertainty. The range is interpreted as the standard deviation of the --- underlying random variable. Might be preferrable over ':+/-' because it --- is more general (doesn't require a 'Floating' constraint) and looks --- a bit nicer. --- --- See 'uStd' for more details. -(+/-) - :: Num a - => a -- ^ The mean or central value - -> a -- ^ The standard deviation of the underlying uncertainty - -> Uncert a -x +/- dx = Un x (dx*dx) -{-# INLINE (+/-) #-} - --- | Create an 'Uncert' around a central value, specifying its uncertainty --- with a given /variance/. The variance is taken to be proportional to --- the square of the range of uncertainty. See 'uStd' for more details. --- --- "Negative variances" are treated as positive. -withVar - :: Num a - => a -- ^ The mean or central value - -> a -- ^ The variance of the underlying uncertainty - -> Uncert a -withVar x vx = Un x (abs vx) -{-# INLINE withVar #-} - -#if __GLASGOW_HASKELL__ >= 708 --- | Pattern match on an 'Uncert' with its central value and its standard --- deviation (see 'uStd' for clarification). --- --- Can also be used to /construct/ an 'Uncert', identically as '+/-'. --- --- /Note:/ Only supported on GHC 7.8 and above. Bidirectional --- functionality (to allow use as a constructor) only supported on GHC --- 7.10 and above. --- -#if __GLASGOW_HASKELL__ >= 710 -pattern (:+/-) :: () => Floating a => a -> a -> Uncert a -#endif -pattern x :+/- dx <- Un x (sqrt->dx) -#if __GLASGOW_HASKELL__ >= 710 - where - x :+/- dx = Un x (dx*dx) -#endif -#endif - --- | Infer an 'Uncert' from a given list of independent /samples/ of an --- underlying uncertain or random distribution. -fromSamples :: Fractional a => [a] -> Uncert a -fromSamples = makeUn . foldStats - where - makeUn (H3 x0 x1 x2) = Un μ v - where - μ = x1/x0 - v = x2/x0 - μ*μ -- maybe use pop var? - foldStats = flip foldl' (H3 0 0 0) $ - \(H3 s0 s1 s2) x -> - H3 (s0 + 1) (s1 + x) (s2 + x*x) -{-# INLINABLE fromSamples #-} - --- | Retrieve both the mean (central) value and the underlying variance of --- an 'Uncert' together. --- --- @uMeanVar ≡ 'uMean' &&& 'uVar'@ -uMeanVar :: Uncert a -> (a, a) -uMeanVar (Un x vx) = (x, vx) -{-# INLINE uMeanVar #-} - --- | Retreve both the mean (central) value and the underlying standard --- deviation of an 'Uncert' together. (See 'uStd' for more details) --- --- @uMeanStd ≡ 'uMean' &&& 'uStd'@ -uMeanStd :: Floating a => Uncert a -> (a, a) -uMeanStd (Un x vx) = (x, sqrt vx) -{-# INLINE uMeanStd #-} - --- | Retrieve the "range" of the underlying distribution of an 'Uncert', --- derived from the standard deviation, where which approximly 68% of --- sampled values are expected to occur (or within which you are 68% --- certain the true value is). --- --- @uRange (x +/- dx) ≡ (x - dx, x + dx)@ -uRange :: Floating a => Uncert a -> (a, a) -uRange (uMeanStd->(x, dx)) = (x - dx, x + dx) -{-# INLINABLE uRange #-} - --- | Like 'withPrecision', except takes a number of "digits" of precision in --- the desired numeric base. For example, in base 2, takes the number of --- /bits/ of precision. --- --- @'withPrecision' ≡ withPrecisionAtBase 10@ -withPrecisionAtBase - :: (Floating a, RealFrac a) - => Int -- ^ The base to determine precision with respect to - -> a -- ^ The approximate value of the 'Uncert' - -> Int -- ^ The number of "digits" of precision to take - -> Uncert a -withPrecisionAtBase b x p = x' +/- dx' - where - leading :: Int - leading = negate . floor . logBase (fromIntegral b) $ x - uncert :: Int - uncert = leading - 1 + fromIntegral p - rounder = fromIntegral b ** fromIntegral uncert - x' = (/ rounder) . fromIntegral . round' . (* rounder) $ x - dx' = 1 / rounder - round' :: RealFrac a => a -> Integer - round' = round -{-# INLINABLE withPrecisionAtBase #-} - --- | Create an 'Uncert' about a given approximate central value, with the --- given number of /digits of precision/ (in decimal notation). --- --- @5.21 `withPrecision` 3 ≡ 5.21 '+/-' 0.01@ -withPrecision - :: (Floating a, RealFrac a) - => a -- ^ The approximate value of the 'Uncert' - -> Int -- ^ The number of "digits" of precision to take - -> Uncert a -withPrecision = withPrecisionAtBase 10 -{-# INLINABLE withPrecision #-} - --- | Like 'uNormalize', but takes a numerical base to round with respect --- to. --- --- @'uNormalize' ≡ uNormalizeAtBase 10@ -uNormalizeAtBase - :: (Floating a, RealFrac a) - => Int -- ^ The base to normalize with respect to - -> Uncert a - -> Uncert a -uNormalizeAtBase b (uMeanStd->(x, dx)) = x' +/- dx' - where - uncert :: Int - uncert = negate . floor . logBase (fromIntegral b) $ dx - rounder = fromIntegral b ** fromIntegral uncert - roundTo = (/ rounder) . fromIntegral . round' . (* rounder) - x' = roundTo x - dx' = roundTo dx - round' :: RealFrac a => a -> Integer - round' = round -{-# INLINABLE uNormalizeAtBase #-} - --- | Attempts to "normalize" an 'Uncert'. Rounds the uncertainty (the --- standard deviation) to one digit of precision, and rounds the central --- moment up to the implied precision. --- --- For example, it makes no real sense to have @542.185433 +/- 83.584@, --- because the extra digits of @542.185433@ past the tens place has no --- meaning because of the overpowering uncertainty. Normalizing this --- results in @540 +/- 80@. --- --- Note that the 'Show' instance for 'Uncert' normalizes values before --- showing them. -uNormalize - :: (Floating a, RealFrac a) - => Uncert a - -> Uncert a -uNormalize = uNormalizeAtBase 10 -{-# INLINABLE uNormalize #-} - -instance (Show a, Floating a, RealFrac a) => Show (Uncert a) where - showsPrec d = uShowsPrec d . uNormalize - --- | Like 'showsPrec' for 'Uncert', but does not normalize the value (see --- 'uNormalize') before showing. See documentation for 'showsPrec' for --- more information on how this is meant to be used. -uShowsPrec :: (Show a, Floating a) => Int -> Uncert a -> ShowS -uShowsPrec d (uMeanStd->(x, dx)) = showParen (d > 5) $ - showsPrec 6 x - . showString " +/- " - . showsPrec 6 dx -{-# INLINABLE uShowsPrec #-} - --- | Like 'show' for 'Uncert', but does not normalize the value (see --- 'uNormalize') before showing. --- --- @'show' ≡ uShow . 'uNormalize'@ -uShow :: (Show a, Floating a) => Uncert a -> String -uShow u = uShowsPrec 0 u "" -{-# INLINABLE uShow #-} - --- | Lifts a multivariate numeric function on a container (given as an @f --- a -> a@) to work on a container of 'Uncert's. Correctly propagates the --- uncertainty according to the second-order (multivariate) taylor --- expansion of the function. Note that if the higher-degree taylor series --- terms are large with respect to the means and variances, this --- approximation may be inaccurate. --- --- Should take any function sufficiently polymorphic over numeric types, so --- you can use things like '*', 'sqrt', 'atan2', etc. --- --- @ --- ghci> liftUF (\[x,y,z] -> x*y+z) [12.2 +/- 0.5, 56 +/- 2, 0.12 +/- 0.08] --- 680 +/- 40 --- @ --- -liftUF - :: (Traversable f, Fractional a) - => (forall s. f (AD s (Sparse a)) -> AD s (Sparse a)) -- ^ Function on container of values to lift - -> f (Uncert a) -- ^ Container of 'Uncert's to apply the function to - -> Uncert a -liftUF f us = Un y vy - where - xs = uMean <$> us - vxs = uVar <$> us - vxsL = toList vxs - (fx, dfxsh) = hessian' f xs - dfxs = fst <$> dfxsh - hess = snd <$> dfxsh - y = fx + hessQuad / 2 - where - hessQuad = dot vxsL - . diag - . toList - $ fmap toList hess - vy = dot vxsL ((^ (2::Int)) <$> dfxs) - dot x = sum . zipWith (*) x . toList - diag = \case [] -> [] - [] :yss -> diag (drop1 <$> yss) - (x:_):yss -> x : diag (drop1 <$> yss) - where - drop1 [] = [] - drop1 (_:zs) = zs -{-# INLINABLE liftUF #-} - --- | Lifts a numeric function over an 'Uncert'. Correctly propagates the --- uncertainty according to the second-order taylor expansion expansion of --- the function. Note that if the higher-degree taylor series terms are --- large with respect to the mean and variance, this approximation may be --- inaccurate. --- --- Should take any function sufficiently polymorphic over numeric types, so --- you can use things like 'sqrt', 'sin', 'negate', etc. --- --- @ --- ghci> liftU (\x -> log x ^ 2) (12.2 +/- 0.5) --- 6.3 +/- 0.2 --- @ -liftU - :: Fractional a - => (forall s. AD s (T.Tower a) -> AD s (T.Tower a)) -- ^ Function on values to lift - -> Uncert a -- ^ 'Uncert' to apply the function to - -> Uncert a -liftU f (Un x vx) = Un y vy - where - fx:dfx:ddfx:_ = T.diffs0 f x - y = fx + ddfx * vx / 2 - vy = dfx*dfx * vx -{-# INLINABLE liftU #-} - --- | Lifts a two-argument (curried) function over two 'Uncert's. Correctly --- propagates the uncertainty according to the second-order (multivariate) --- taylor expansion expansion of the function. Note that if the --- higher-degree taylor series terms are large with respect to the mean and --- variance, this approximation may be inaccurate. --- --- Should take any function sufficiently polymorphic over numeric types, so --- you can use things like '*', 'atan2', '**', etc. --- --- @ --- ghci> liftU2 (\x y -> x**y) (13.5 +/- 0.1) (1.64 +/- 0.08) --- 70 +/- 10 --- @ -liftU2 - :: Fractional a - => (forall s. AD s (Sparse a) -> AD s (Sparse a) -> AD s (Sparse a)) - -> Uncert a - -> Uncert a - -> Uncert a -liftU2 f = curryH2 $ liftUF (uncurryH2 f) -{-# INLINABLE liftU2 #-} - --- | Lifts a three-argument (curried) function over three 'Uncert's. See --- 'liftU2' and 'liftUF' for more details. -liftU3 - :: Fractional a - => (forall s. AD s (Sparse a) -> AD s (Sparse a) -> AD s (Sparse a) -> AD s (Sparse a)) - -> Uncert a - -> Uncert a - -> Uncert a - -> Uncert a -liftU3 f = curryH3 $ liftUF (uncurryH3 f) -{-# INLINABLE liftU3 #-} - --- | Lifts a four-argument (curried) function over four 'Uncert's. See --- 'liftU2' and 'liftUF' for more details. -liftU4 - :: Fractional a - => (forall s. AD s (Sparse a) -> AD s (Sparse a) -> AD s (Sparse a) -> AD s (Sparse a) -> AD s (Sparse a)) - -> Uncert a - -> Uncert a - -> Uncert a - -> Uncert a - -> Uncert a -liftU4 f = curryH4 $ liftUF (uncurryH4 f) -{-# INLINABLE liftU4 #-} - --- | Lifts a five-argument (curried) function over five 'Uncert's. See --- 'liftU2' and 'liftUF' for more details. -liftU5 - :: Fractional a - => (forall s. AD s (Sparse a) -> AD s (Sparse a) -> AD s (Sparse a) -> AD s (Sparse a) -> AD s (Sparse a) -> AD s (Sparse a)) - -> Uncert a - -> Uncert a - -> Uncert a - -> Uncert a - -> Uncert a - -> Uncert a -liftU5 f = curryH5 $ liftUF (uncurryH5 f) -{-# INLINABLE liftU5 #-} - -instance Fractional a => Num (Uncert a) where - (+) = liftU2 (+) - {-# INLINE (+) #-} - (*) = liftU2 (*) - {-# INLINE (*) #-} - (-) = liftU2 (-) - {-# INLINE (-) #-} - negate = liftU negate - {-# INLINE negate #-} - abs = liftU abs - {-# INLINE abs #-} - signum = liftU signum - {-# INLINE signum #-} - fromInteger = exact . fromInteger - {-# INLINE fromInteger #-} - -instance Fractional a => Fractional (Uncert a) where - recip = liftU recip - {-# INLINE recip #-} - (/) = liftU2 (/) - {-# INLINE (/) #-} - fromRational = exact . fromRational - {-# INLINE fromRational #-} - -instance Floating a => Floating (Uncert a) where - pi = exact pi - {-# INLINE pi #-} - exp = liftU exp - {-# INLINE exp #-} - log = liftU log - {-# INLINE log #-} - sqrt = liftU sqrt - {-# INLINE sqrt #-} - (**) = liftU2 (**) - {-# INLINE (**) #-} - logBase = liftU2 logBase - {-# INLINE logBase #-} - sin = liftU sin - {-# INLINE sin #-} - cos = liftU cos - {-# INLINE cos #-} - asin = liftU asin - {-# INLINE asin #-} - acos = liftU acos - {-# INLINE acos #-} - atan = liftU atan - {-# INLINE atan #-} - sinh = liftU sinh - {-# INLINE sinh #-} - cosh = liftU cosh - {-# INLINE cosh #-} - asinh = liftU asinh - {-# INLINE asinh #-} - acosh = liftU acosh - {-# INLINE acosh #-} - atanh = liftU atanh - {-# INLINE atanh #-} - -instance Eq a => Eq (Uncert a) where - (==) = (==) `on` uMean - {-# INLINE (==) #-} - (/=) = (/=) `on` uMean - {-# INLINE (/=) #-} - -instance Ord a => Ord (Uncert a) where - compare = comparing uMean - {-# INLINE compare #-} - -instance (Fractional a, Real a) => Real (Uncert a) where - toRational = toRational . uMean - {-# INLINE toRational #-} - -instance RealFrac a => RealFrac (Uncert a) where - properFraction x = (n, d) - where - d = liftU (snd' . properFraction) x - n = fst . properFraction $ uMean x - snd' :: (Int, b) -> b - snd' = snd - {-# INLINABLE properFraction #-} - truncate = truncate . uMean - {-# INLINE truncate #-} - round = round . uMean - {-# INLINE round #-} - ceiling = ceiling . uMean - {-# INLINE ceiling #-} - floor = floor . uMean - {-# INLINE floor #-} - -instance RealFloat a => RealFloat (Uncert a) where - floatRadix = floatRadix . uMean - {-# INLINE floatRadix #-} - floatDigits = floatDigits . uMean - {-# INLINE floatDigits #-} - floatRange = floatRange . uMean - {-# INLINE floatRange #-} - decodeFloat = decodeFloat . uMean - {-# INLINE decodeFloat #-} - exponent = exponent . uMean - {-# INLINE exponent #-} - isNaN = isNaN . uMean - {-# INLINE isNaN #-} - isInfinite = isInfinite . uMean - {-# INLINE isInfinite #-} - isDenormalized = isDenormalized . uMean - {-# INLINE isDenormalized #-} - isNegativeZero = isNegativeZero . uMean - {-# INLINE isNegativeZero #-} - isIEEE = isIEEE . uMean - {-# INLINE isIEEE #-} - encodeFloat a b = exact (encodeFloat a b) - {-# INLINE encodeFloat #-} - significand = liftU significand - {-# INLINE significand #-} - atan2 = liftU2 atan2 - {-# INLINE atan2 #-} - diff --git a/src/Data/Uncertain/Correlated.hs b/src/Data/Uncertain/Correlated.hs deleted file mode 100644 --- a/src/Data/Uncertain/Correlated.hs +++ /dev/null @@ -1,90 +0,0 @@ -{-# LANGUAGE RankNTypes #-} - --- | --- Module : Data.Uncertain.Correlated --- Copyright : (c) Justin Le 2016 --- License : BSD3 --- --- Maintainer : justin@jle.im --- Stability : experimental --- Portability : non-portable --- --- Provides the 'Corr' monad, which allows one to describe complex --- relationships between random variables and evaluate their propagated --- uncertainties /respecting/ their inter-correlations. --- --- See the "Data.Uncertain.Correlated.Interactive" module for an --- "interactive" and exploratory interface for this module's functionality. --- - -module Data.Uncertain.Correlated - ( -- * 'Corr' - Corr, evalCorr - -- * Uncertain and Correlated Values - , CVar - -- ** Sampling - , sampleUncert, sampleExact, constC - -- ** Resolving - , resolveUncert - -- * Applying arbitrary functions - , liftC, liftC2, liftC3, liftC4, liftC5, liftCF - ) - where - -import Control.Monad.Free -import Control.Monad.Trans.State -import Data.Uncertain -import Data.Uncertain.Correlated.Internal -import qualified Data.IntMap.Strict as M - --- | Evaluates the value described by a 'Corr' monad, taking into account --- inter-correlations between samples. --- --- Takes a universally qualified 'Corr', which should not affect usage. --- See the examples in the documentation for 'Corr'. The univeral --- qualification is mostly a type system trick to ensure that you aren't --- allowed to ever use 'evalCorr' to evaluate a 'CVar'. -evalCorr :: Fractional a => (forall s. Corr s a b) -> b -evalCorr c = evalState (corrToState c) (0, M.empty) -{-# INLINABLE evalCorr #-} - --- | Generate a sample in 'Corr' from an 'Uncert' value, independently from --- all other samples. --- --- Note that you can only sample @'Uncert' a@s within a @'Corr' s a@, meaning --- that all other "sampled" values are also @a@s. -sampleUncert :: Uncert a -> Corr s a (CVar s a) -sampleUncert u = liftF $ Gen u id -{-# INLINE sampleUncert #-} - --- | Generate an exact sample in 'Corr' with zero uncertainty, --- independently from all other samples. --- --- Not super useful, since you can do something equivalent with 'constC' --- or the numeric instances: --- --- @ --- sampleExact x ≡ return ('constC' x) --- sampleExact 10 ≡ return 10 --- @ --- --- But is provided for completeness alongside 'sampleUncert'. --- --- Note that you can exactly sample an @a@ within a @'Corr' s a@, meaning --- that all other "sampled" values are also @a@s. --- -sampleExact :: a -> Corr s a (CVar s a) -sampleExact = return . constC -{-# INLINE sampleExact #-} - --- | "Resolve" an 'Uncert' from a 'CVar' using its potential multiple --- samples and sample sources, taking into account inter-correlations --- between 'CVar's and samples. --- --- Note that if you use 'sampleUncert' on the result, the new sample will --- be treated as something completely independent. Usually this should --- only be used as the "exit point" of a 'Corr' description. -resolveUncert :: CVar s a -> Corr s a (Uncert a) -resolveUncert v = liftF $ Rei v id -{-# INLINE resolveUncert #-} - diff --git a/src/Data/Uncertain/Correlated/Interactive.hs b/src/Data/Uncertain/Correlated/Interactive.hs deleted file mode 100644 --- a/src/Data/Uncertain/Correlated/Interactive.hs +++ /dev/null @@ -1,119 +0,0 @@ --- | --- Module : Data.Uncertain.Correlated.Interactive --- Copyright : (c) Justin Le 2016 --- License : BSD3 --- --- Maintainer : justin@jle.im --- Stability : experimental --- Portability : non-portable --- --- Exports all of the interface of "Data.Uncertain.Correlated", except --- meant to be run in a /ghci/ session "interactively" for exploratory --- purposes, or in a plain 'IO' action (instead of inside a 'Corr' monad). --- --- For example, with the "Data.Uncertain.Correlated" interface: --- --- @ --- ghci> evalCorr $ do --- x <- sampleUncert $ 12.5 +/- 0.8 --- y <- sampleUncert $ 15.9 +/- 0.5 --- z <- sampleUncert $ 1.52 +/- 0.07 --- let k = y**x --- resolveUncert $ (x+z) * logBase z k --- 1200 +/- 200 --- @ --- --- And with the interface from this "interactive" module: --- --- @ --- ghci> x <- 'sampleUncert' $ 12.5 +/- 0.8 --- ghci> y <- sampleUncert $ 15.9 +/- 0.5 --- ghci> z <- sampleUncert $ 1.52 +/- 0.07 --- ghci> let k = y**x --- ghci> 'resolveUncert' $ (x+z) * logBase z k --- 1200 +/- 200 --- @ --- --- The main purpose of this module is to allow one to use /ghci/ as a fancy --- "calculator" for computing and exploring propagated uncertainties of --- complex and potentially correlated samples with uncertainty. --- --- Because many of the names overlap with the names from the --- "Data.Uncertain.Correlated" module, it is recommended that you never --- have both imported at the same time in /ghci/ or in a file, or import --- them qualified if you must. --- --- Also note that all of these methods only work with @'Uncertain' --- 'Double'@s, and are not polymorphic over different numeric types. --- --- Be aware that this module is not robustly tested in heavily concurrent --- situations/applications. --- -module Data.Uncertain.Correlated.Interactive - ( -- * Uncertain and Correlated Values - CVar, CVarIO - -- ** Sampling - , sampleUncert, sampleExact, constC - -- ** Resolving - , resolveUncert - -- * Applying arbitrary functions - , liftC, liftC2, liftC3, liftC4, liftC5, liftCF - ) - where - -import Control.Monad.ST -import Control.Monad.Trans.State -import Data.IORef -import Data.Tuple -import Data.Uncertain -import Data.Uncertain.Correlated.Internal -import System.IO.Unsafe (unsafePerformIO) -import qualified Data.IntMap.Strict as M -import qualified Data.Uncertain.Correlated as C - --- | A 'CVar' specialized to work in an "interactive" context, in /ghci/ or --- 'IO'. -type CVarIO = CVar RealWorld Double - --- ssh, don't tell anyone we're using 'unsafePerformIO' -globalCorrMap :: IORef (M.Key, M.IntMap (Uncert Double)) -{-# NOINLINE globalCorrMap #-} -globalCorrMap = unsafePerformIO $ newIORef (0, M.empty) - -runCorrIO :: Corr RealWorld Double a -> IO a -runCorrIO c = atomicModifyIORef' globalCorrMap - (swap . runState (corrToState c)) -{-# INLINE runCorrIO #-} - --- | Generate a sample in 'IO' from an @'Uncert' 'Double'@ value, --- independently from all other samples. -sampleUncert :: Uncert Double -> IO CVarIO -sampleUncert u = runCorrIO $ C.sampleUncert u -{-# INLINABLE sampleUncert #-} - --- | Generate an exact sample in 'IO' with zero uncertainty, --- independently from all other samples. --- --- Not super useful, since you can do something equivalent with 'constC' --- or the numeric instances: --- --- @ --- sampleExact x ≡ return ('constC' x) --- sampleExact 10 ≡ return 10 --- @ --- --- But is provided for completeness alongside 'sampleUncert'. -sampleExact :: Double -> IO CVarIO -sampleExact d = runCorrIO $ C.sampleExact d -{-# INLINABLE sampleExact #-} - --- | "Resolve" an 'Uncert' from a 'CVarIO' using its potential multiple --- samples and sample sources, taking into account inter-correlations --- between 'CVarIO's and samples. --- --- Note that if you use 'sampleUncert' on the result, the new sample will --- be treated as something completely independent. Usually this should --- only be used as the "final value" of your computation or exploration. -resolveUncert :: CVarIO -> IO (Uncert Double) -resolveUncert v = runCorrIO $ C.resolveUncert v -{-# INLINABLE resolveUncert #-} diff --git a/src/Data/Uncertain/Correlated/Internal.hs b/src/Data/Uncertain/Correlated/Internal.hs deleted file mode 100644 --- a/src/Data/Uncertain/Correlated/Internal.hs +++ /dev/null @@ -1,374 +0,0 @@ -{-# LANGUAGE CPP #-} -{-# LANGUAGE DeriveFunctor #-} -{-# LANGUAGE GADTs #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} -{-# LANGUAGE KindSignatures #-} -{-# LANGUAGE LambdaCase #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE RankNTypes #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE StandaloneDeriving #-} -{-# OPTIONS_HADDOCK hide #-} -{-# OPTIONS_HADDOCK prune #-} - --- | --- Module : Data.Uncertain.Correlated.Internal --- Copyright : (c) Justin Le 2016 --- License : BSD3 --- --- Maintainer : justin@jle.im --- Stability : experimental --- Portability : non-portable --- --- Internal utility functions for functionality shared by --- "Data.Uncertain.Correlated" and "Data.Uncertain.Correlated.Interactive". --- - -module Data.Uncertain.Correlated.Internal - ( CVar, dephantom - , CorrF(..), Corr - , liftCF - , constC, liftC, liftC2, liftC3, liftC4, liftC5 - , corrToState - ) - where - -import Control.Arrow ((***)) -import Control.Monad.Free -import Control.Monad.Trans.State -import Data.Hople -import Data.Uncertain -import Numeric.AD.Mode.Sparse -import qualified Data.IntMap.Strict as M - -#if __GLASGOW_HASKELL__ < 710 -import Control.Applicative (Applicative) -import Data.Functor ((<$>)) -#endif - --- | Represents a single sample (or a value calculated from samples) within --- the 'Corr' monad. These can be created with 'sampleUncert', --- 'sampleExact', and 'constC', or made by combinining others with its --- numeric typeclass instances (like 'Num') or its functions lifting --- arbitrary numeric functions (like 'liftC2'). These keep track of --- inter-correlations between sources, and if you add together two 'CVar's --- that are correlated, their results will reflect this. --- --- Can be "resolved" into the uncertain value they represent using --- 'resolveUncert'. --- --- Note that these are parameterized by a dummy phantom parameter 's' so --- that they can't be "evaluated" out of the 'Corr' they live in with --- 'evalCorr'. --- --- Note that a @'CVar' s a@ can only ever meaningfully "exist" in a @'Corr' --- s a@, meaning that the all samples within that 'Corr' are of the same --- type. -data CVar s a where - CK :: a -> CVar s a - CV :: M.Key -> CVar s a - CF :: Functor f - => (forall t. f (AD t (Sparse a)) -> AD t (Sparse a)) - -> f (CVar s a) - -> CVar s a - --- | Unsafe function to bypass the universal qualification guard for --- returning 'CVar's from 'Corr's. -dephantom :: CVar s a -> CVar t a -dephantom = \case CK x -> CK x - CV k -> CV k - CF f xs -> CF f (dephantom <$> xs) - -data CorrF :: * -> * -> * -> * where - Gen :: Uncert a -> (CVar s a -> b) -> CorrF s a b - Fun :: Functor f - => (forall t. f (AD t (Sparse a)) -> AD t (Sparse a)) - -> f (CVar s a) - -> (CVar s a -> b) - -> CorrF s a b - Rei :: CVar s a - -> (Uncert a -> b) - -> CorrF s a b - -instance Functor (CorrF s a) where - fmap f = \case Gen u next -> Gen u (f . next) - Fun g us next -> Fun g us (f . next) - Rei v next -> Rei v (f . next) - - --- | The 'Corr' monad allows us to keep track of correlated and --- non-independent samples. It fixes a basic "failure" of the 'Uncert' --- type, which can't describe correlated samples. --- --- For example, consider the difference between: --- --- @ --- ghci> sum $ replicate 10 (12.5 '+/-' 0.8) --- 125 +/- 3 --- ghci> 10 * (12.5 +/- 0.8) --- 125 +/- 8 --- @ --- --- The first one represents the addition of ten independent samples, whose --- errors will in general cancel eachother out. The second one represents --- sampling once and multiplying it by ten, which will amplify any error by --- a full factor of 10. --- --- See how the 'Corr' monad expresses the above computations: --- --- @ --- ghci> 'evalCorr' $ do --- x <- 'sampleUncert' $ 12.5 '+/-' 0.8 --- y1 <- 'resolveUncert' $ sum (replicate 10 x) --- y2 <- resolveUncert $ 10 * x --- return (y1, y2) --- (125 +/- 8, 125 +/- 8) --- --- ghci> 'evalCorr' $ do --- xs <- replicateM 10 ('sampleUncert' (12.5 +/- 0.8)) --- 'resolveUncert' $ sum xs --- 125 +/- 3 --- @ --- --- The first example samples once and describes operations on the single --- sample; the second example samples 10 times with 'replicateM' and sums --- all of the results. --- --- Things are more interesting when you sample multiple variables: --- --- @ --- ghci> 'evalCorr' $ do --- x <- 'sampleUncert' $ 12.5 '+/-' 0.8 --- y <- sampleUncert $ 15.9 +/- 0.5 --- z <- sampleUncert $ 1.52 +/- 0.07 --- let k = y ** x --- 'resolveUncert' $ (x+z) * logBase z k --- 1200 +/- 200 --- @ --- --- The first parameter is a dummy phantom parameter used to prevent 'CVar's --- from leaking out of the computation (see 'evalCorr'). The second --- parameter is the numeric type of all samples within the description (for --- example, if you ever sample an 'Uncert Double', the second parameter wil --- be 'Double'). The third parameter is the result type of the --- computation -- the value the 'Corr' is describing. -newtype Corr s a b = Corr { corrFree :: Free (CorrF s a) b - } - deriving (Functor, Applicative, Monad) - -deriving instance MonadFree (CorrF s a) (Corr s a) - -corrToState - :: (Monad m, Fractional a) - => Corr s a b - -> StateT (M.Key, M.IntMap (Uncert a)) m b -corrToState = iterM go . corrFree - where - go = \case - Gen u next -> do - i <- gets fst - modify $ succ *** M.insert i u - next (CV i) - Fun f us next -> - next $ CF f us - Rei v next -> do - u <- gets (getCVar v . snd) - next u - getCVar - :: forall a s. Fractional a - => CVar s a - -> M.IntMap (Uncert a) - -> Uncert a - getCVar cv = liftUF (cVarToF cv) - where - cVarToF - :: CVar s a - -> (forall t. M.IntMap (AD t (Sparse a)) -> AD t (Sparse a)) - cVarToF (CK x) _ = auto x - cVarToF (CV k) us = us M.! k - cVarToF (CF f cs) us = f (flip cVarToF us <$> cs) -{-# INLINABLE corrToState #-} - --- | Lifts a multivariate numeric function on a container (given as an @f --- a -> a@) to work on a container of 'CVar's. Correctly propagates the --- uncertainty according to the second-order (multivariate) taylor --- expansion of the function, and properly takes into account and keeps --- track of all inter-correlations between the 'CVar' samples. Note that --- if the higher-degree taylor series terms are large with respect to the --- means and variances, this approximation may be inaccurate. --- --- Should take any function sufficiently polymorphic over numeric types, so --- you can use things like '*', 'sqrt', 'atan2', etc. --- --- @ --- ghci> evalCorr $ do --- x <- sampleUncert $ 12.5 +/- 0.8 --- y <- sampleUncert $ 15.9 +/- 0.5 --- z <- sampleUncert $ 1.52 +/- 0.07 --- resolveUncert $ liftCF (\[a,b,c] -> (a+c) * logBase c (b**a)) x y z --- 1200 +/- 200 --- @ --- -liftCF - :: (Functor f, Fractional a) - => (forall t. f (AD t (Sparse a)) -> AD t (Sparse a)) -- ^ Function on container of values to lift - -> f (CVar s a) -- ^ Container of 'CVar' samples to apply the function to - -> CVar s a -liftCF f cs = CF f cs -{-# INLINE liftCF #-} - --- | Creates a 'CVar' representing a completely independent sample from all --- other 'CVar's containing the exact value given. -constC :: a -> CVar s a -constC = CK -{-# INLINE constC #-} - --- | Lifts a numeric function over the sample represented by a 'CVar'. --- Correctly propagates the uncertainty according to the second-order --- taylor expansion expansion of the function. Note that if the --- higher-degree taylor series terms are large with respect to the mean and --- variance, this approximation may be inaccurate. --- --- Should take any function sufficiently polymorphic over numeric types, so --- you can use things like 'sqrt', 'sin', 'negate', etc. --- --- @ --- ghci> evalCorr $ do --- x <- sampleUncert $ 12.5 +/- 0.8 --- y <- sampleUncert $ 15.9 +/- 0.5 --- resolveUncert $ liftC (\z -> log z ^ 2) (x + y) --- 11.2 +/- 0.2 --- @ --- -liftC - :: Fractional a - => (forall t. AD t (Sparse a) -> AD t (Sparse a)) -- ^ Function on values to lift - -> CVar s a -- ^ 'CVar' sample to apply the function to - -> CVar s a -liftC f = curryH1 $ liftCF (uncurryH1 f) -{-# INLINABLE liftC #-} - --- | Lifts a two-argument (curried) function over the samples represented --- by two 'CVar's. Correctly propagates the uncertainty according to the --- second-order (multivariate) taylor expansion expansion of the function, --- and properly takes into account and keeps track of all --- inter-correlations between the 'CVar' samples. Note that if the --- higher-degree taylor series terms are large with respect to the mean and --- variance, this approximation may be inaccurate. --- --- Should take any function sufficiently polymorphic over numeric types, so --- you can use things like '*', 'atan2', '**', etc. --- --- @ --- ghci> evalCorr $ do --- x <- sampleUncert $ 12.5 +/- 0.8 --- y <- sampleUncert $ 15.9 +/- 0.5 --- resolveUncert $ liftC2 (\a b -> log (a + b) ^ 2) x y --- 11.2 +/- 0.2 --- @ --- -liftC2 - :: Fractional a - => (forall t. AD t (Sparse a) -> AD t (Sparse a) -> AD t (Sparse a)) - -> CVar s a - -> CVar s a - -> CVar s a -liftC2 f = curryH2 $ liftCF (uncurryH2 f) -{-# INLINABLE liftC2 #-} - --- | Lifts a three-argument (curried) function over the samples represented --- by three 'CVar's. See 'liftC2' and 'liftCF' for more details. -liftC3 - :: Fractional a - => (forall t. AD t (Sparse a) -> AD t (Sparse a) -> AD t (Sparse a) -> AD t (Sparse a)) - -> CVar s a - -> CVar s a - -> CVar s a - -> CVar s a -liftC3 f = curryH3 $ liftCF (uncurryH3 f) -{-# INLINABLE liftC3 #-} - --- | Lifts a four-argument (curried) function over the samples represented --- by four 'CVar's. See 'liftC2' and 'liftCF' for more details. -liftC4 - :: Fractional a - => (forall t. AD t (Sparse a) -> AD t (Sparse a) -> AD t (Sparse a) -> AD t (Sparse a) -> AD t (Sparse a)) - -> CVar s a - -> CVar s a - -> CVar s a - -> CVar s a - -> CVar s a -liftC4 f = curryH4 $ liftCF (uncurryH4 f) -{-# INLINABLE liftC4 #-} - --- | Lifts a five-argument (curried) function over the samples represented --- by five 'CVar's. See 'liftC2' and 'liftCF' for more details. -liftC5 - :: Fractional a - => (forall t. AD t (Sparse a) -> AD t (Sparse a) -> AD t (Sparse a) -> AD t (Sparse a) -> AD t (Sparse a) -> AD t (Sparse a)) - -> CVar s a - -> CVar s a - -> CVar s a - -> CVar s a - -> CVar s a - -> CVar s a -liftC5 f = curryH5 $ liftCF (uncurryH5 f) -{-# INLINABLE liftC5 #-} - -instance Fractional a => Num (CVar s a) where - (+) = liftC2 (+) - {-# INLINE (+) #-} - (*) = liftC2 (*) - {-# INLINE (*) #-} - (-) = liftC2 (-) - {-# INLINE (-) #-} - negate = liftC negate - {-# INLINE negate #-} - abs = liftC abs - {-# INLINE abs #-} - signum = liftC signum - {-# INLINE signum #-} - fromInteger = constC . fromInteger - {-# INLINE fromInteger #-} - -instance Fractional a => Fractional (CVar s a) where - recip = liftC recip - {-# INLINE recip #-} - (/) = liftC2 (/) - {-# INLINE (/) #-} - fromRational = constC . fromRational - {-# INLINE fromRational #-} - -instance Floating a => Floating (CVar s a) where - pi = constC pi - {-# INLINE pi #-} - exp = liftC exp - {-# INLINE exp #-} - log = liftC log - {-# INLINE log #-} - sqrt = liftC sqrt - {-# INLINE sqrt #-} - (**) = liftC2 (**) - {-# INLINE (**) #-} - logBase = liftC2 logBase - {-# INLINE logBase #-} - sin = liftC sin - {-# INLINE sin #-} - cos = liftC cos - {-# INLINE cos #-} - asin = liftC asin - {-# INLINE asin #-} - acos = liftC acos - {-# INLINE acos #-} - atan = liftC atan - {-# INLINE atan #-} - sinh = liftC sinh - {-# INLINE sinh #-} - cosh = liftC cosh - {-# INLINE cosh #-} - asinh = liftC asinh - {-# INLINE asinh #-} - acosh = liftC acosh - {-# INLINE acosh #-} - atanh = liftC atanh - {-# INLINE atanh #-} diff --git a/src/Data/Uncertain/MonteCarlo.hs b/src/Data/Uncertain/MonteCarlo.hs deleted file mode 100644 --- a/src/Data/Uncertain/MonteCarlo.hs +++ /dev/null @@ -1,329 +0,0 @@ -{-# LANGUAGE CPP #-} -{-# LANGUAGE ImplicitParams #-} -{-# LANGUAGE ViewPatterns #-} - --- | --- Module : Data.Uncertain.MonteCarlo --- Copyright : (c) Justin Le 2016 --- License : BSD3 --- --- Maintainer : justin@jle.im --- Stability : experimental --- Portability : non-portable --- --- Provides an interface for computing and propagating uncertainty by using --- . --- --- Basically simulates sampling from the distribution represented by the given --- 'Uncert's, applying the function of interest, and aggregating the mean --- and standard deviation of the results. @x '+/-' dx@ is treated as --- a random variable whose probability density is the normal distribution --- with mean @x@ and standard deviation @dx@. --- --- This module attempts to duplicate the API offered by "Data.Uncertain" --- and is meant to be imported qualified alongside "Data.Uncertain" --- --- @ --- import Data.Uncertain --- import qualified Data.Uncertain.MonteCarlo as MC --- @ --- --- Actions are parameterized over all 'PrimMonad' instances, so can be run --- under both 'ST' and 'IO', making it suitable for exploratory purposes. --- All functions require a 'Gen' from "System.Random.MWC" for random value --- generation purposes. --- --- @ --- ghci> import qualified Data.Uncertain.MonteCarlo as MC --- ghci> import System.Random.MWC --- ghci> let x = 1.52 '+/-' 0.07 --- ghci> let y = 781.4 +/- 0.3 --- ghci> let z = 1.53e-1 `'withPrecision'` 3 --- ghci> g <- 'create' --- ghci> cosh x --- 2.4 +/- 0.2 --- ghci> MC.liftU cosh x g --- 2.4 +/- 0.2 --- ghci> exp x / z * sin (y ** z) --- 10.9 +/- 0.9 --- ghci> MC.liftU3 (\a b c -> exp a / c * sin (b**c)) x y z g --- 10.8 +/- 1.0 --- ghci> pi + 3 * logBase x y --- 52 +/- 5 --- ghci> MC.liftU2 (\a b -> pi + 3 * logBase a b) x y g --- 51 +/- 5 --- @ --- - -module Data.Uncertain.MonteCarlo - ( -- * Sampling from an 'Uncert' - sampleUncert - -- * Lifting functions via Monte Carlo simulation - -- ** Fixed iterations - , liftU, liftU2, liftU3, liftU4, liftU5, liftUF - -- ** Variable iterations - , liftU', liftU2', liftU3', liftU4', liftU5', liftUF' - ) - where - -import Control.Monad -import Control.Monad.Primitive -import Data.Hople -import Data.Uncertain (Uncert, fromSamples, uMeanStd) -import System.Random.MWC -import System.Random.MWC.Distributions - -#if __GLASGOW_HASKELL__ < 710 -import Control.Applicative (Applicative) -import Data.Functor ((<$>)) -import Data.Traversable (Traversable(..)) -#endif - --- | Sample a random 'Double' from the distribution specified by an --- @'Uncert' 'Double'@. @x '+/-' dx@ is treated as a random variable whose --- probability density is the normal distribution with mean @x@ and --- standard deviation @dx@. --- -sampleUncert -#if __GLASGOW_HASKELL__ < 710 - :: PrimMonad m -#else - :: (PrimMonad m, Functor m) -#endif - => Uncert Double - -> Gen (PrimState m) - -> m Double -sampleUncert (uMeanStd->(x, dx)) g = normal x dx g -{-# INLINABLE sampleUncert #-} - --- | Lifts a numeric function over an 'Uncert' using a Monte Carlo --- simulation with 1000 samples. --- --- @ --- ghci> g <- 'create' --- ghci> MC.liftU (\x -> log x ^ 2) (12.2 +/- 0.5) g --- 6.3 +/- 0.2 --- @ --- -liftU -#if __GLASGOW_HASKELL__ >= 710 - :: PrimMonad m -#else - :: (PrimMonad m, Functor m) -#endif - => (Double -> Double) - -> Uncert Double - -> Gen (PrimState m) - -> m (Uncert Double) -liftU = liftU' 1000 -{-# INLINE liftU #-} - --- | Lifts a multivariate numeric function on a container (given as an @f --- a -> a@) to work on a container of 'Uncert's using a Monte Carlo --- simulation with 1000 samples. --- --- @ --- ghci> g <- 'create' --- ghci> M.liftUF (\[x,y,z] -> x*y+z) [12.2 +/- 0.5, 56 +/- 2, 0.12 +/- 0.08] g --- 680 +/- 40 --- @ --- -liftUF -#if __GLASGOW_HASKELL__ >= 710 - :: (Traversable f, PrimMonad m) -#else - :: (Traversable f, PrimMonad m, Applicative m) -#endif - => (f Double -> Double) - -> f (Uncert Double) - -> Gen (PrimState m) - -> m (Uncert Double) -liftUF = liftUF' 1000 -{-# INLINE liftUF #-} - --- | Lifts a two-argument (curried) function over two 'Uncert's using --- a Monte Carlo simulation with 1000 samples. --- --- @ --- ghci> g <- 'create' --- ghci> MC.liftU2 (\x y -> x**y) (13.5 +/- 0.1) (1.64 +/- 0.08) --- 70 +/- 20 --- @ --- -liftU2 -#if __GLASGOW_HASKELL__ >= 710 - :: PrimMonad m -#else - :: (PrimMonad m, Applicative m) -#endif - => (Double -> Double -> Double) - -> Uncert Double - -> Uncert Double - -> Gen (PrimState m) - -> m (Uncert Double) -liftU2 = liftU2' 1000 -{-# INLINE liftU2 #-} - --- | Lifts a three-argument (curried) function over three 'Uncert's. See --- 'liftU2' and 'liftUF' for more details. -liftU3 -#if __GLASGOW_HASKELL__ >= 710 - :: PrimMonad m -#else - :: (PrimMonad m, Applicative m) -#endif - => (Double -> Double -> Double -> Double) - -> Uncert Double - -> Uncert Double - -> Uncert Double - -> Gen (PrimState m) - -> m (Uncert Double) -liftU3 = liftU3' 1000 -{-# INLINE liftU3 #-} - --- | Lifts a four-argument (curried) function over four 'Uncert's. See --- 'liftU2' and 'liftUF' for more details. -liftU4 -#if __GLASGOW_HASKELL__ >= 710 - :: PrimMonad m -#else - :: (PrimMonad m, Applicative m) -#endif - => (Double -> Double -> Double -> Double -> Double) - -> Uncert Double - -> Uncert Double - -> Uncert Double - -> Uncert Double - -> Gen (PrimState m) - -> m (Uncert Double) -liftU4 = liftU4' 1000 -{-# INLINE liftU4 #-} - --- | Lifts a five-argument (curried) function over five 'Uncert's. See --- 'liftU2' and 'liftUF' for more details. -liftU5 -#if __GLASGOW_HASKELL__ >= 710 - :: PrimMonad m -#else - :: (PrimMonad m, Applicative m) -#endif - => (Double -> Double -> Double -> Double -> Double -> Double) - -> Uncert Double - -> Uncert Double - -> Uncert Double - -> Uncert Double - -> Uncert Double - -> Gen (PrimState m) - -> m (Uncert Double) -liftU5 = liftU5' 1000 -{-# INLINE liftU5 #-} - --- | Like 'liftU', but allows you to specify the number of samples to run --- the Monte Carlo simulation with. -liftU' -#if __GLASGOW_HASKELL__ >= 710 - :: PrimMonad m -#else - :: (PrimMonad m, Functor m) -#endif - => Int - -> (Double -> Double) - -> Uncert Double - -> Gen (PrimState m) - -> m (Uncert Double) -liftU' n f u g = fromSamples <$> replicateM n samp - where - samp = f <$> sampleUncert u g -{-# INLINABLE liftU' #-} - --- | Like 'liftUF', but allows you to specify the number of samples to run --- the Monte Carlo simulation with. -liftUF' -#if __GLASGOW_HASKELL__ >= 710 - :: (Traversable f, PrimMonad m) -#else - :: (Traversable f, PrimMonad m, Applicative m) -#endif - => Int - -> (f Double -> Double) - -> f (Uncert Double) - -> Gen (PrimState m) - -> m (Uncert Double) -liftUF' n f us g = fromSamples <$> replicateM n samp - where - samp = f <$> traverse (flip sampleUncert g) us -{-# INLINABLE liftUF' #-} - --- | Like 'liftU2', but allows you to specify the number of samples to run --- the Monte Carlo simulation with. -liftU2' -#if __GLASGOW_HASKELL__ >= 710 - :: PrimMonad m -#else - :: (PrimMonad m, Applicative m) -#endif - => Int - -> (Double -> Double -> Double) - -> Uncert Double - -> Uncert Double - -> Gen (PrimState m) - -> m (Uncert Double) -liftU2' n f x y = liftUF' n (uncurryH2 f) (H2 x y) -{-# INLINABLE liftU2' #-} - --- | Like 'liftU3', but allows you to specify the number of samples to run --- the Monte Carlo simulation with. -liftU3' -#if __GLASGOW_HASKELL__ >= 710 - :: PrimMonad m -#else - :: (PrimMonad m, Applicative m) -#endif - => Int - -> (Double -> Double -> Double -> Double) - -> Uncert Double - -> Uncert Double - -> Uncert Double - -> Gen (PrimState m) - -> m (Uncert Double) -liftU3' n f x y z = liftUF' n (uncurryH3 f) (H3 x y z) -{-# INLINABLE liftU3' #-} - --- | Like 'liftU4', but allows you to specify the number of samples to run --- the Monte Carlo simulation with. -liftU4' -#if __GLASGOW_HASKELL__ >= 710 - :: PrimMonad m -#else - :: (PrimMonad m, Applicative m) -#endif - => Int - -> (Double -> Double -> Double -> Double -> Double) - -> Uncert Double - -> Uncert Double - -> Uncert Double - -> Uncert Double - -> Gen (PrimState m) - -> m (Uncert Double) -liftU4' n f x y z a = liftUF' n (uncurryH4 f) (H4 x y z a) -{-# INLINABLE liftU4' #-} - --- | Like 'liftU5', but allows you to specify the number of samples to run --- the Monte Carlo simulation with. -liftU5' -#if __GLASGOW_HASKELL__ >= 710 - :: PrimMonad m -#else - :: (PrimMonad m, Applicative m) -#endif - => Int - -> (Double -> Double -> Double -> Double -> Double -> Double) - -> Uncert Double - -> Uncert Double - -> Uncert Double - -> Uncert Double - -> Uncert Double - -> Gen (PrimState m) - -> m (Uncert Double) -liftU5' n f x y z a b = liftUF' n (uncurryH5 f) (H5 x y z a b) -{-# INLINEABLE liftU5' #-} diff --git a/src/Numeric/Uncertain.hs b/src/Numeric/Uncertain.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Uncertain.hs @@ -0,0 +1,600 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE DeriveGeneric #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE LambdaCase #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE ViewPatterns #-} + +#if __GLASGOW_HASKELL__ >= 708 +{-# LANGUAGE PatternSynonyms #-} +#endif + +-- | +-- Module : Numeric.Uncertain +-- Copyright : (c) Justin Le 2016 +-- License : BSD3 +-- +-- Maintainer : justin@jle.im +-- Stability : experimental +-- Portability : non-portable +-- + +module Numeric.Uncertain + ( -- * 'Uncert' + Uncert +#if __GLASGOW_HASKELL__ >= 708 + , pattern (:+/-) +#endif + -- ** Creating 'Uncert' values + , (+/-), exact, withPrecision, withPrecisionAtBase, withVar, fromSamples + -- ** Inspecting properties + , uMean, uVar, uStd, uMeanVar, uMeanStd, uRange + -- * Applying arbitrary functions + , liftU + , liftU2, liftU3, liftU4, liftU5, liftUF + -- * Utility functions + , uNormalize, uNormalizeAtBase + , uShow, uShowsPrec + ) + where + +import Data.Data +import Data.Foldable (toList, foldl') +import Data.Function +import Data.Hople +import Data.Ord +import GHC.Generics +import Numeric.AD.Mode.Sparse +import Prelude.Compat +import qualified Numeric.AD.Mode.Tower as T + + +-- | Represents an independent experimental value centered around a mean +-- value with "inherent" and independent uncertainty. +-- +-- Mostly useful due to its instances of numeric typeclasses like `Num`, +-- `Fractional`, etc., which allows you to add and multiply and apply +-- arbitrary numerical functions to them and have the uncertainty +-- propagate appropriately. You can also lift arbitrary (sufficiently +-- polymorphic) functions with 'liftU', 'liftUF', 'liftU2' and family. +-- +-- @ +-- ghci> let x = 1.52 '+/-' 0.07 +-- ghci> let y = 781.4 +/- 0.3 +-- ghci> let z = 1.53e-1 `'withPrecision'` 3 +-- ghci> cosh x +-- 2.4 +/- 0.2 +-- ghci> exp x / z * sin (y ** z) +-- 10.9 +/- 0.9 +-- ghci> pi + 3 * logBase x y +-- 52 +/- 5 +-- @ +-- +-- Uncertaintly is properly propagated according to the second-degree +-- taylor series approximations of the applied functions. However, if the +-- higher-degree terms are large with respect to to the means and +-- variances of the uncertain values, these approximations may be +-- inaccurate. +-- +-- Can be created with 'exact' to represent an "exact" measurement with no +-- uncertainty, '+/-' and ':+/-' to specify a standard deviation as +-- a range, 'withPrecision' to specify through decimal precision, and +-- 'withVar' to specify with a variance. Can also be inferred from a list +-- of samples with 'fromSamples' +-- +-- @ +-- 7.13 '+/-' 0.05 +-- 91800 +/- 100 +-- 12.5 `'withVar'` 0.36 +-- 'exact' 7.9512 +-- 81.42 `'withPrecision'` 4 +-- 7 :: Uncertain Double +-- 9.18 :: Uncertain Double +-- 'fromSamples' [12.5, 12.7, 12.6, 12.6, 12.5] +-- @ +-- +-- Can be deconstructed with ':+/-', the pattern synonym/pseudo-constructor +-- which matches on the mean and a standard deviation (supported on GHC +-- 7.8+, with bidirectional constructor functionality supported on GHC +-- 7.10+). You can also access properties with 'uMean', 'uStd', 'uVar', +-- 'uMeanStd', 'uMeanVar', 'uRange', etc. +-- +-- It's important to remember that each "occurrence" represents a unique +-- independent sample, so: +-- +-- @ +-- ghci> let x = 15 '+/-' 2 in x + x +-- 30 +/- 3 +-- +-- ghci> let x = 15 +/- 2 in x*2 +-- 30 +/- 4 +-- @ +-- +-- @x + x@ does not represent adding the same sample to itself twice, it +-- represents /independently/ sampling two values within the range @15 +/- 2@ +-- and adding them together. In general, errors and deviations will cancel +-- each-other out, leading to a smaller uncertainty. +-- +-- However, @x*2@ represents taking /one/ sample and multiplying it by two. +-- This yields a greater uncertainty, because errors and deviations are +-- amplified. +-- +-- Also be aware that the 'Show' instance "normalizes" the result, and +-- won't show any mean/central point to a decimal precision smaller than +-- the uncertainty, rounding off the excess. +-- +data Uncert a = Un { _uMean :: !a + , _uVar :: !a -- ^ maintained to be positive! + } + deriving (Data, Typeable, Generic, Generic1) + +-- | Get the mean/central value/expected value of an 'Uncert'. +uMean :: Uncert a -> a +uMean = _uMean +{-# INLINE uMean #-} + +-- | Get the /variance/ of the uncertainty of an 'Uncert', proportional to +-- the square of "how uncertain" a value is. Is the square of 'uStd'. +uVar :: Uncert a -> a +uVar = _uVar +{-# INLINE uVar #-} + +-- | Get the /standard deviation/ of the uncertainty of an 'Uncert', +-- proportional to "how uncertain" a value is. +-- +-- Very informally, it can be thought of as the interval above and below +-- the mean that about 68% of sampled values will fall under after repeated +-- sampling, or as the range that one is 68% sure the true value is within. +-- +-- Is the square root of 'uVar'. +uStd :: Floating a => Uncert a -> a +uStd = sqrt . uVar +{-# INLINE uStd #-} + +-- | Create an 'Uncert' with an exact value and 0 uncertainty. +exact + :: Num a + => a -- ^ The exact value + -> Uncert a +exact x = Un x 0 +{-# INLINE exact #-} + +infixl 6 +/- +#if __GLASGOW_HASKELL__ >= 708 +infixl 6 :+/- +#endif + +-- | Create an 'Uncert' around a central value and a given "range" of +-- uncertainty. The range is interpreted as the standard deviation of the +-- underlying random variable. Might be preferrable over ':+/-' because it +-- is more general (doesn't require a 'Floating' constraint) and looks +-- a bit nicer. +-- +-- See 'uStd' for more details. +(+/-) + :: Num a + => a -- ^ The mean or central value + -> a -- ^ The standard deviation of the underlying uncertainty + -> Uncert a +x +/- dx = Un x (dx*dx) +{-# INLINE (+/-) #-} + +-- | Create an 'Uncert' around a central value, specifying its uncertainty +-- with a given /variance/. The variance is taken to be proportional to +-- the square of the range of uncertainty. See 'uStd' for more details. +-- +-- "Negative variances" are treated as positive. +withVar + :: Num a + => a -- ^ The mean or central value + -> a -- ^ The variance of the underlying uncertainty + -> Uncert a +withVar x vx = Un x (abs vx) +{-# INLINE withVar #-} + +#if __GLASGOW_HASKELL__ >= 708 +-- | Pattern match on an 'Uncert' with its central value and its standard +-- deviation (see 'uStd' for clarification). +-- +-- Can also be used to /construct/ an 'Uncert', identically as '+/-'. +-- +-- /Note:/ Only supported on GHC 7.8 and above. Bidirectional +-- functionality (to allow use as a constructor) only supported on GHC +-- 7.10 and above. +-- +#if __GLASGOW_HASKELL__ >= 710 +pattern (:+/-) :: () => Floating a => a -> a -> Uncert a +#endif +pattern x :+/- dx <- Un x (sqrt->dx) +#if __GLASGOW_HASKELL__ >= 710 + where + x :+/- dx = Un x (dx*dx) +#endif +#endif + +-- | Infer an 'Uncert' from a given list of independent /samples/ of an +-- underlying uncertain or random distribution. +fromSamples :: Fractional a => [a] -> Uncert a +fromSamples = makeUn . foldStats + where + makeUn (H3 x0 x1 x2) = Un μ v + where + μ = x1/x0 + v = x2/x0 - μ*μ -- maybe use pop var? + foldStats = flip foldl' (H3 0 0 0) $ + \(H3 s0 s1 s2) x -> + H3 (s0 + 1) (s1 + x) (s2 + x*x) +{-# INLINABLE fromSamples #-} + +-- | Retrieve both the mean (central) value and the underlying variance of +-- an 'Uncert' together. +-- +-- @uMeanVar ≡ 'uMean' &&& 'uVar'@ +uMeanVar :: Uncert a -> (a, a) +uMeanVar (Un x vx) = (x, vx) +{-# INLINE uMeanVar #-} + +-- | Retreve both the mean (central) value and the underlying standard +-- deviation of an 'Uncert' together. (See 'uStd' for more details) +-- +-- @uMeanStd ≡ 'uMean' &&& 'uStd'@ +uMeanStd :: Floating a => Uncert a -> (a, a) +uMeanStd (Un x vx) = (x, sqrt vx) +{-# INLINE uMeanStd #-} + +-- | Retrieve the "range" of the underlying distribution of an 'Uncert', +-- derived from the standard deviation, where which approximly 68% of +-- sampled values are expected to occur (or within which you are 68% +-- certain the true value is). +-- +-- @uRange (x +/- dx) ≡ (x - dx, x + dx)@ +uRange :: Floating a => Uncert a -> (a, a) +uRange (uMeanStd->(x, dx)) = (x - dx, x + dx) +{-# INLINABLE uRange #-} + +-- | Like 'withPrecision', except takes a number of "digits" of precision in +-- the desired numeric base. For example, in base 2, takes the number of +-- /bits/ of precision. +-- +-- @'withPrecision' ≡ withPrecisionAtBase 10@ +withPrecisionAtBase + :: (Floating a, RealFrac a) + => Int -- ^ The base to determine precision with respect to + -> a -- ^ The approximate value of the 'Uncert' + -> Int -- ^ The number of "digits" of precision to take + -> Uncert a +withPrecisionAtBase b x p = x' +/- dx' + where + leading :: Int + leading = negate . floor . logBase (fromIntegral b) $ x + uncert :: Int + uncert = leading - 1 + fromIntegral p + rounder = fromIntegral b ** fromIntegral uncert + x' = (/ rounder) . fromIntegral . round' . (* rounder) $ x + dx' = 1 / rounder + round' :: RealFrac a => a -> Integer + round' = round +{-# INLINABLE withPrecisionAtBase #-} + +-- | Create an 'Uncert' about a given approximate central value, with the +-- given number of /digits of precision/ (in decimal notation). +-- +-- @5.21 `withPrecision` 3 ≡ 5.21 '+/-' 0.01@ +withPrecision + :: (Floating a, RealFrac a) + => a -- ^ The approximate value of the 'Uncert' + -> Int -- ^ The number of "digits" of precision to take + -> Uncert a +withPrecision = withPrecisionAtBase 10 +{-# INLINABLE withPrecision #-} + +-- | Like 'uNormalize', but takes a numerical base to round with respect +-- to. +-- +-- @'uNormalize' ≡ uNormalizeAtBase 10@ +uNormalizeAtBase + :: (Floating a, RealFrac a) + => Int -- ^ The base to normalize with respect to + -> Uncert a + -> Uncert a +uNormalizeAtBase b (uMeanStd->(x, dx)) = x' +/- dx' + where + uncert :: Int + uncert = negate . floor . logBase (fromIntegral b) $ dx + rounder = fromIntegral b ** fromIntegral uncert + roundTo = (/ rounder) . fromIntegral . round' . (* rounder) + x' = roundTo x + dx' = roundTo dx + round' :: RealFrac a => a -> Integer + round' = round +{-# INLINABLE uNormalizeAtBase #-} + +-- | Attempts to "normalize" an 'Uncert'. Rounds the uncertainty (the +-- standard deviation) to one digit of precision, and rounds the central +-- moment up to the implied precision. +-- +-- For example, it makes no real sense to have @542.185433 +/- 83.584@, +-- because the extra digits of @542.185433@ past the tens place has no +-- meaning because of the overpowering uncertainty. Normalizing this +-- results in @540 +/- 80@. +-- +-- Note that the 'Show' instance for 'Uncert' normalizes values before +-- showing them. +uNormalize + :: (Floating a, RealFrac a) + => Uncert a + -> Uncert a +uNormalize = uNormalizeAtBase 10 +{-# INLINABLE uNormalize #-} + +instance (Show a, Floating a, RealFrac a) => Show (Uncert a) where + showsPrec d = uShowsPrec d . uNormalize + +-- | Like 'showsPrec' for 'Uncert', but does not normalize the value (see +-- 'uNormalize') before showing. See documentation for 'showsPrec' for +-- more information on how this is meant to be used. +uShowsPrec :: (Show a, Floating a) => Int -> Uncert a -> ShowS +uShowsPrec d (uMeanStd->(x, dx)) = showParen (d > 5) $ + showsPrec 6 x + . showString " +/- " + . showsPrec 6 dx +{-# INLINABLE uShowsPrec #-} + +-- | Like 'show' for 'Uncert', but does not normalize the value (see +-- 'uNormalize') before showing. +-- +-- @'show' ≡ uShow . 'uNormalize'@ +uShow :: (Show a, Floating a) => Uncert a -> String +uShow u = uShowsPrec 0 u "" +{-# INLINABLE uShow #-} + +-- | Lifts a multivariate numeric function on a container (given as an @f +-- a -> a@) to work on a container of 'Uncert's. Correctly propagates the +-- uncertainty according to the second-order (multivariate) taylor +-- expansion of the function. Note that if the higher-degree taylor series +-- terms are large with respect to the means and variances, this +-- approximation may be inaccurate. +-- +-- Should take any function sufficiently polymorphic over numeric types, so +-- you can use things like '*', 'sqrt', 'atan2', etc. +-- +-- @ +-- ghci> liftUF (\[x,y,z] -> x*y+z) [12.2 +/- 0.5, 56 +/- 2, 0.12 +/- 0.08] +-- 680 +/- 40 +-- @ +-- +liftUF + :: (Traversable f, Fractional a) + => (forall s. f (AD s (Sparse a)) -> AD s (Sparse a)) -- ^ Function on container of values to lift + -> f (Uncert a) -- ^ Container of 'Uncert's to apply the function to + -> Uncert a +liftUF f us = Un y vy + where + xs = uMean <$> us + vxs = uVar <$> us + vxsL = toList vxs + (fx, dfxsh) = hessian' f xs + dfxs = fst <$> dfxsh + hess = snd <$> dfxsh + y = fx + hessQuad / 2 + where + hessQuad = dot vxsL + . diag + . toList + $ fmap toList hess + vy = dot vxsL ((^ (2::Int)) <$> dfxs) + dot x = sum . zipWith (*) x . toList + diag = \case [] -> [] + [] :yss -> diag (drop1 <$> yss) + (x:_):yss -> x : diag (drop1 <$> yss) + where + drop1 [] = [] + drop1 (_:zs) = zs +{-# INLINABLE liftUF #-} + +-- | Lifts a numeric function over an 'Uncert'. Correctly propagates the +-- uncertainty according to the second-order taylor expansion expansion of +-- the function. Note that if the higher-degree taylor series terms are +-- large with respect to the mean and variance, this approximation may be +-- inaccurate. +-- +-- Should take any function sufficiently polymorphic over numeric types, so +-- you can use things like 'sqrt', 'sin', 'negate', etc. +-- +-- @ +-- ghci> liftU (\x -> log x ^ 2) (12.2 +/- 0.5) +-- 6.3 +/- 0.2 +-- @ +liftU + :: Fractional a + => (forall s. AD s (T.Tower a) -> AD s (T.Tower a)) -- ^ Function on values to lift + -> Uncert a -- ^ 'Uncert' to apply the function to + -> Uncert a +liftU f (Un x vx) = Un y vy + where + fx:dfx:ddfx:_ = T.diffs0 f x + y = fx + ddfx * vx / 2 + vy = dfx*dfx * vx +{-# INLINABLE liftU #-} + +-- | Lifts a two-argument (curried) function over two 'Uncert's. Correctly +-- propagates the uncertainty according to the second-order (multivariate) +-- taylor expansion expansion of the function. Note that if the +-- higher-degree taylor series terms are large with respect to the mean and +-- variance, this approximation may be inaccurate. +-- +-- Should take any function sufficiently polymorphic over numeric types, so +-- you can use things like '*', 'atan2', '**', etc. +-- +-- @ +-- ghci> liftU2 (\x y -> x**y) (13.5 +/- 0.1) (1.64 +/- 0.08) +-- 70 +/- 10 +-- @ +liftU2 + :: Fractional a + => (forall s. AD s (Sparse a) -> AD s (Sparse a) -> AD s (Sparse a)) + -> Uncert a + -> Uncert a + -> Uncert a +liftU2 f = curryH2 $ liftUF (uncurryH2 f) +{-# INLINABLE liftU2 #-} + +-- | Lifts a three-argument (curried) function over three 'Uncert's. See +-- 'liftU2' and 'liftUF' for more details. +liftU3 + :: Fractional a + => (forall s. AD s (Sparse a) -> AD s (Sparse a) -> AD s (Sparse a) -> AD s (Sparse a)) + -> Uncert a + -> Uncert a + -> Uncert a + -> Uncert a +liftU3 f = curryH3 $ liftUF (uncurryH3 f) +{-# INLINABLE liftU3 #-} + +-- | Lifts a four-argument (curried) function over four 'Uncert's. See +-- 'liftU2' and 'liftUF' for more details. +liftU4 + :: Fractional a + => (forall s. AD s (Sparse a) -> AD s (Sparse a) -> AD s (Sparse a) -> AD s (Sparse a) -> AD s (Sparse a)) + -> Uncert a + -> Uncert a + -> Uncert a + -> Uncert a + -> Uncert a +liftU4 f = curryH4 $ liftUF (uncurryH4 f) +{-# INLINABLE liftU4 #-} + +-- | Lifts a five-argument (curried) function over five 'Uncert's. See +-- 'liftU2' and 'liftUF' for more details. +liftU5 + :: Fractional a + => (forall s. AD s (Sparse a) -> AD s (Sparse a) -> AD s (Sparse a) -> AD s (Sparse a) -> AD s (Sparse a) -> AD s (Sparse a)) + -> Uncert a + -> Uncert a + -> Uncert a + -> Uncert a + -> Uncert a + -> Uncert a +liftU5 f = curryH5 $ liftUF (uncurryH5 f) +{-# INLINABLE liftU5 #-} + +instance Fractional a => Num (Uncert a) where + (+) = liftU2 (+) + {-# INLINE (+) #-} + (*) = liftU2 (*) + {-# INLINE (*) #-} + (-) = liftU2 (-) + {-# INLINE (-) #-} + negate = liftU negate + {-# INLINE negate #-} + abs = liftU abs + {-# INLINE abs #-} + signum = liftU signum + {-# INLINE signum #-} + fromInteger = exact . fromInteger + {-# INLINE fromInteger #-} + +instance Fractional a => Fractional (Uncert a) where + recip = liftU recip + {-# INLINE recip #-} + (/) = liftU2 (/) + {-# INLINE (/) #-} + fromRational = exact . fromRational + {-# INLINE fromRational #-} + +instance Floating a => Floating (Uncert a) where + pi = exact pi + {-# INLINE pi #-} + exp = liftU exp + {-# INLINE exp #-} + log = liftU log + {-# INLINE log #-} + sqrt = liftU sqrt + {-# INLINE sqrt #-} + (**) = liftU2 (**) + {-# INLINE (**) #-} + logBase = liftU2 logBase + {-# INLINE logBase #-} + sin = liftU sin + {-# INLINE sin #-} + cos = liftU cos + {-# INLINE cos #-} + asin = liftU asin + {-# INLINE asin #-} + acos = liftU acos + {-# INLINE acos #-} + atan = liftU atan + {-# INLINE atan #-} + sinh = liftU sinh + {-# INLINE sinh #-} + cosh = liftU cosh + {-# INLINE cosh #-} + asinh = liftU asinh + {-# INLINE asinh #-} + acosh = liftU acosh + {-# INLINE acosh #-} + atanh = liftU atanh + {-# INLINE atanh #-} + +instance Eq a => Eq (Uncert a) where + (==) = (==) `on` uMean + {-# INLINE (==) #-} + (/=) = (/=) `on` uMean + {-# INLINE (/=) #-} + +instance Ord a => Ord (Uncert a) where + compare = comparing uMean + {-# INLINE compare #-} + +instance (Fractional a, Real a) => Real (Uncert a) where + toRational = toRational . uMean + {-# INLINE toRational #-} + +instance RealFrac a => RealFrac (Uncert a) where + properFraction x = (n, d) + where + d = liftU (snd' . properFraction) x + n = fst . properFraction $ uMean x + snd' :: (Int, b) -> b + snd' = snd + {-# INLINABLE properFraction #-} + truncate = truncate . uMean + {-# INLINE truncate #-} + round = round . uMean + {-# INLINE round #-} + ceiling = ceiling . uMean + {-# INLINE ceiling #-} + floor = floor . uMean + {-# INLINE floor #-} + +instance RealFloat a => RealFloat (Uncert a) where + floatRadix = floatRadix . uMean + {-# INLINE floatRadix #-} + floatDigits = floatDigits . uMean + {-# INLINE floatDigits #-} + floatRange = floatRange . uMean + {-# INLINE floatRange #-} + decodeFloat = decodeFloat . uMean + {-# INLINE decodeFloat #-} + exponent = exponent . uMean + {-# INLINE exponent #-} + isNaN = isNaN . uMean + {-# INLINE isNaN #-} + isInfinite = isInfinite . uMean + {-# INLINE isInfinite #-} + isDenormalized = isDenormalized . uMean + {-# INLINE isDenormalized #-} + isNegativeZero = isNegativeZero . uMean + {-# INLINE isNegativeZero #-} + isIEEE = isIEEE . uMean + {-# INLINE isIEEE #-} + encodeFloat a b = exact (encodeFloat a b) + {-# INLINE encodeFloat #-} + significand = liftU significand + {-# INLINE significand #-} + atan2 = liftU2 atan2 + {-# INLINE atan2 #-} + diff --git a/src/Numeric/Uncertain/Correlated.hs b/src/Numeric/Uncertain/Correlated.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Uncertain/Correlated.hs @@ -0,0 +1,90 @@ +{-# LANGUAGE RankNTypes #-} + +-- | +-- Module : Numeric.Uncertain.Correlated +-- Copyright : (c) Justin Le 2016 +-- License : BSD3 +-- +-- Maintainer : justin@jle.im +-- Stability : experimental +-- Portability : non-portable +-- +-- Provides the 'Corr' monad, which allows one to describe complex +-- relationships between random variables and evaluate their propagated +-- uncertainties /respecting/ their inter-correlations. +-- +-- See the "Numeric.Uncertain.Correlated.Interactive" module for an +-- "interactive" and exploratory interface for this module's functionality. +-- + +module Numeric.Uncertain.Correlated + ( -- * 'Corr' + Corr, evalCorr + -- * Uncertain and Correlated Values + , CVar + -- ** Sampling + , sampleUncert, sampleExact, constC + -- ** Resolving + , resolveUncert + -- * Applying arbitrary functions + , liftC, liftC2, liftC3, liftC4, liftC5, liftCF + ) + where + +import Control.Monad.Free +import Control.Monad.Trans.State +import Numeric.Uncertain +import Numeric.Uncertain.Correlated.Internal +import qualified Data.IntMap.Strict as M + +-- | Evaluates the value described by a 'Corr' monad, taking into account +-- inter-correlations between samples. +-- +-- Takes a universally qualified 'Corr', which should not affect usage. +-- See the examples in the documentation for 'Corr'. The univeral +-- qualification is mostly a type system trick to ensure that you aren't +-- allowed to ever use 'evalCorr' to evaluate a 'CVar'. +evalCorr :: Fractional a => (forall s. Corr s a b) -> b +evalCorr c = evalState (corrToState c) (0, M.empty) +{-# INLINABLE evalCorr #-} + +-- | Generate a sample in 'Corr' from an 'Uncert' value, independently from +-- all other samples. +-- +-- Note that you can only sample @'Uncert' a@s within a @'Corr' s a@, meaning +-- that all other "sampled" values are also @a@s. +sampleUncert :: Uncert a -> Corr s a (CVar s a) +sampleUncert u = liftF $ Gen u id +{-# INLINE sampleUncert #-} + +-- | Generate an exact sample in 'Corr' with zero uncertainty, +-- independently from all other samples. +-- +-- Not super useful, since you can do something equivalent with 'constC' +-- or the numeric instances: +-- +-- @ +-- sampleExact x ≡ return ('constC' x) +-- sampleExact 10 ≡ return 10 +-- @ +-- +-- But is provided for completeness alongside 'sampleUncert'. +-- +-- Note that you can exactly sample an @a@ within a @'Corr' s a@, meaning +-- that all other "sampled" values are also @a@s. +-- +sampleExact :: a -> Corr s a (CVar s a) +sampleExact = return . constC +{-# INLINE sampleExact #-} + +-- | "Resolve" an 'Uncert' from a 'CVar' using its potential multiple +-- samples and sample sources, taking into account inter-correlations +-- between 'CVar's and samples. +-- +-- Note that if you use 'sampleUncert' on the result, the new sample will +-- be treated as something completely independent. Usually this should +-- only be used as the "exit point" of a 'Corr' description. +resolveUncert :: CVar s a -> Corr s a (Uncert a) +resolveUncert v = liftF $ Rei v id +{-# INLINE resolveUncert #-} + diff --git a/src/Numeric/Uncertain/Correlated/Interactive.hs b/src/Numeric/Uncertain/Correlated/Interactive.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Uncertain/Correlated/Interactive.hs @@ -0,0 +1,119 @@ +-- | +-- Module : Numeric.Uncertain.Correlated.Interactive +-- Copyright : (c) Justin Le 2016 +-- License : BSD3 +-- +-- Maintainer : justin@jle.im +-- Stability : experimental +-- Portability : non-portable +-- +-- Exports all of the interface of "Numeric.Uncertain.Correlated", except +-- meant to be run in a /ghci/ session "interactively" for exploratory +-- purposes, or in a plain 'IO' action (instead of inside a 'Corr' monad). +-- +-- For example, with the "Numeric.Uncertain.Correlated" interface: +-- +-- @ +-- ghci> evalCorr $ do +-- x <- sampleUncert $ 12.5 +/- 0.8 +-- y <- sampleUncert $ 15.9 +/- 0.5 +-- z <- sampleUncert $ 1.52 +/- 0.07 +-- let k = y**x +-- resolveUncert $ (x+z) * logBase z k +-- 1200 +/- 200 +-- @ +-- +-- And with the interface from this "interactive" module: +-- +-- @ +-- ghci> x <- 'sampleUncert' $ 12.5 +/- 0.8 +-- ghci> y <- sampleUncert $ 15.9 +/- 0.5 +-- ghci> z <- sampleUncert $ 1.52 +/- 0.07 +-- ghci> let k = y**x +-- ghci> 'resolveUncert' $ (x+z) * logBase z k +-- 1200 +/- 200 +-- @ +-- +-- The main purpose of this module is to allow one to use /ghci/ as a fancy +-- "calculator" for computing and exploring propagated uncertainties of +-- complex and potentially correlated samples with uncertainty. +-- +-- Because many of the names overlap with the names from the +-- "Numeric.Uncertain.Correlated" module, it is recommended that you never +-- have both imported at the same time in /ghci/ or in a file, or import +-- them qualified if you must. +-- +-- Also note that all of these methods only work with @'Uncertain' +-- 'Double'@s, and are not polymorphic over different numeric types. +-- +-- Be aware that this module is not robustly tested in heavily concurrent +-- situations/applications. +-- +module Numeric.Uncertain.Correlated.Interactive + ( -- * Uncertain and Correlated Values + CVar, CVarIO + -- ** Sampling + , sampleUncert, sampleExact, constC + -- ** Resolving + , resolveUncert + -- * Applying arbitrary functions + , liftC, liftC2, liftC3, liftC4, liftC5, liftCF + ) + where + +import Control.Monad.ST +import Control.Monad.Trans.State +import Data.IORef +import Data.Tuple +import Numeric.Uncertain +import Numeric.Uncertain.Correlated.Internal +import System.IO.Unsafe (unsafePerformIO) +import qualified Data.IntMap.Strict as M +import qualified Numeric.Uncertain.Correlated as C + +-- | A 'CVar' specialized to work in an "interactive" context, in /ghci/ or +-- 'IO'. +type CVarIO = CVar RealWorld Double + +-- ssh, don't tell anyone we're using 'unsafePerformIO' +globalCorrMap :: IORef (M.Key, M.IntMap (Uncert Double)) +{-# NOINLINE globalCorrMap #-} +globalCorrMap = unsafePerformIO $ newIORef (0, M.empty) + +runCorrIO :: Corr RealWorld Double a -> IO a +runCorrIO c = atomicModifyIORef' globalCorrMap + (swap . runState (corrToState c)) +{-# INLINE runCorrIO #-} + +-- | Generate a sample in 'IO' from an @'Uncert' 'Double'@ value, +-- independently from all other samples. +sampleUncert :: Uncert Double -> IO CVarIO +sampleUncert u = runCorrIO $ C.sampleUncert u +{-# INLINABLE sampleUncert #-} + +-- | Generate an exact sample in 'IO' with zero uncertainty, +-- independently from all other samples. +-- +-- Not super useful, since you can do something equivalent with 'constC' +-- or the numeric instances: +-- +-- @ +-- sampleExact x ≡ return ('constC' x) +-- sampleExact 10 ≡ return 10 +-- @ +-- +-- But is provided for completeness alongside 'sampleUncert'. +sampleExact :: Double -> IO CVarIO +sampleExact d = runCorrIO $ C.sampleExact d +{-# INLINABLE sampleExact #-} + +-- | "Resolve" an 'Uncert' from a 'CVarIO' using its potential multiple +-- samples and sample sources, taking into account inter-correlations +-- between 'CVarIO's and samples. +-- +-- Note that if you use 'sampleUncert' on the result, the new sample will +-- be treated as something completely independent. Usually this should +-- only be used as the "final value" of your computation or exploration. +resolveUncert :: CVarIO -> IO (Uncert Double) +resolveUncert v = runCorrIO $ C.resolveUncert v +{-# INLINABLE resolveUncert #-} diff --git a/src/Numeric/Uncertain/Correlated/Internal.hs b/src/Numeric/Uncertain/Correlated/Internal.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Uncertain/Correlated/Internal.hs @@ -0,0 +1,372 @@ +{-# LANGUAGE DeriveFunctor #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE KindSignatures #-} +{-# LANGUAGE LambdaCase #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# OPTIONS_HADDOCK hide #-} +{-# OPTIONS_HADDOCK prune #-} + +-- | +-- Module : Numeric.Uncertain.Correlated.Internal +-- Copyright : (c) Justin Le 2016 +-- License : BSD3 +-- +-- Maintainer : justin@jle.im +-- Stability : experimental +-- Portability : non-portable +-- +-- Internal utility functions for functionality shared by +-- "Numeric.Uncertain.Correlated" and +-- "Numeric.Uncertain.Correlated.Interactive". +-- + +module Numeric.Uncertain.Correlated.Internal + ( CVar, dephantom + , CorrF(..), Corr + , liftCF + , constC, liftC, liftC2, liftC3, liftC4, liftC5 + , corrToState + ) + where + +import Control.Arrow ((***)) +import Control.Monad.Free +import Control.Monad.Trans.State +import Prelude.Compat +import Data.Hople +import Numeric.Uncertain +import Numeric.AD.Mode.Sparse +import qualified Data.IntMap.Strict as M + + +-- | Represents a single sample (or a value calculated from samples) within +-- the 'Corr' monad. These can be created with 'sampleUncert', +-- 'sampleExact', and 'constC', or made by combinining others with its +-- numeric typeclass instances (like 'Num') or its functions lifting +-- arbitrary numeric functions (like 'liftC2'). These keep track of +-- inter-correlations between sources, and if you add together two 'CVar's +-- that are correlated, their results will reflect this. +-- +-- Can be "resolved" into the uncertain value they represent using +-- 'resolveUncert'. +-- +-- Note that these are parameterized by a dummy phantom parameter 's' so +-- that they can't be "evaluated" out of the 'Corr' they live in with +-- 'evalCorr'. +-- +-- Note that a @'CVar' s a@ can only ever meaningfully "exist" in a @'Corr' +-- s a@, meaning that the all samples within that 'Corr' are of the same +-- type. +data CVar s a where + CK :: a -> CVar s a + CV :: M.Key -> CVar s a + CF :: Functor f + => (forall t. f (AD t (Sparse a)) -> AD t (Sparse a)) + -> f (CVar s a) + -> CVar s a + +-- | Unsafe function to bypass the universal qualification guard for +-- returning 'CVar's from 'Corr's. +dephantom :: CVar s a -> CVar t a +dephantom = \case CK x -> CK x + CV k -> CV k + CF f xs -> CF f (dephantom <$> xs) + +data CorrF :: * -> * -> * -> * where + Gen :: Uncert a -> (CVar s a -> b) -> CorrF s a b + Fun :: Functor f + => (forall t. f (AD t (Sparse a)) -> AD t (Sparse a)) + -> f (CVar s a) + -> (CVar s a -> b) + -> CorrF s a b + Rei :: CVar s a + -> (Uncert a -> b) + -> CorrF s a b + +instance Functor (CorrF s a) where + fmap f = \case Gen u next -> Gen u (f . next) + Fun g us next -> Fun g us (f . next) + Rei v next -> Rei v (f . next) + + +-- | The 'Corr' monad allows us to keep track of correlated and +-- non-independent samples. It fixes a basic "failure" of the 'Uncert' +-- type, which can't describe correlated samples. +-- +-- For example, consider the difference between: +-- +-- @ +-- ghci> sum $ replicate 10 (12.5 '+/-' 0.8) +-- 125 +/- 3 +-- ghci> 10 * (12.5 +/- 0.8) +-- 125 +/- 8 +-- @ +-- +-- The first one represents the addition of ten independent samples, whose +-- errors will in general cancel eachother out. The second one represents +-- sampling once and multiplying it by ten, which will amplify any error by +-- a full factor of 10. +-- +-- See how the 'Corr' monad expresses the above computations: +-- +-- @ +-- ghci> 'evalCorr' $ do +-- x <- 'sampleUncert' $ 12.5 '+/-' 0.8 +-- y1 <- 'resolveUncert' $ sum (replicate 10 x) +-- y2 <- resolveUncert $ 10 * x +-- return (y1, y2) +-- (125 +/- 8, 125 +/- 8) +-- +-- ghci> 'evalCorr' $ do +-- xs <- replicateM 10 ('sampleUncert' (12.5 +/- 0.8)) +-- 'resolveUncert' $ sum xs +-- 125 +/- 3 +-- @ +-- +-- The first example samples once and describes operations on the single +-- sample; the second example samples 10 times with 'replicateM' and sums +-- all of the results. +-- +-- Things are more interesting when you sample multiple variables: +-- +-- @ +-- ghci> 'evalCorr' $ do +-- x <- 'sampleUncert' $ 12.5 '+/-' 0.8 +-- y <- sampleUncert $ 15.9 +/- 0.5 +-- z <- sampleUncert $ 1.52 +/- 0.07 +-- let k = y ** x +-- 'resolveUncert' $ (x+z) * logBase z k +-- 1200 +/- 200 +-- @ +-- +-- The first parameter is a dummy phantom parameter used to prevent 'CVar's +-- from leaking out of the computation (see 'evalCorr'). The second +-- parameter is the numeric type of all samples within the description (for +-- example, if you ever sample an 'Uncert Double', the second parameter wil +-- be 'Double'). The third parameter is the result type of the +-- computation -- the value the 'Corr' is describing. +newtype Corr s a b = Corr { corrFree :: Free (CorrF s a) b + } + deriving (Functor, Applicative, Monad) + +deriving instance MonadFree (CorrF s a) (Corr s a) + +corrToState + :: (Monad m, Fractional a) + => Corr s a b + -> StateT (M.Key, M.IntMap (Uncert a)) m b +corrToState = iterM go . corrFree + where + go = \case + Gen u next -> do + i <- gets fst + modify $ succ *** M.insert i u + next (CV i) + Fun f us next -> + next $ CF f us + Rei v next -> do + u <- gets (getCVar v . snd) + next u + getCVar + :: forall a s. Fractional a + => CVar s a + -> M.IntMap (Uncert a) + -> Uncert a + getCVar cv = liftUF (cVarToF cv) + where + cVarToF + :: CVar s a + -> (forall t. M.IntMap (AD t (Sparse a)) -> AD t (Sparse a)) + cVarToF (CK x) _ = auto x + cVarToF (CV k) us = us M.! k + cVarToF (CF f cs) us = f (flip cVarToF us <$> cs) +{-# INLINABLE corrToState #-} + +-- | Lifts a multivariate numeric function on a container (given as an @f +-- a -> a@) to work on a container of 'CVar's. Correctly propagates the +-- uncertainty according to the second-order (multivariate) taylor +-- expansion of the function, and properly takes into account and keeps +-- track of all inter-correlations between the 'CVar' samples. Note that +-- if the higher-degree taylor series terms are large with respect to the +-- means and variances, this approximation may be inaccurate. +-- +-- Should take any function sufficiently polymorphic over numeric types, so +-- you can use things like '*', 'sqrt', 'atan2', etc. +-- +-- @ +-- ghci> evalCorr $ do +-- x <- sampleUncert $ 12.5 +/- 0.8 +-- y <- sampleUncert $ 15.9 +/- 0.5 +-- z <- sampleUncert $ 1.52 +/- 0.07 +-- resolveUncert $ liftCF (\[a,b,c] -> (a+c) * logBase c (b**a)) x y z +-- 1200 +/- 200 +-- @ +-- +liftCF + :: (Functor f, Fractional a) + => (forall t. f (AD t (Sparse a)) -> AD t (Sparse a)) -- ^ Function on container of values to lift + -> f (CVar s a) -- ^ Container of 'CVar' samples to apply the function to + -> CVar s a +liftCF f cs = CF f cs +{-# INLINE liftCF #-} + +-- | Creates a 'CVar' representing a completely independent sample from all +-- other 'CVar's containing the exact value given. +constC :: a -> CVar s a +constC = CK +{-# INLINE constC #-} + +-- | Lifts a numeric function over the sample represented by a 'CVar'. +-- Correctly propagates the uncertainty according to the second-order +-- taylor expansion expansion of the function. Note that if the +-- higher-degree taylor series terms are large with respect to the mean and +-- variance, this approximation may be inaccurate. +-- +-- Should take any function sufficiently polymorphic over numeric types, so +-- you can use things like 'sqrt', 'sin', 'negate', etc. +-- +-- @ +-- ghci> evalCorr $ do +-- x <- sampleUncert $ 12.5 +/- 0.8 +-- y <- sampleUncert $ 15.9 +/- 0.5 +-- resolveUncert $ liftC (\z -> log z ^ 2) (x + y) +-- 11.2 +/- 0.2 +-- @ +-- +liftC + :: Fractional a + => (forall t. AD t (Sparse a) -> AD t (Sparse a)) -- ^ Function on values to lift + -> CVar s a -- ^ 'CVar' sample to apply the function to + -> CVar s a +liftC f = curryH1 $ liftCF (uncurryH1 f) +{-# INLINABLE liftC #-} + +-- | Lifts a two-argument (curried) function over the samples represented +-- by two 'CVar's. Correctly propagates the uncertainty according to the +-- second-order (multivariate) taylor expansion expansion of the function, +-- and properly takes into account and keeps track of all +-- inter-correlations between the 'CVar' samples. Note that if the +-- higher-degree taylor series terms are large with respect to the mean and +-- variance, this approximation may be inaccurate. +-- +-- Should take any function sufficiently polymorphic over numeric types, so +-- you can use things like '*', 'atan2', '**', etc. +-- +-- @ +-- ghci> evalCorr $ do +-- x <- sampleUncert $ 12.5 +/- 0.8 +-- y <- sampleUncert $ 15.9 +/- 0.5 +-- resolveUncert $ liftC2 (\a b -> log (a + b) ^ 2) x y +-- 11.2 +/- 0.2 +-- @ +-- +liftC2 + :: Fractional a + => (forall t. AD t (Sparse a) -> AD t (Sparse a) -> AD t (Sparse a)) + -> CVar s a + -> CVar s a + -> CVar s a +liftC2 f = curryH2 $ liftCF (uncurryH2 f) +{-# INLINABLE liftC2 #-} + +-- | Lifts a three-argument (curried) function over the samples represented +-- by three 'CVar's. See 'liftC2' and 'liftCF' for more details. +liftC3 + :: Fractional a + => (forall t. AD t (Sparse a) -> AD t (Sparse a) -> AD t (Sparse a) -> AD t (Sparse a)) + -> CVar s a + -> CVar s a + -> CVar s a + -> CVar s a +liftC3 f = curryH3 $ liftCF (uncurryH3 f) +{-# INLINABLE liftC3 #-} + +-- | Lifts a four-argument (curried) function over the samples represented +-- by four 'CVar's. See 'liftC2' and 'liftCF' for more details. +liftC4 + :: Fractional a + => (forall t. AD t (Sparse a) -> AD t (Sparse a) -> AD t (Sparse a) -> AD t (Sparse a) -> AD t (Sparse a)) + -> CVar s a + -> CVar s a + -> CVar s a + -> CVar s a + -> CVar s a +liftC4 f = curryH4 $ liftCF (uncurryH4 f) +{-# INLINABLE liftC4 #-} + +-- | Lifts a five-argument (curried) function over the samples represented +-- by five 'CVar's. See 'liftC2' and 'liftCF' for more details. +liftC5 + :: Fractional a + => (forall t. AD t (Sparse a) -> AD t (Sparse a) -> AD t (Sparse a) -> AD t (Sparse a) -> AD t (Sparse a) -> AD t (Sparse a)) + -> CVar s a + -> CVar s a + -> CVar s a + -> CVar s a + -> CVar s a + -> CVar s a +liftC5 f = curryH5 $ liftCF (uncurryH5 f) +{-# INLINABLE liftC5 #-} + +instance Fractional a => Num (CVar s a) where + (+) = liftC2 (+) + {-# INLINE (+) #-} + (*) = liftC2 (*) + {-# INLINE (*) #-} + (-) = liftC2 (-) + {-# INLINE (-) #-} + negate = liftC negate + {-# INLINE negate #-} + abs = liftC abs + {-# INLINE abs #-} + signum = liftC signum + {-# INLINE signum #-} + fromInteger = constC . fromInteger + {-# INLINE fromInteger #-} + +instance Fractional a => Fractional (CVar s a) where + recip = liftC recip + {-# INLINE recip #-} + (/) = liftC2 (/) + {-# INLINE (/) #-} + fromRational = constC . fromRational + {-# INLINE fromRational #-} + +instance Floating a => Floating (CVar s a) where + pi = constC pi + {-# INLINE pi #-} + exp = liftC exp + {-# INLINE exp #-} + log = liftC log + {-# INLINE log #-} + sqrt = liftC sqrt + {-# INLINE sqrt #-} + (**) = liftC2 (**) + {-# INLINE (**) #-} + logBase = liftC2 logBase + {-# INLINE logBase #-} + sin = liftC sin + {-# INLINE sin #-} + cos = liftC cos + {-# INLINE cos #-} + asin = liftC asin + {-# INLINE asin #-} + acos = liftC acos + {-# INLINE acos #-} + atan = liftC atan + {-# INLINE atan #-} + sinh = liftC sinh + {-# INLINE sinh #-} + cosh = liftC cosh + {-# INLINE cosh #-} + asinh = liftC asinh + {-# INLINE asinh #-} + acosh = liftC acosh + {-# INLINE acosh #-} + atanh = liftC atanh + {-# INLINE atanh #-} diff --git a/src/Numeric/Uncertain/MonteCarlo.hs b/src/Numeric/Uncertain/MonteCarlo.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Uncertain/MonteCarlo.hs @@ -0,0 +1,326 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE ImplicitParams #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE ViewPatterns #-} + +-- | +-- Module : Numeric.Uncertain.MonteCarlo +-- Copyright : (c) Justin Le 2016 +-- License : BSD3 +-- +-- Maintainer : justin@jle.im +-- Stability : experimental +-- Portability : non-portable +-- +-- Provides an interface for computing and propagating uncertainty by using +-- . +-- +-- Basically simulates sampling from the distribution represented by the given +-- 'Uncert's, applying the function of interest, and aggregating the mean +-- and standard deviation of the results. @x '+/-' dx@ is treated as +-- a random variable whose probability density is the normal distribution +-- with mean @x@ and standard deviation @dx@. +-- +-- This module attempts to duplicate the API offered by "Numeric.Uncertain" +-- and is meant to be imported qualified alongside "Numeric.Uncertain" +-- +-- @ +-- import Numeric.Uncertain +-- import qualified Numeric.Uncertain.MonteCarlo as MC +-- @ +-- +-- Actions are parameterized over all 'PrimMonad' instances, so can be run +-- under both 'ST' and 'IO', making it suitable for exploratory purposes. +-- All functions require a 'Gen' from "System.Random.MWC" for random value +-- generation purposes. +-- +-- @ +-- ghci> import qualified Numeric.Uncertain.MonteCarlo as MC +-- ghci> import System.Random.MWC +-- ghci> let x = 1.52 '+/-' 0.07 +-- ghci> let y = 781.4 +/- 0.3 +-- ghci> let z = 1.53e-1 `'withPrecision'` 3 +-- ghci> g <- 'create' +-- ghci> cosh x +-- 2.4 +/- 0.2 +-- ghci> MC.liftU cosh x g +-- 2.4 +/- 0.2 +-- ghci> exp x / z * sin (y ** z) +-- 10.9 +/- 0.9 +-- ghci> MC.liftU3 (\a b c -> exp a / c * sin (b**c)) x y z g +-- 10.8 +/- 1.0 +-- ghci> pi + 3 * logBase x y +-- 52 +/- 5 +-- ghci> MC.liftU2 (\a b -> pi + 3 * logBase a b) x y g +-- 51 +/- 5 +-- @ +-- + +module Numeric.Uncertain.MonteCarlo + ( -- * Sampling from an 'Uncert' + sampleUncert + -- * Lifting functions via Monte Carlo simulation + -- ** Fixed iterations + , liftU, liftU2, liftU3, liftU4, liftU5, liftUF + -- ** Variable iterations + , liftU', liftU2', liftU3', liftU4', liftU5', liftUF' + ) + where + +import Control.Monad +import Control.Monad.Primitive +import Data.Hople +import Numeric.Uncertain (Uncert, fromSamples, uMeanStd) +import Prelude.Compat +import System.Random.MWC +import System.Random.MWC.Distributions + + +-- | Sample a random 'Double' from the distribution specified by an +-- @'Uncert' 'Double'@. @x '+/-' dx@ is treated as a random variable whose +-- probability density is the normal distribution with mean @x@ and +-- standard deviation @dx@. +-- +sampleUncert +#if __GLASGOW_HASKELL__ < 710 + :: PrimMonad m +#else + :: (PrimMonad m, Functor m) +#endif + => Uncert Double + -> Gen (PrimState m) + -> m Double +sampleUncert (uMeanStd->(x, dx)) g = normal x dx g +{-# INLINABLE sampleUncert #-} + +-- | Lifts a numeric function over an 'Uncert' using a Monte Carlo +-- simulation with 1000 samples. +-- +-- @ +-- ghci> g <- 'create' +-- ghci> MC.liftU (\x -> log x ^ 2) (12.2 +/- 0.5) g +-- 6.3 +/- 0.2 +-- @ +-- +liftU +#if __GLASGOW_HASKELL__ >= 710 + :: PrimMonad m +#else + :: (PrimMonad m, Functor m) +#endif + => (Double -> Double) + -> Uncert Double + -> Gen (PrimState m) + -> m (Uncert Double) +liftU = liftU' 1000 +{-# INLINE liftU #-} + +-- | Lifts a multivariate numeric function on a container (given as an @f +-- a -> a@) to work on a container of 'Uncert's using a Monte Carlo +-- simulation with 1000 samples. +-- +-- @ +-- ghci> g <- 'create' +-- ghci> M.liftUF (\[x,y,z] -> x*y+z) [12.2 +/- 0.5, 56 +/- 2, 0.12 +/- 0.08] g +-- 680 +/- 40 +-- @ +-- +liftUF +#if __GLASGOW_HASKELL__ >= 710 + :: (Traversable f, PrimMonad m) +#else + :: (Traversable f, PrimMonad m, Applicative m) +#endif + => (f Double -> Double) + -> f (Uncert Double) + -> Gen (PrimState m) + -> m (Uncert Double) +liftUF = liftUF' 1000 +{-# INLINE liftUF #-} + +-- | Lifts a two-argument (curried) function over two 'Uncert's using +-- a Monte Carlo simulation with 1000 samples. +-- +-- @ +-- ghci> g <- 'create' +-- ghci> MC.liftU2 (\x y -> x**y) (13.5 +/- 0.1) (1.64 +/- 0.08) +-- 70 +/- 20 +-- @ +-- +liftU2 +#if __GLASGOW_HASKELL__ >= 710 + :: PrimMonad m +#else + :: (PrimMonad m, Applicative m) +#endif + => (Double -> Double -> Double) + -> Uncert Double + -> Uncert Double + -> Gen (PrimState m) + -> m (Uncert Double) +liftU2 = liftU2' 1000 +{-# INLINE liftU2 #-} + +-- | Lifts a three-argument (curried) function over three 'Uncert's. See +-- 'liftU2' and 'liftUF' for more details. +liftU3 +#if __GLASGOW_HASKELL__ >= 710 + :: PrimMonad m +#else + :: (PrimMonad m, Applicative m) +#endif + => (Double -> Double -> Double -> Double) + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Gen (PrimState m) + -> m (Uncert Double) +liftU3 = liftU3' 1000 +{-# INLINE liftU3 #-} + +-- | Lifts a four-argument (curried) function over four 'Uncert's. See +-- 'liftU2' and 'liftUF' for more details. +liftU4 +#if __GLASGOW_HASKELL__ >= 710 + :: PrimMonad m +#else + :: (PrimMonad m, Applicative m) +#endif + => (Double -> Double -> Double -> Double -> Double) + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Gen (PrimState m) + -> m (Uncert Double) +liftU4 = liftU4' 1000 +{-# INLINE liftU4 #-} + +-- | Lifts a five-argument (curried) function over five 'Uncert's. See +-- 'liftU2' and 'liftUF' for more details. +liftU5 +#if __GLASGOW_HASKELL__ >= 710 + :: PrimMonad m +#else + :: (PrimMonad m, Applicative m) +#endif + => (Double -> Double -> Double -> Double -> Double -> Double) + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Gen (PrimState m) + -> m (Uncert Double) +liftU5 = liftU5' 1000 +{-# INLINE liftU5 #-} + +-- | Like 'liftU', but allows you to specify the number of samples to run +-- the Monte Carlo simulation with. +liftU' +#if __GLASGOW_HASKELL__ >= 710 + :: PrimMonad m +#else + :: (PrimMonad m, Functor m) +#endif + => Int + -> (Double -> Double) + -> Uncert Double + -> Gen (PrimState m) + -> m (Uncert Double) +liftU' n f u g = fromSamples <$> replicateM n samp + where + samp = f <$> sampleUncert u g +{-# INLINABLE liftU' #-} + +-- | Like 'liftUF', but allows you to specify the number of samples to run +-- the Monte Carlo simulation with. +liftUF' +#if __GLASGOW_HASKELL__ >= 710 + :: (Traversable f, PrimMonad m) +#else + :: (Traversable f, PrimMonad m, Applicative m) +#endif + => Int + -> (f Double -> Double) + -> f (Uncert Double) + -> Gen (PrimState m) + -> m (Uncert Double) +liftUF' n f us g = fromSamples <$> replicateM n samp + where + samp = f <$> traverse (flip sampleUncert g) us +{-# INLINABLE liftUF' #-} + +-- | Like 'liftU2', but allows you to specify the number of samples to run +-- the Monte Carlo simulation with. +liftU2' +#if __GLASGOW_HASKELL__ >= 710 + :: PrimMonad m +#else + :: (PrimMonad m, Applicative m) +#endif + => Int + -> (Double -> Double -> Double) + -> Uncert Double + -> Uncert Double + -> Gen (PrimState m) + -> m (Uncert Double) +liftU2' n f x y = liftUF' n (uncurryH2 f) (H2 x y) +{-# INLINABLE liftU2' #-} + +-- | Like 'liftU3', but allows you to specify the number of samples to run +-- the Monte Carlo simulation with. +liftU3' +#if __GLASGOW_HASKELL__ >= 710 + :: PrimMonad m +#else + :: (PrimMonad m, Applicative m) +#endif + => Int + -> (Double -> Double -> Double -> Double) + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Gen (PrimState m) + -> m (Uncert Double) +liftU3' n f x y z = liftUF' n (uncurryH3 f) (H3 x y z) +{-# INLINABLE liftU3' #-} + +-- | Like 'liftU4', but allows you to specify the number of samples to run +-- the Monte Carlo simulation with. +liftU4' +#if __GLASGOW_HASKELL__ >= 710 + :: PrimMonad m +#else + :: (PrimMonad m, Applicative m) +#endif + => Int + -> (Double -> Double -> Double -> Double -> Double) + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Gen (PrimState m) + -> m (Uncert Double) +liftU4' n f x y z a = liftUF' n (uncurryH4 f) (H4 x y z a) +{-# INLINABLE liftU4' #-} + +-- | Like 'liftU5', but allows you to specify the number of samples to run +-- the Monte Carlo simulation with. +liftU5' +#if __GLASGOW_HASKELL__ >= 710 + :: PrimMonad m +#else + :: (PrimMonad m, Applicative m) +#endif + => Int + -> (Double -> Double -> Double -> Double -> Double -> Double) + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Gen (PrimState m) + -> m (Uncert Double) +liftU5' n f x y z a b = liftUF' n (uncurryH5 f) (H5 x y z a b) +{-# INLINEABLE liftU5' #-} diff --git a/uncertain.cabal b/uncertain.cabal --- a/uncertain.cabal +++ b/uncertain.cabal @@ -2,7 +2,7 @@ -- documentation, see http://haskell.org/cabal/users-guide/ name: uncertain -version: 0.2.0.0 +version: 0.3.0.0 synopsis: Manipulating numbers with inherent experimental/measurement uncertainty description: See . . @@ -25,20 +25,21 @@ location: git://github.com/mstksg/uncertain.git library - exposed-modules: Data.Uncertain - Data.Uncertain.Correlated - Data.Uncertain.Correlated.Interactive - Data.Uncertain.MonteCarlo + exposed-modules: Numeric.Uncertain + Numeric.Uncertain.Correlated + Numeric.Uncertain.Correlated.Interactive + Numeric.Uncertain.MonteCarlo other-modules: Data.Hople - Data.Uncertain.Correlated.Internal + Numeric.Uncertain.Correlated.Internal -- other-extensions: - build-depends: base >=4.6 && < 5 - , ad >= 4 - , containers - , free - , mwc-random - , primitive - , transformers + build-depends: base >= 4.6 && < 5 + , ad >= 4 + , containers >= 0.5 + , free >= 4 + , mwc-random >= 0.10 + , primitive >= 0.1 + , transformers >= 0.2 + , base-compat hs-source-dirs: src default-language: Haskell2010 ghc-options: -Wall