working code

.gitignore

@@ -0,0 +1,2 @@
+.stack-work/
+*~

README.md

@@ -0,0 +1,12 @@
+# mud
+
+
+Run the demo:
+```sh
+stack run
+```
+
+## tests
+```
+stack test
+```

Setup.hs

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

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
+

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 <https://github.com/mathematicalmichael/mud#readme>
+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

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 <https://github.com/mathematicalmichael/mud#readme>
+
+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
+

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
+

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

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
+

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

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

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
+

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
+

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
+

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
+