diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,13 @@ +Changelog +========= + +Version 0.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 + diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,148 @@ +Uncertain +========= + +[![Build Status](https://travis-ci.org/mstksg/uncertain.svg?branch=master)](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 . + +## 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 diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/src/Data/Hople.hs b/src/Data/Hople.hs new file mode 100644 --- /dev/null +++ b/src/Data/Hople.hs @@ -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 #-} + diff --git a/src/Data/Uncertain.hs b/src/Data/Uncertain.hs new file mode 100644 --- /dev/null +++ b/src/Data/Uncertain.hs @@ -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 #-} + diff --git a/src/Data/Uncertain/Correlated.hs b/src/Data/Uncertain/Correlated.hs new file mode 100644 --- /dev/null +++ b/src/Data/Uncertain/Correlated.hs @@ -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 #-} + diff --git a/src/Data/Uncertain/Correlated/Interactive.hs b/src/Data/Uncertain/Correlated/Interactive.hs new file mode 100644 --- /dev/null +++ b/src/Data/Uncertain/Correlated/Interactive.hs @@ -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 #-} diff --git a/src/Data/Uncertain/Correlated/Internal.hs b/src/Data/Uncertain/Correlated/Internal.hs new file mode 100644 --- /dev/null +++ b/src/Data/Uncertain/Correlated/Internal.hs @@ -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 #-} diff --git a/src/Data/Uncertain/MonteCarlo.hs b/src/Data/Uncertain/MonteCarlo.hs new file mode 100644 --- /dev/null +++ b/src/Data/Uncertain/MonteCarlo.hs @@ -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 +-- . +-- +-- Basically simulates sampling from the distribution represented by the given +-- 'Uncert's, applying the function of interest, and aggregating the mean +-- and standard deviation of the results. @x '+/-' dx@ is treated as +-- a random variable whose probability density is the normal distribution +-- with mean @x@ and standard deviation @dx@. +-- +-- This module attempts to duplicate the API offered by "Data.Uncertain" +-- and is meant to be imported qualified alongside "Data.Uncertain" +-- +-- @ +-- import Data.Uncertain +-- import qualified Data.Uncertain.MonteCarlo as MC +-- @ +-- +-- Actions are parameterized over all 'PrimMonad' instances, so can be run +-- under both 'ST' and 'IO', making it suitable for exploratory purposes. +-- All functions require a 'Gen' from "System.Random.MWC" for random value +-- generation purposes. +-- +-- @ +-- ghci> import qualified Data.Uncertain.MonteCarlo as MC +-- ghci> import System.Random.MWC +-- ghci> let x = 1.52 '+/-' 0.07 +-- ghci> let y = 781.4 +/- 0.3 +-- ghci> let z = 1.53e-1 `'withPrecision'` 3 +-- ghci> g <- 'create' +-- ghci> cosh x +-- 2.4 +/- 0.2 +-- ghci> MC.liftU cosh x g +-- 2.4 +/- 0.2 +-- ghci> exp x / z * sin (y ** z) +-- 10.9 +/- 0.9 +-- ghci> MC.liftU3 (\a b c -> exp a / c * sin (b**c)) x y z g +-- 10.8 +/- 1.0 +-- ghci> pi + 3 * logBase x y +-- 52 +/- 5 +-- ghci> MC.liftU2 (\a b -> pi + 3 * logBase a b) x y g +-- 51 +/- 5 +-- @ +-- + +module Data.Uncertain.MonteCarlo + ( -- * Sampling from an 'Uncert' + sampleUncert + -- * Lifting functions via Monte Carlo simulation + -- ** Fixed iterations + , liftU, liftU2, liftU3, liftU4, liftU5, liftUF + -- ** Variable iterations + , liftU', liftU2', liftU3', liftU4', liftU5', liftUF' + ) + where + +import Control.Monad +import Control.Monad.Primitive +import Data.Hople +import Data.Uncertain (Uncert, fromSamples, uMeanStd) +import System.Random.MWC +import System.Random.MWC.Distributions + +#if __GLASGOW_HASKELL__ < 710 +import Control.Applicative (Applicative) +import Data.Functor ((<$>)) +import Data.Traversable (Traversable(..)) +#endif + +-- | Sample a random 'Double' from the distribution specified by an +-- @'Uncert' 'Double'@. @x '+/-' dx@ is treated as a random variable whose +-- probability density is the normal distribution with mean @x@ and +-- standard deviation @dx@. +-- +sampleUncert +#if __GLASGOW_HASKELL__ < 710 + :: PrimMonad m +#else + :: (PrimMonad m, Functor m) +#endif + => Uncert Double + -> Gen (PrimState m) + -> m Double +sampleUncert (uMeanStd->(x, dx)) g = normal x dx g +{-# INLINABLE sampleUncert #-} + +-- | Lifts a numeric function over an 'Uncert' using a Monte Carlo +-- simulation with 1000 samples. +-- +-- @ +-- ghci> g <- 'create' +-- ghci> MC.liftU (\x -> log x ^ 2) (12.2 +/- 0.5) g +-- 6.3 +/- 0.2 +-- @ +-- +liftU +#if __GLASGOW_HASKELL__ >= 710 + :: PrimMonad m +#else + :: (PrimMonad m, Functor m) +#endif + => (Double -> Double) + -> Uncert Double + -> Gen (PrimState m) + -> m (Uncert Double) +liftU = liftU' 1000 +{-# INLINE liftU #-} + +-- | Lifts a multivariate numeric function on a container (given as an @f +-- a -> a@) to work on a container of 'Uncert's using a Monte Carlo +-- simulation with 1000 samples. +-- +-- @ +-- ghci> g <- 'create' +-- ghci> M.liftUF (\[x,y,z] -> x*y+z) [12.2 +/- 0.5, 56 +/- 2, 0.12 +/- 0.08] g +-- 680 +/- 40 +-- @ +-- +liftUF +#if __GLASGOW_HASKELL__ >= 710 + :: (Traversable f, PrimMonad m) +#else + :: (Traversable f, PrimMonad m, Applicative m) +#endif + => (f Double -> Double) + -> f (Uncert Double) + -> Gen (PrimState m) + -> m (Uncert Double) +liftUF = liftUF' 1000 +{-# INLINE liftUF #-} + +-- | Lifts a two-argument (curried) function over two 'Uncert's using +-- a Monte Carlo simulation with 1000 samples. +-- +-- @ +-- ghci> g <- 'create' +-- ghci> MC.liftU2 (\x y -> x**y) (13.5 +/- 0.1) (1.64 +/- 0.08) +-- 70 +/- 20 +-- @ +-- +liftU2 +#if __GLASGOW_HASKELL__ >= 710 + :: PrimMonad m +#else + :: (PrimMonad m, Applicative m) +#endif + => (Double -> Double -> Double) + -> Uncert Double + -> Uncert Double + -> Gen (PrimState m) + -> m (Uncert Double) +liftU2 = liftU2' 1000 +{-# INLINE liftU2 #-} + +-- | Lifts a three-argument (curried) function over three 'Uncert's. See +-- 'liftU2' and 'liftUF' for more details. +liftU3 +#if __GLASGOW_HASKELL__ >= 710 + :: PrimMonad m +#else + :: (PrimMonad m, Applicative m) +#endif + => (Double -> Double -> Double -> Double) + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Gen (PrimState m) + -> m (Uncert Double) +liftU3 = liftU3' 1000 +{-# INLINE liftU3 #-} + +-- | Lifts a four-argument (curried) function over four 'Uncert's. See +-- 'liftU2' and 'liftUF' for more details. +liftU4 +#if __GLASGOW_HASKELL__ >= 710 + :: PrimMonad m +#else + :: (PrimMonad m, Applicative m) +#endif + => (Double -> Double -> Double -> Double -> Double) + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Gen (PrimState m) + -> m (Uncert Double) +liftU4 = liftU4' 1000 +{-# INLINE liftU4 #-} + +-- | Lifts a five-argument (curried) function over five 'Uncert's. See +-- 'liftU2' and 'liftUF' for more details. +liftU5 +#if __GLASGOW_HASKELL__ >= 710 + :: PrimMonad m +#else + :: (PrimMonad m, Applicative m) +#endif + => (Double -> Double -> Double -> Double -> Double -> Double) + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Gen (PrimState m) + -> m (Uncert Double) +liftU5 = liftU5' 1000 +{-# INLINE liftU5 #-} + +-- | Like 'liftU', but allows you to specify the number of samples to run +-- the Monte Carlo simulation with. +liftU' +#if __GLASGOW_HASKELL__ >= 710 + :: PrimMonad m +#else + :: (PrimMonad m, Functor m) +#endif + => Int + -> (Double -> Double) + -> Uncert Double + -> Gen (PrimState m) + -> m (Uncert Double) +liftU' n f u g = fromSamples <$> replicateM n samp + where + samp = f <$> sampleUncert u g +{-# INLINABLE liftU' #-} + +-- | Like 'liftUF', but allows you to specify the number of samples to run +-- the Monte Carlo simulation with. +liftUF' +#if __GLASGOW_HASKELL__ >= 710 + :: (Traversable f, PrimMonad m) +#else + :: (Traversable f, PrimMonad m, Applicative m) +#endif + => Int + -> (f Double -> Double) + -> f (Uncert Double) + -> Gen (PrimState m) + -> m (Uncert Double) +liftUF' n f us g = fromSamples <$> replicateM n samp + where + samp = f <$> traverse (flip sampleUncert g) us +{-# INLINABLE liftUF' #-} + +-- | Like 'liftU2', but allows you to specify the number of samples to run +-- the Monte Carlo simulation with. +liftU2' +#if __GLASGOW_HASKELL__ >= 710 + :: PrimMonad m +#else + :: (PrimMonad m, Applicative m) +#endif + => Int + -> (Double -> Double -> Double) + -> Uncert Double + -> Uncert Double + -> Gen (PrimState m) + -> m (Uncert Double) +liftU2' n f x y = liftUF' n (uncurryH2 f) (H2 x y) +{-# INLINABLE liftU2' #-} + +-- | Like 'liftU3', but allows you to specify the number of samples to run +-- the Monte Carlo simulation with. +liftU3' +#if __GLASGOW_HASKELL__ >= 710 + :: PrimMonad m +#else + :: (PrimMonad m, Applicative m) +#endif + => Int + -> (Double -> Double -> Double -> Double) + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Gen (PrimState m) + -> m (Uncert Double) +liftU3' n f x y z = liftUF' n (uncurryH3 f) (H3 x y z) +{-# INLINABLE liftU3' #-} + +-- | Like 'liftU4', but allows you to specify the number of samples to run +-- the Monte Carlo simulation with. +liftU4' +#if __GLASGOW_HASKELL__ >= 710 + :: PrimMonad m +#else + :: (PrimMonad m, Applicative m) +#endif + => Int + -> (Double -> Double -> Double -> Double -> Double) + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Gen (PrimState m) + -> m (Uncert Double) +liftU4' n f x y z a = liftUF' n (uncurryH4 f) (H4 x y z a) +{-# INLINABLE liftU4' #-} + +-- | Like 'liftU5', but allows you to specify the number of samples to run +-- the Monte Carlo simulation with. +liftU5' +#if __GLASGOW_HASKELL__ >= 710 + :: PrimMonad m +#else + :: (PrimMonad m, Applicative m) +#endif + => Int + -> (Double -> Double -> Double -> Double -> Double -> Double) + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Uncert Double + -> Gen (PrimState m) + -> m (Uncert Double) +liftU5' n f x y z a b = liftUF' n (uncurryH5 f) (H5 x y z a b) +{-# INLINEABLE liftU5' #-} diff --git a/uncertain.cabal b/uncertain.cabal new file mode 100644 --- /dev/null +++ b/uncertain.cabal @@ -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 . + . + Documentation maintained at +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