uncertain 0.2.0.0 → 0.3.0.0
raw patch · 14 files changed
+1537/−1538 lines, 14 filesdep +base-compat
Dependencies added: base-compat
Files
- CHANGELOG.md +10/−1
- README.md +4/−4
- src/Data/Hople.hs +2/−5
- src/Data/Uncertain.hs +0/−603
- src/Data/Uncertain/Correlated.hs +0/−90
- src/Data/Uncertain/Correlated/Interactive.hs +0/−119
- src/Data/Uncertain/Correlated/Internal.hs +0/−374
- src/Data/Uncertain/MonteCarlo.hs +0/−329
- src/Numeric/Uncertain.hs +600/−0
- src/Numeric/Uncertain/Correlated.hs +90/−0
- src/Numeric/Uncertain/Correlated/Interactive.hs +119/−0
- src/Numeric/Uncertain/Correlated/Internal.hs +372/−0
- src/Numeric/Uncertain/MonteCarlo.hs +326/−0
- uncertain.cabal +14/−13
CHANGELOG.md view
@@ -1,6 +1,15 @@ Changelog ========= +Version 0.3.0.0+---------------++<https://github.com/mstksg/uncertain/releases/tag/v0.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
README.md view
@@ -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
src/Data/Hople.hs view
@@ -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
− src/Data/Uncertain.hs
@@ -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 #-}-
@@ -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 #-}-
@@ -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 #-}
@@ -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 #-}
− src/Data/Uncertain/MonteCarlo.hs
@@ -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--- <https://en.wikipedia.org/wiki/Monte_Carlo_method Monte Carlo simulations>.------ 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' #-}
+ src/Numeric/Uncertain.hs view
@@ -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 #-}+
@@ -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 #-}+
@@ -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 #-}
@@ -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 #-}
+ src/Numeric/Uncertain/MonteCarlo.hs view
@@ -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+-- <https://en.wikipedia.org/wiki/Monte_Carlo_method Monte Carlo simulations>.+--+-- 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' #-}
uncertain.cabal view
@@ -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 <https://github.com/mstksg/uncertain/blob/master/README.md README.md>. .@@ -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