buffon (empty) → 0.1.0.0

raw patch · 4 files changed

+437/−0 lines, 4 filesdep +basedep +monad-primitivedep +mwc-randomsetup-changed

Dependencies added: base, monad-primitive, mwc-random, mwc-random-monad, primitive, transformers

Files

Data/Distribution/Buffon.hs +373/−0
LICENSE +26/−0
Setup.hs +2/−0
buffon.cabal +36/−0

+ Data/Distribution/Buffon.hs view

@@ -0,0 +1,373 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE TypeFamilies #-}++-- From "On Buffon Machines and Numbers" by Flajolet, Pelletier, and Soria++module Data.Distribution.Buffon (+    Buffon, runBuffon, runBuffonWithSystemRandomT,+    toss, third, toNumberWith, toNumberM, toNumber,+    expectationWith, expectationM, expectation,+    bernoulli, if_, mean, evenParity, geometric,+    vonNeumann, polylogarithmic, polylogarithmic',+    poisson, poisson', anotherPoisson, logarithmic, logarithmic',+    alternating, evenAlternating, oddAlternating,+    isAlternating, cosine, sine, cotangent, bump, erf,+    ternary, binary, ternaryBistoch, Bistoch, Bistoch',+    fromBistoch, fromBistoch', bistoch, integrate', arcsin,+    squareRoot, ramanujan, arctan, integrate, createReal,+    pi8, pi4, zeta3+  ) where+import Control.Monad.IO.Class ( MonadIO )+import Control.Monad.Primitive ( PrimMonad )+import Control.Monad.Primitive.Class ( MonadPrim(..) )+import Control.Monad.ST (ST, runST)+import Control.Monad.Trans.Class ( MonadTrans )+import Data.Bits ( testBit )+import Data.Function ( fix )+import Data.Word ( Word64 )+import System.Random.MWC ( Variate )+import System.Random.MWC.Monad ( uniform, Rand, runWithCreate, runWithSystemRandomT )++newtype Buffon m a = Buffon { unBuffon :: Rand m a }+    deriving (Functor, Applicative, Monad, MonadIO, MonadTrans)++instance (PrimMonad m, MonadPrim m) => MonadPrim (Buffon m) where+    type BasePrimMonad (Buffon m) = BasePrimMonad m+    liftPrim = Buffon . liftPrim++runBuffon :: MonadPrim m => Buffon m a -> m a+runBuffon = runWithCreate . unBuffon++runBuffonWithSystemRandomT :: (MonadPrim m, BasePrimMonad m ~ IO) => Buffon m a -> m a+runBuffonWithSystemRandomT = runWithSystemRandomT . unBuffon++----++-- | Toss a coin.  P(toss) = 1/2+toss :: (MonadPrim m) => Buffon m Bool+toss = Buffon uniform++-- | A biased coin. P(third) = 1/3+third :: (MonadPrim m) => Buffon m Bool+third = do+    a <- toss+    b <- toss+    if a && b then third else return (not (a||b))++toNumberWith :: (Fractional n, MonadPrim m) => Int -> Buffon m Bool -> m n+toNumberWith n p = runBuffon (go n 0)+    where go 0 !acc = return (acc/fromIntegral n)+          go c acc = if_ p (go (c-1) (acc+1)) (go (c-1) acc)++toNumberWithSystemRandomT :: (Fractional n) => Int -> Buffon IO Bool -> IO n+toNumberWithSystemRandomT n p = runBuffonWithSystemRandomT (go n 0)+    where go 0 !acc = return (acc/fromIntegral n)+          go c acc = if_ p (go (c-1) (acc+1)) (go (c-1) acc)++toNumberM :: (Fractional n, MonadPrim m) => Buffon m Bool -> m n+toNumberM = toNumberWith 1000000++-- | Estimate the probability of getting True+toNumber :: (Fractional n) => (forall s. Buffon (ST s) Bool) -> n+toNumber p = runST (toNumberM p)++expectationWith :: (Integral a, Fractional n, MonadPrim m) => Int -> Buffon m a -> m n+expectationWith n p = runBuffon (go n 0)+    where go 0 !acc = return (fromIntegral acc/fromIntegral n)+          go c acc = do x <- p; go (c-1) (acc+x)++expectationM :: (Integral a, Fractional n, MonadPrim m) => Buffon m a -> m n+expectationM = expectationWith 1000000++-- | Estimate the expectation+expectation :: (Integral a, Fractional n) => (forall s. Buffon (ST s) a) -> n+expectation p = runST (expectationM p)++----++--  If P(p) = x then p = bernoulli x+-- bernoulli x = fmap (\n -> bits x !! (n+1) == 1) (geometric toss)+--  where bits x returns the bits of the binary expansion of x++-- Cheating bernoulli+-- | P(bernoulli x) = x+bernoulli :: (Ord n, Variate n, Fractional n, MonadPrim m) => n -> Buffon m Bool+bernoulli p = do+    x <- Buffon uniform+    if x <= p then 1 else 0++-- | P(if_ (bernoulli r) (bernoulli p) (bernoulli q)) = rp + (1-r)q+if_ :: (MonadPrim m) => Buffon m Bool -> Buffon m a -> Buffon m a -> Buffon m a+if_ r p q = do b <- r; if b then p else q++-- | P(mean (bernoulli p) (bernoulli q)) = (p+q)/2+mean :: (MonadPrim m) => Buffon m a -> Buffon m a -> Buffon m a+mean p q = if_ toss p q++instance Eq (Buffon m a) where (==) = error "(==) not defined for Buffon"++instance Show (Buffon m a) where show _ = "<buffon>"++-- Warning: These don't obey the laws you'd expect.+-- They are left-biased insofar as:+-- fix (\p -> p*q) == undefined regardless of q (same for +)+-- Practically, this means you want left associated uses of * and +.+instance (MonadPrim m) => Num (Buffon m Bool) where+    p * q = if_ q p 0            -- | P(bernoulli p*bernoulli q) = pq+    p + q = if_ q 1 p            -- | P(bernoulli p+bernoulli q) = p + q - pq+    negate p = not <$> p         -- | P(-bernoulli p)  = 1 - p+    fromInteger 0 = return False -- | P(0)   = 0+    fromInteger 1 = return True  -- | P(1)   = 1+    abs = error "abs not defined for Buffon"+    signum = error "signum not defined for Buffon"++----++-- | P(evenParity (bernoulli p)) = 1/(1+p)+evenParity :: (MonadPrim m) => Buffon m Bool -> Buffon m Bool+evenParity p = fix (\q -> if_ p (if_ p q 0) 1)++-- | von Neumann schema+-- |+-- | Let Perms be a subset of the class of all permutations.+-- | Let F_n be the number of permutations of n items in Perms.+-- | Define F(z) = sum_(n>=0) F_n z^n/n!, i.e. the exponential generating function.+-- | Let U = (U_1, ..., U_n) be a vector of bit streams.+-- | Let type(U) be the permutation that sorts U.  type(U) is the order type of U.+-- | The von Neumann schema is:+-- |+-- | vonNeumann[Perms] p = go 1+-- |    where go k = do+-- |            n <- geometric p+-- |            let U be an n-vector of uniformly distributed bit streams+-- |            if type(U) in Perms then return (n, k) else go (k+1)+-- |+-- | P(fst (vonNeumann[Perms] (bernoulli p)) == n) = F_n/F(p) p^n/n!+-- | s = (1-p)F(p)+-- | P(snd (vonNeumann[Perms] (bernoulli p)) == k) = s(1-s)^(k-1)+-- | sum_(k>=1) k P(snd (vonNeumann[Perms] (bernoulli p)) == k) = 1/s+-- | The actual function takes a function inClass such that P(inClass n) = F_n/n!+vonNeumann :: (MonadPrim m) => (Int -> Buffon m Bool) -> Buffon m Bool -> Buffon m (Int, Int)+vonNeumann inClass p = go 1+    where go !k = do+            n <- geometric p+            if_ (inClass n) (return (n,k)) (go (k+1))++-- | P(geometric (bernoulli p) == r) = (1-p)p^r+-- | Fits the von Neumann schema with F(z) = 1/(1-z) (i.e. Perms = all permutations)+geometric :: (MonadPrim m) => Buffon m Bool -> Buffon m Int+geometric p = go 0+    where go !acc = if_ p (go (acc+1)) (return acc)++-- | P(fst (poisson (bernoulli p)) == r) = exp(-p)p^r/r!+-- | Fits the von Neumann schema with F(z) = exp(z) (i.e. Perms = sorted permutations)+-- Note: P(liftM2 (+) (fst (poisson (bernoulli p))) (fst (poisson (bernoulli q))) == r) = exp(-p-q)(p+q)^r/r!+poisson :: (MonadPrim m) => Buffon m Bool -> Buffon m (Int, Int)+poisson = vonNeumann isSorted+    where isSorted 0 = return True+          isSorted 1 = 1+          isSorted i = loopFalse 0+            where loopFalse j | j < i     = mean (loopTrue j (j+1)) (loopFalse (j+1))+                              | otherwise = isSorted i+                  loopTrue cut j | j < i     = mean (loopTrue cut (j+1)) 0+                                 | otherwise = isSorted cut * isSorted (i-cut)++anotherPoisson :: (MonadPrim m) => Buffon m Bool -> Buffon m (Int, Int)+anotherPoisson = vonNeumann isSorted+    where bernoullis = map (bernoulli . recip) [2 :: Double ..]+          isSorted n = product (take (n-1) bernoullis)++-- | P(poisson' (bernoulli p)) = exp(-p) = P(fst (poisson (bernoulli p)) == 0)+poisson' :: (MonadPrim m) => Buffon m Bool -> Buffon m Bool+poisson' = fmap ((0==) . fst) . poisson++-- | P(fst (logarithmic (bernoulli p)) == r) = -p^r/(rlog(1-p))+-- | Fits the von Neumann schema with F(z) = -log(1-z) (i.e. Perms = cyclic permutations)+logarithmic :: (MonadPrim m) => Buffon m Bool -> Buffon m (Int, Int)+logarithmic = polylogarithmic 1++-- | P(logarithmic' (bernoulli p)) = -p/log(1-p) = P(fst (logarithmic (bernoulli p)) == 1)+logarithmic' :: (MonadPrim m) => Buffon m Bool -> Buffon m Bool+logarithmic' = fmap ((1==) . fst) . logarithmic++-- | P(fst (polylogarthmic k (bernoulli p)) == r) = p^r/(r^k li_k(p))+polylogarithmic :: (MonadPrim m) => Int -> Buffon m Bool -> Buffon m (Int, Int)+polylogarithmic r = vonNeumann (\n -> (isCyclic n)^r)+    where isCyclic 0 = return False+          isCyclic 1 = 1+          isCyclic i = do+              b <- toss+              let inner 0 !acc = isCyclic acc+                  inner j acc = do+                      b' <- toss+                      case compare b b' of+                          LT -> 0+                          EQ -> inner (j-1) (acc+1)+                          GT -> inner (j-1) acc+              inner (i-1) 1++-- | P(polylogarithmic' k (bernoulli p)) = p/li_k(p) = P(fst (polylogarithmic k (bernoulli p)) == 1)+polylogarithmic' :: (MonadPrim m) => Int -> Buffon m Bool -> Buffon m Bool+polylogarithmic' r = fmap ((1==) . fst) . polylogarithmic r++-- | P(fst (alternating (bernoulli p)) == r) = A_r/r! p^r/(sec(p)+tan(p)) = A_r/r! p^r/tan(z/2+pi/4)+-- | Fits the von Neumann schema with F(z) = sec(z) + tan(z) = tan(z/2+pi/4) (i.e. Perms = alternating permutations)+-- | For A_n see OEIS A000111 +alternating :: (MonadPrim m) => Buffon m Bool -> Buffon m (Int, Int)+alternating = vonNeumann (\n -> isAlternating False n [])++-- | P(fst (evenAlternating (bernoulli p)) == r) = A_r/r! p^r/sec(p) for even r, 0 otherwise+-- | Fits the von Neumann schema with F(z) = sec(z) (i.e. Perms = even length alternating permutations)+evenAlternating :: (MonadPrim m) => Buffon m Bool -> Buffon m (Int, Int)+evenAlternating = vonNeumann (\n -> if even n then isAlternating False n [] else 0)++-- | P(fst (oddAlternating (bernoulli p)) == r) = A_r/r! p^r/tan(p) for odd r, 0 otherwise+-- | Fits the von Neumann schema with F(z) = tan(z) (i.e. Perms = odd length alternating permutations)+oddAlternating :: (MonadPrim m) => Buffon m Bool -> Buffon m (Int, Int)+oddAlternating = vonNeumann (\n -> if odd n then isAlternating False n [] else 0)++isAlternating :: (MonadPrim m) => Bool -> Int -> [Bool] -> Buffon m Bool+isAlternating !_  0  _ = 1+isAlternating !_  1  _ = 1+isAlternating !p !i bs = walk bs [] -- is it better to reverse bs here or append in the calls?+    where walk [] acc = walk' acc+          walk (b:bs) acc = do+            b' <- toss+            case if p then compare b b' else compare b' b of+                LT -> 0+                EQ -> walk bs (b:acc)+                GT -> isAlternating (not p) (i-1) (acc++[b'])+          walk' acc = do+            b <- toss+            b' <- toss+            case compare b b' of+                LT -> 0+                EQ -> walk' (b':acc)+                GT -> isAlternating (not p) (i-1) (acc++[p==b'])++-- | P(cosine (bernoulli p)) = cos(p) = P(fst (evenAlternating (bernoulli p)) == 0)+cosine :: (MonadPrim m) => Buffon m Bool -> Buffon m Bool+cosine = fmap ((0==) . fst) . evenAlternating++-- | P(sine (bernoulli p)) = sin(p)+sine :: (MonadPrim m) => Buffon m Bool -> Buffon m Bool+sine p = squareRoot (-cosine p^2)++-- | P(cotangent (bernoulli p)) = pcot(p) = P(fst (oddAlternating (bernoulli p)) == 1)+cotangent :: (MonadPrim m) => Buffon m Bool -> Buffon m Bool+cotangent = fmap ((1==) . fst) . oddAlternating++-- | P(squareRoot (bernoulli p)) = sqrt(p)+squareRoot :: (MonadPrim m) => Buffon m Bool -> Buffon m Bool+squareRoot = bump 1 . negate++-- | P(bump t (bernoulli p)) = (1-p)S(p/2)+-- |    where S(z) = sum_(n>=0) ((t+1)n choose n) z^(tn+n)+bump :: (MonadPrim m) => Int -> Buffon m Bool -> Buffon m Bool+bump t p = go 0+    where go !c = if_ p (bump' (bump' go) c) (return (c==0))+          bump' k !c = mean (k (c+t)) (k (c-1))++-- | Bi(nary )stoch(astic) (grammar)+-- | Let G be a Bistoch X.  It represents the following grammar:+-- | X_i = T (G X_i False) | H (G X_i True)+-- | where X_i is in X, the set of non-terminals and we sequence the list of results from g.+-- | See ternaryBistoch for an example.+type Bistoch k = k -> Bool -> [k]++-- | Bistoch ~ Bistoch'+type Bistoch' k = k -> ([k],[k])+                                       +fromBistoch' :: Bistoch' k -> Bistoch k+fromBistoch' g k False = fst (g k)+fromBistoch' g k True  = snd (g k)++fromBistoch :: Bistoch k -> Bistoch' k+fromBistoch g k = (g k False, g k True)++-- | ternaryBistoch = fromBistoch' $ \() -> ([], [(),(),()])+-- | Represents the grammar: X = T | H X X X+-- | This grammar corresponds to ternary trees.+ternaryBistoch :: Bistoch ()+ternaryBistoch = fromBistoch' $ \() -> ([], [(),(),()])++-- | P(ternary (bernoulli p)) = T(p/2) +-- |    where 2T(p/2) = p(1+T(p/2)^3)+ternary :: (MonadPrim m) => Buffon m Bool -> Buffon m Bool+ternary = bistoch ternaryBistoch ()++-- | P(binary (bernoulli p)) = 1/p - sqrt(1/p^2 - 1)+binary :: (MonadPrim m) => Buffon m Bool -> Buffon m Bool+binary = bistoch (fromBistoch' $ \() -> ([], [(),()])) ()++-- | Sample from a binary stochastic grammar, g, with given starting non-terminal s.+-- | L(g;s) is the language generated by the grammar starting from non-terminal s.+-- | See the Chomsky-Schutzenberger theorem.+-- | (-p)*bistoch g s p = do+-- |    n <- geometric p+-- |    ws <- replicateM n toss+-- |    return (L(g;s) matches ws) +-- | P(bistoch g s (bernoulli p)) = S(p/2)+-- |    where S(z) = sum_(n>=0) S_n z^n+-- |      and S_n are the number of words of length n matched by L(g;s).+bistoch :: (MonadPrim m) => Bistoch k -> k -> Buffon m Bool -> Buffon m Bool+bistoch g s p = matches [s] 1+    where matches     [] c = c+          matches (k:ks) c = (do+            ks' <- g k <$> toss+            matches ks' (matches ks c)) * p++-- | P(ramanujan) = 1/pi+ramanujan :: (MonadPrim m) => Buffon m Bool+ramanujan = do+    x1 <- geometric (toss*toss)+    x2 <- geometric (toss*toss)+    b <- bernoulli (5/9 :: Double)+    let t = 2*(if b then x1 + x2 + 1 else x1 + x2)+    let go 0 !c = return $ c == 0+        go i  c = mean (go (i-1) (c+1)) (go (i-1) (c-1))+    (go t 0 * go t 0) * go t 0++-- Close enough to a real...+-- We could save random bits by flipping on demand and caching the results but it doesn't seem worth it...+createReal :: (MonadPrim m) => Buffon m (Buffon m Bool)+createReal = do+    u <- Buffon uniform+    return (do+        n <- geometric toss+        return (testBit (u :: Word64) n))++-- | P(integrate' f (bernoulli p)) = int_0^p f(w)dw/p+integrate' :: (MonadPrim m) => (Buffon m Bool -> Buffon m Bool) -> Buffon m Bool -> Buffon m Bool+integrate' f p = createReal >>= f . (p*)++-- | P(integrate f (bernoulli p)) = int_0^p f(w)dw+integrate :: (MonadPrim m) => (Buffon m Bool -> Buffon m Bool) -> Buffon m Bool -> Buffon m Bool+integrate f p = integrate' f p * p++-- | P(arctan (bernoulli p)) = atan(p)+arctan :: (MonadPrim m) => Buffon m Bool -> Buffon m Bool+arctan p = (do u <- createReal; evenParity (p*(p*(u*u)))) * p++-- | P(arcsin (bernoulli p)) = asin(p)/2+arcsin :: (MonadPrim m) => Buffon m Bool -> Buffon m Bool+arcsin p = mean (integrate (\w -> evenParity w * squareRoot (-w*w)) p) (-squareRoot (-p*p))++-- | P(erf (bernoulli p)) = sqrt(pi)erf(p)/2+erf :: (MonadPrim m) => Buffon m Bool -> Buffon m Bool+erf p = integrate (\w -> poisson' (w*w)) p++-- | P(pi8) = pi/8 (via (atan(1/2) + atan(1/3))/2)+pi8 :: (MonadPrim m) => Buffon m Bool+pi8 = mean (arctan toss) (arctan third)++-- | P(pi4) = pi/4 (via atan(1))+pi4 :: (MonadPrim m) => Buffon m Bool+pi4 = do u <- createReal; evenParity (u*u)++-- | P(zeta3) = 3zeta(3)/4+zeta3 :: (MonadPrim m) => Buffon m Bool+zeta3 = integrate' (\x -> integrate' (\y -> integrate' (\z -> evenParity (x*y*z)) 1) 1) 1

+ LICENSE view

@@ -0,0 +1,26 @@+Copyright (c) 2015, Derek Elkins
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+1. Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+2. 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.
+
+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.

+ Setup.hs view

@@ -0,0 +1,2 @@+import Distribution.Simple
+main = defaultMain

+ buffon.cabal view

@@ -0,0 +1,36 @@+-- Initial buffon.cabal generated by cabal init.  For further +-- documentation, see http://haskell.org/cabal/users-guide/++name:                buffon+version:             0.1.0.0+synopsis:            An implementation of Buffon machines.+description:         An implementation of everything in "On Buffon Machines and Numbers".+license:             BSD2+license-file:        LICENSE+author:              Derek Elkins+maintainer:          derek.a.elkins+github@gmail.com+homepage:            https://github.com/derekelkins/buffon+bug-reports:         https://github.com/derekelkins/buffon/issues+copyright:           2015 Derek Elkins+category:            Math+build-type:          Simple+-- extra-source-files:  +cabal-version:       >=1.10++source-repository head+  type: git+  location: git://github.com/derekelkins/buffon.git++library+  exposed-modules: Data.Distribution.Buffon+  -- other-modules:       +  -- other-extensions:    +  build-depends:       +    base >=4.8 && <4.9,+    mwc-random >= 0.13 && <0.14,+    mwc-random-monad >=0.7 && <0.8,+    monad-primitive >= 0.1 && <0.2,+    transformers >= 0.4 && <0.5,+    primitive >= 0.6 && <0.7+  -- hs-source-dirs:      +  default-language:    Haskell2010