uncertain (empty) → 0.2.0.0
raw patch · 11 files changed
+1842/−0 lines, 11 filesdep +addep +basedep +containerssetup-changed
Dependencies added: ad, base, containers, free, mwc-random, primitive, transformers
Files
- CHANGELOG.md +13/−0
- LICENSE +30/−0
- README.md +148/−0
- Setup.hs +2/−0
- src/Data/Hople.hs +90/−0
- src/Data/Uncertain.hs +603/−0
- src/Data/Uncertain/Correlated.hs +90/−0
- src/Data/Uncertain/Correlated/Interactive.hs +119/−0
- src/Data/Uncertain/Correlated/Internal.hs +374/−0
- src/Data/Uncertain/MonteCarlo.hs +329/−0
- uncertain.cabal +44/−0
+ CHANGELOG.md view
@@ -0,0 +1,13 @@+Changelog+=========++Version 0.2.0.0+---------------++<https://github.com/mstksg/uncertain/releases/tag/v0.2.0.0>++* 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+
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2014, Justin Le++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Justin Le nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,148 @@+Uncertain+=========++[](https://travis-ci.org/mstksg/uncertain)++Provides tools to manipulate numbers with inherent experimental/measurement+uncertainty, and propagates them through functions based on principles from+statistics.++Documentation maintained at <https://mstksg.github.io/uncertain>.++## Usage++```haskell+import Data.Uncertain+```++### Create numbers++```haskell+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 descontructed/analyzed with `:+/-` (pattern synonym/pseudo-constructor+matching on the mean and standard deviation), `uMean`, `uStd`, `uVar`, etc.++### Manipulate with error propagation++```haskell+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+```++Propagates uncertainty using second-order multivariate Taylor expansions of+functions, computed using the *[ad][]* library.++[ad]: https://hackage.haskell.org/package/ad++#### Arbitrary numeric functions++```haskell+ghci> liftUF (\[x,y,z] -> x*y+z)+ [ 12.2 +/- 0.5+ , 56 +/- 2+ , 0.12 +/- 0.08+ ]+680 +/- 40+```++## Correlated samples++Can propagate uncertainty on complex functions take from potentially correlated+samples.++```haskell+ghci> import Data.Uncertain.Correlated+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+```++### "Interactive" Exploratory Mode++*Correlated* module functionality can be used in *ghci* or `IO` or `ST`, for+"interactive" exploration.++```haskell+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+```++## Monte Carlo-based propagation of uncertainty++Provides a module for propagating uncertainty using [Monte Carlo+simulations][]++[Monte Carlo simulations]: https://en.wikipedia.org/wiki/Monte_Carlo_method++```haskell+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+```++## Comparisons++Note that this is very different from other libraries with similar data types+(like from [intervals][] and [rounding][]); these do not attempt to maintain intervals or+simply digit precisions; they instead are intended to model actual+experimental and measurement data with their uncertainties, and apply+functions to the data with the uncertainties and properly propagating the+errors with sound statistical principles.++[data-interval]: https://hackage.haskell.org/package/intervals+[rounding]: https://hackage.haskell.org/package/rounding++For a clear example, take++```haskell+> (52 +/- 6) + (39 +/- 4)+91. +/- 7.+```++In a library like [interval], this would result in `91 +/- 10` (that is, a+lower bound of 46 + 35 and an upper bound of 58 + 43). However, with+experimental data, errors in two independent samples tend to "cancel out", and+result in an overall aggregate uncertainty in the sum of approximately 7.++## Copyright++Copyright (c) Justin Le 2016
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ src/Data/Hople.hs view
@@ -0,0 +1,90 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE DeriveFoldable #-}+{-# LANGUAGE DeriveFunctor #-}+{-# LANGUAGE DeriveTraversable #-}+{-# OPTIONS_HADDOCK hide #-}+{-# OPTIONS_HADDOCK prune #-}++-- |+-- Module : Data.Hople+-- Copyright : (c) Justin Le 2016+-- License : BSD3+--+-- Maintainer : justin@jle.im+-- Stability : experimental+-- Portability : non-portable+--+-- Homogeneous strict tuples used for implementing 'liftU2', etc.++module Data.Hople+ ( H1(..)+ , H2(..)+ , H3(..)+ , H4(..)+ , H5(..)+ , curryH1, curryH2, curryH3, curryH4, curryH5+ , uncurryH1, uncurryH2, uncurryH3, uncurryH4, uncurryH5+ )+ where++#if __GLASGOW_HASKELL__ < 710+import Data.Foldable+import Data.Traversable+#endif+++data H1 a = H1 !a+ deriving (Functor, Foldable, Traversable, Show)++data H2 a = H2 !a !a+ deriving (Functor, Foldable, Traversable, Show)++data H3 a = H3 !a !a !a+ deriving (Functor, Foldable, Traversable, Show)++data H4 a = H4 !a !a !a !a+ deriving (Functor, Foldable, Traversable, Show)++data H5 a = H5 !a !a !a !a !a+ deriving (Functor, Foldable, Traversable, Show)++curryH1 :: (H1 a -> a) -> a -> a+curryH1 f x = f (H1 x)+{-# INLINE curryH1 #-}++curryH2 :: (H2 a -> a) -> a -> a -> a+curryH2 f x y = f (H2 x y)+{-# INLINE curryH2 #-}++curryH3 :: (H3 a -> a) -> a -> a -> a -> a+curryH3 f x y z = f (H3 x y z)+{-# INLINE curryH3 #-}++curryH4 :: (H4 a -> a) -> a -> a -> a -> a -> a+curryH4 f x y z a = f (H4 x y z a)+{-# INLINE curryH4 #-}++curryH5 :: (H5 a -> a) -> a -> a -> a -> a -> a -> a+curryH5 f x y z a b = f (H5 x y z a b)+{-# INLINE curryH5 #-}++uncurryH1 :: (a -> a) -> H1 a -> a+uncurryH1 f (H1 x) = f x+{-# INLINE uncurryH1 #-}++uncurryH2 :: (a -> a -> a) -> H2 a -> a+uncurryH2 f (H2 x y) = f x y+{-# INLINE uncurryH2 #-}++uncurryH3 :: (a -> a -> a -> a) -> H3 a -> a+uncurryH3 f (H3 x y z) = f x y z+{-# INLINE uncurryH3 #-}++uncurryH4 :: (a -> a -> a -> a -> a) -> H4 a -> a+uncurryH4 f (H4 x y z a) = f x y z a+{-# INLINE uncurryH4 #-}++uncurryH5 :: (a -> a -> a -> a -> a -> a) -> H5 a -> a+uncurryH5 f (H5 x y z a b) = f x y z a b+{-# INLINE uncurryH5 #-}+
+ src/Data/Uncertain.hs view
@@ -0,0 +1,603 @@+{-# 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 #-}+
@@ -0,0 +1,90 @@+{-# 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 #-}+
@@ -0,0 +1,119 @@+-- |+-- 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 #-}
@@ -0,0 +1,374 @@+{-# 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 view
@@ -0,0 +1,329 @@+{-# 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' #-}
+ uncertain.cabal view
@@ -0,0 +1,44 @@+-- Initial uncertain.cabal generated by cabal init. For further +-- documentation, see http://haskell.org/cabal/users-guide/++name: uncertain+version: 0.2.0.0+synopsis: Manipulating numbers with inherent experimental/measurement uncertainty+description: See <https://github.com/mstksg/uncertain/blob/master/README.md README.md>.+ .+ Documentation maintained at <https://mstksg.github.io/uncertain>+homepage: https://github.com/mstksg/uncertain+bug-reports: https://github.com/mstksg/uncertain/issues+license: BSD3+license-file: LICENSE+author: Justin Le+maintainer: justin@jle.im+copyright: (c) Justin Le 2016+category: Math+build-type: Simple+extra-source-files: README.md+ CHANGELOG.md+cabal-version: >=1.10++source-repository head+ type: git+ location: git://github.com/mstksg/uncertain.git++library+ exposed-modules: Data.Uncertain+ Data.Uncertain.Correlated+ Data.Uncertain.Correlated.Interactive+ Data.Uncertain.MonteCarlo+ other-modules: Data.Hople+ Data.Uncertain.Correlated.Internal+ -- other-extensions: + build-depends: base >=4.6 && < 5+ , ad >= 4+ , containers+ , free+ , mwc-random+ , primitive+ , transformers+ hs-source-dirs: src+ default-language: Haskell2010+ ghc-options: -Wall