From 32cdc463254b7d8e35ba71da2d6238e1fcc6e70e Mon Sep 17 00:00:00 2001 From: Michael Pilosov Date: Sat, 18 May 2024 01:47:18 +0000 Subject: [PATCH] working code --- .gitignore | 2 ++ README.md | 12 +++++++ Setup.hs | 2 ++ app/Main.hs | 29 ++++++++++++++++ mud.cabal | 73 +++++++++++++++++++++++++++++++++++++++++ package.yaml | 64 ++++++++++++++++++++++++++++++++++++ src/AcceptReject.hs | 33 +++++++++++++++++++ src/MUD.hs | 51 ++++++++++++++++++++++++++++ src/RandomUtils.hs | 25 ++++++++++++++ stack.yaml | 67 +++++++++++++++++++++++++++++++++++++ stack.yaml.lock | 13 ++++++++ test/Main.hs | 12 +++++++ test/TestKDE.hs | 26 +++++++++++++++ test/TestMUD.hs | 17 ++++++++++ test/TestRandomUtils.hs | 18 ++++++++++ 15 files changed, 444 insertions(+) create mode 100644 .gitignore create mode 100644 README.md create mode 100644 Setup.hs create mode 100644 app/Main.hs create mode 100644 mud.cabal create mode 100644 package.yaml create mode 100644 src/AcceptReject.hs create mode 100644 src/MUD.hs create mode 100644 src/RandomUtils.hs create mode 100644 stack.yaml create mode 100644 stack.yaml.lock create mode 100644 test/Main.hs create mode 100644 test/TestKDE.hs create mode 100644 test/TestMUD.hs create mode 100644 test/TestRandomUtils.hs diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..c368d45 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +.stack-work/ +*~ \ No newline at end of file diff --git a/README.md b/README.md new file mode 100644 index 0000000..d147468 --- /dev/null +++ b/README.md @@ -0,0 +1,12 @@ +# mud + + +Run the demo: +```sh +stack run +``` + +## tests +``` +stack test +``` diff --git a/Setup.hs b/Setup.hs new file mode 100644 index 0000000..9a994af --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/app/Main.hs b/app/Main.hs new file mode 100644 index 0000000..e877a49 --- /dev/null +++ b/app/Main.hs @@ -0,0 +1,29 @@ +-- app/Main.hs +import Text.Printf (printf) +import MUD (updated, getPF, wme) +import RandomUtils (generateUniformPoints, generateNormalPoints) +import Data.List (elemIndex) +import Data.Maybe (fromJust) + +main :: IO () +main = do + let f = id + let numData = 100 + let mesh = [-1, -0.995 .. 1] + lam_lg <- generateUniformPoints 1000 (-1.0, 1.0) + let trueValue = 0.1 + let noiseLevel = 0.121 + noise <- generateNormalPoints numData 0 noiseLevel + let d = replicate numData trueValue + let noisyData = zipWith (+) d noise + let g x = wme (f x) noisyData noiseLevel + let q = getPF g lam_lg + let v = updated g q mesh + + -- print $ v + + let maxIndex = fromJust $ elemIndex (maximum v) v + let mudPoint = mesh !! maxIndex + + putStrLn $ printf "For (true = %f with sigma = %f and N = %d), the MUD point from mesh is: %f" trueValue noiseLevel numData mudPoint + diff --git a/mud.cabal b/mud.cabal new file mode 100644 index 0000000..cdb5401 --- /dev/null +++ b/mud.cabal @@ -0,0 +1,73 @@ +cabal-version: 2.2 + +-- This file has been generated from package.yaml by hpack version 0.36.0. +-- +-- see: https://github.com/sol/hpack + +name: mud +version: 0.1.0.0 +description: Please see the README on GitHub at +homepage: https://github.com/mathematicalmichael/mud#readme +bug-reports: https://github.com/mathematicalmichael/mud/issues +author: Michael Pilosov +maintainer: consistentbayes@gmail.com +copyright: 2024 Michael Pilosov +license: BSD-3-Clause +build-type: Simple +extra-source-files: + README.md + +source-repository head + type: git + location: https://github.com/mathematicalmichael/mud + +library + exposed-modules: + MUD + RandomUtils + other-modules: + AcceptReject + Paths_mud + autogen-modules: + Paths_mud + hs-source-dirs: + src + ghc-options: -Wall -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wmissing-export-lists -Wmissing-home-modules -Wpartial-fields -Wredundant-constraints + build-depends: + base >=4.7 && <5 + , random + default-language: Haskell2010 + +executable mud-exe + main-is: Main.hs + other-modules: + Paths_mud + autogen-modules: + Paths_mud + hs-source-dirs: + app + ghc-options: -Wall -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wmissing-export-lists -Wmissing-home-modules -Wpartial-fields -Wredundant-constraints -threaded -rtsopts -with-rtsopts=-N + build-depends: + base >=4.7 && <5 + , mud + , random + default-language: Haskell2010 + +test-suite test-suite + type: exitcode-stdio-1.0 + main-is: Main.hs + other-modules: + TestKDE + TestMUD + TestRandomUtils + Paths_mud + autogen-modules: + Paths_mud + hs-source-dirs: + test + ghc-options: -Wall -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wmissing-export-lists -Wmissing-home-modules -Wpartial-fields -Wredundant-constraints -threaded -rtsopts -with-rtsopts=-N + build-depends: + base >=4.7 && <5 + , mud + , random + default-language: Haskell2010 diff --git a/package.yaml b/package.yaml new file mode 100644 index 0000000..3d82966 --- /dev/null +++ b/package.yaml @@ -0,0 +1,64 @@ +name: mud +version: 0.1.0.0 +github: "mathematicalmichael/mud" +license: BSD-3-Clause +author: "Michael Pilosov" +maintainer: "consistentbayes@gmail.com" +copyright: "2024 Michael Pilosov" + +extra-source-files: +- README.md + +# Metadata used when publishing your package +# synopsis: Short description of your package +# category: Uncategorized + +# To avoid duplicated efforts in documentation and dealing with the +# complications of embedding Haddock markup inside cabal files, it is +# common to point users to the README.md file. +description: Please see the README on GitHub at + +dependencies: +- base >= 4.7 && < 5 +- random + +ghc-options: +- -Wall +- -Wcompat +- -Widentities +- -Wincomplete-record-updates +- -Wincomplete-uni-patterns +- -Wmissing-export-lists +- -Wmissing-home-modules +- -Wpartial-fields +- -Wredundant-constraints + +library: + source-dirs: src + exposed-modules: + - MUD + - RandomUtils + +executables: + mud-exe: + main: Main.hs + source-dirs: app + ghc-options: + - -threaded + - -rtsopts + - -with-rtsopts=-N + dependencies: + - mud + +tests: + test-suite: + main: Main.hs + source-dirs: test + ghc-options: + - -threaded + - -rtsopts + - -with-rtsopts=-N + dependencies: + - mud + - random + diff --git a/src/AcceptReject.hs b/src/AcceptReject.hs new file mode 100644 index 0000000..7527739 --- /dev/null +++ b/src/AcceptReject.hs @@ -0,0 +1,33 @@ +module AcceptReject (main) where + +import System.Random +import Control.Monad (replicateM) + +-- Acceptance-Rejection sampling for normal distribution +generateNormalPointsAR :: Int -> Double -> Double -> IO [Double] +generateNormalPointsAR n mu sigma = do + gen <- newStdGen + let samples = take n $ filter acceptCandidate $ randomPairs gen + return $ map ((\x -> mu + sigma * x) . fst) samples + where + acceptCandidate (u1, u2) = + let x = sqrt (-2 * log u1) * cos (2 * pi * u2) + in u2 < exp (-0.5 * x * x) + randomPairs g = zip (randoms g) (randoms g) + +-- Generate random points from a uniform distribution in the range [a, b] +generateUniformPoints :: Int -> Double -> Double -> IO [Double] +generateUniformPoints n a b = replicateM n (randomRIO (a, b)) + +main :: IO () +main = do + -- Generate 10 random points from a uniform distribution in the range [-1, 1] + uniformPoints <- generateUniformPoints 10 (-1) 1 + putStrLn "Uniform distribution points:" + print uniformPoints + + -- Generate 10 random points from a normal distribution using acceptance-rejection sampling + normalPointsAR <- generateNormalPointsAR 10 0 1 + putStrLn "Normal distribution points (Acceptance-Rejection):" + print normalPointsAR + diff --git a/src/MUD.hs b/src/MUD.hs new file mode 100644 index 0000000..3da8fd6 --- /dev/null +++ b/src/MUD.hs @@ -0,0 +1,51 @@ +module MUD (uniformPDF, gaussianPDF, kde, wme, getPF, updated) where + +-- Function to evaluate Uniform PDF at a point x with range [a, b] +uniformPDF :: Double -> Double -> Double -> Double +uniformPDF a b x + | x >= a && x <= b = 1 / (b - a) + | otherwise = 0 + +-- Function to evaluate Gaussian PDF at a point x with mean mu and standard deviation sigma +gaussianPDF :: Double -> Double -> Double -> Double +gaussianPDF x mu sigma = (1 / sqrt (2 * pi * sigma**2)) * exp (-(x - mu)**2 / (2 * sigma**2)) + +-- Kernel Density Estimation function +kde :: [Double] -> Double -> (Double -> Double) +kde ctlpoints bandwidth = \x -> (1 / (n * actualBandwidth)) * sum [gaussianPDF x xi actualBandwidth | xi <- ctlpoints] + where + n = fromIntegral (length ctlpoints) + actualBandwidth + | bandwidth == 0 = 1.06 * stdDev ctlpoints * n ** (-1/5) + | otherwise = bandwidth + +-- Calculate standard deviation +stdDev :: [Double] -> Double +stdDev xs = sqrt $ sum [(x - mean) ** 2 | x <- xs] / n + where + n = fromIntegral (length xs) + mean = sum xs / n + +-- Estimate the push-forward distribution +getPF :: (Double -> Double) -> [Double] -> (Double -> Double) +getPF f x = kde y 0 + where + y = map f x + +-- Weighted Mean Error Functional +wme :: Double -> [Double] -> Double -> Double +wme y d s = (1 / sqrt n) * sum [ (y - di) / std | di <- d ] + where + n = fromIntegral (length d) + std + | s == 0 = stdDev d + | otherwise = s + +-- Updated density computation. +updated :: (Double -> Double) -> (Double -> Double) -> [Double] -> [Double] +updated f d x = [ i * (o / p) | (i, o, p) <- zip3 a b c] + where + y = map f x + a = map (\xi -> gaussianPDF xi 0 1) x + b = map (\yi -> gaussianPDF yi 0 1) y + c = map d y diff --git a/src/RandomUtils.hs b/src/RandomUtils.hs new file mode 100644 index 0000000..0ab5924 --- /dev/null +++ b/src/RandomUtils.hs @@ -0,0 +1,25 @@ +-- src/RandomUtils.hs +module RandomUtils (generateUniformPoints, generateNormalPoints) where + +import System.Random +import Control.Monad (replicateM) + +-- Generate random points from a uniform distribution in the range [a, b] +generateUniformPoints :: Int -> (Double, Double) -> IO [Double] +generateUniformPoints n (a, b) = replicateM n (randomRIO (a, b)) + +-- Generate random points from a normal distribution using the Box-Muller transform +generateNormalPoints :: Int -> Double -> Double -> IO [Double] +generateNormalPoints n mu sigma = do + gen <- newStdGen + let randomFloats = take (2 * n) $ randomRs (0, 1) gen + return $ boxMullerTransform mu sigma randomFloats + +-- Box-Muller transform to generate normally distributed random values +boxMullerTransform :: Double -> Double -> [Double] -> [Double] +boxMullerTransform mu sigma (u1:u2:us) = + let z0 = sqrt (-2 * log u1) * cos (2 * pi * u2) + z1 = sqrt (-2 * log u1) * sin (2 * pi * u2) + in (mu + z0 * sigma) : (mu + z1 * sigma) : boxMullerTransform mu sigma us +boxMullerTransform _ _ _ = [] -- If input list is not in pairs, return empty list + diff --git a/stack.yaml b/stack.yaml new file mode 100644 index 0000000..654544d --- /dev/null +++ b/stack.yaml @@ -0,0 +1,67 @@ +# This file was automatically generated by 'stack init' +# +# Some commonly used options have been documented as comments in this file. +# For advanced use and comprehensive documentation of the format, please see: +# https://docs.haskellstack.org/en/stable/yaml_configuration/ + +# Resolver to choose a 'specific' stackage snapshot or a compiler version. +# A snapshot resolver dictates the compiler version and the set of packages +# to be used for project dependencies. For example: +# +# resolver: lts-22.21 +# resolver: nightly-2024-05-06 +# resolver: ghc-9.6.5 +# +# The location of a snapshot can be provided as a file or url. Stack assumes +# a snapshot provided as a file might change, whereas a url resource does not. +# +# resolver: ./custom-snapshot.yaml +# resolver: https://example.com/snapshots/2023-01-01.yaml +resolver: + url: https://raw.githubusercontent.com/commercialhaskell/stackage-snapshots/master/lts/22/22.yaml + +# User packages to be built. +# Various formats can be used as shown in the example below. +# +# packages: +# - some-directory +# - https://example.com/foo/bar/baz-0.0.2.tar.gz +# subdirs: +# - auto-update +# - wai +packages: +- . +# Dependency packages to be pulled from upstream that are not in the resolver. +# These entries can reference officially published versions as well as +# forks / in-progress versions pinned to a git hash. For example: +# +# extra-deps: +# - acme-missiles-0.3 +# - git: https://github.com/commercialhaskell/stack.git +# commit: e7b331f14bcffb8367cd58fbfc8b40ec7642100a +# +# extra-deps: [] + +# Override default flag values for local packages and extra-deps +# flags: {} + +# Extra package databases containing global packages +# extra-package-dbs: [] + +# Control whether we use the GHC we find on the path +# system-ghc: true +# +# Require a specific version of Stack, using version ranges +# require-stack-version: -any # Default +# require-stack-version: ">=2.15" +# +# Override the architecture used by Stack, especially useful on Windows +# arch: i386 +# arch: x86_64 +# +# Extra directories used by Stack for building +# extra-include-dirs: [/path/to/dir] +# extra-lib-dirs: [/path/to/dir] +# +# Allow a newer minor version of GHC than the snapshot specifies +# compiler-check: newer-minor diff --git a/stack.yaml.lock b/stack.yaml.lock new file mode 100644 index 0000000..e690081 --- /dev/null +++ b/stack.yaml.lock @@ -0,0 +1,13 @@ +# This file was autogenerated by Stack. +# You should not edit this file by hand. +# For more information, please see the documentation at: +# https://docs.haskellstack.org/en/stable/lock_files + +packages: [] +snapshots: +- completed: + sha256: 4be1ca5d31689b524a7f0f17a439bbe9136465213edc498e9a395899a670f2aa + size: 718486 + url: https://raw.githubusercontent.com/commercialhaskell/stackage-snapshots/master/lts/22/22.yaml + original: + url: https://raw.githubusercontent.com/commercialhaskell/stackage-snapshots/master/lts/22/22.yaml diff --git a/test/Main.hs b/test/Main.hs new file mode 100644 index 0000000..605807d --- /dev/null +++ b/test/Main.hs @@ -0,0 +1,12 @@ +module Main (main) where + +import TestKDE (tests) +import TestMUD (tests) +import TestRandomUtils (tests) + +main :: IO () +main = do + TestKDE.tests + TestRandomUtils.tests + TestMUD.tests + diff --git a/test/TestKDE.hs b/test/TestKDE.hs new file mode 100644 index 0000000..33368db --- /dev/null +++ b/test/TestKDE.hs @@ -0,0 +1,26 @@ +-- TestKDE.hs +module TestKDE (tests) where + +import MUD (kde) +import RandomUtils (generateUniformPoints) +import Data.List (intercalate) + +tests :: IO () +tests = do + -- List of 21 equispaced points from -5 to 5 + let grdpoints :: [Double] + grdpoints = [-5, -4.5 .. 5] + + -- Generate 10 random points from a normal distribution + controlPoints <- generateUniformPoints 10 (-1.0, 1.0) + + -- Fit the KDE to the random points with a specified bandwidth + let bandwidth = 0.5 + let kdeEstimator = kde controlPoints bandwidth + + -- Evaluate the KDE on the grid from -5 to 5 + let evaluations = map kdeEstimator grdpoints + + putStrLn "KDE evaluations at points from -5 to 5:" + putStrLn $ intercalate "\n" $ zipWith (\x y -> "x = " ++ show x ++ ", KDE = " ++ show y) grdpoints evaluations + diff --git a/test/TestMUD.hs b/test/TestMUD.hs new file mode 100644 index 0000000..05ecf8f --- /dev/null +++ b/test/TestMUD.hs @@ -0,0 +1,17 @@ +-- TestMUD.hs +module TestMUD (tests) where + +import MUD (updated, getPF) +import RandomUtils (generateUniformPoints) + +tests :: IO () +tests = do + -- Example of Updated density + let grdpoints = [-5, -4.5 .. 5] + -- let f = (\x -> x) + let f = id + lam_sm <- generateUniformPoints 100 (-1.0, 1.0) + let p = getPF f lam_sm + let u = updated f p grdpoints + print u + diff --git a/test/TestRandomUtils.hs b/test/TestRandomUtils.hs new file mode 100644 index 0000000..c480057 --- /dev/null +++ b/test/TestRandomUtils.hs @@ -0,0 +1,18 @@ +-- TestRandomUtils.hs +module TestRandomUtils (tests) where + +import RandomUtils (generateUniformPoints, generateNormalPoints) + +-- Define your tests here +tests :: IO () +tests = do + -- Generate 10 random points from a uniform distribution in the range [-1, 1] + uniformPoints <- generateUniformPoints 10 (-1, 1) + putStrLn "Uniform distribution points:" + print uniformPoints + + -- Generate 10 random points from a normal distribution using the Box-Muller transform + normalPoints <- generateNormalPoints 10 0 1 + putStrLn "Normal distribution points (Box-Muller):" + print normalPoints +