MCMC function

Author: mrtkp9993

Project.toml

@@ -1,7 +1,7 @@
 name = "NumericalAlgorithms"
 uuid = "3b47cccd-2e0f-43d2-aaa1-0b61cbee3591"
 authors = ["Murat Koptur <[email protected]>"]
-version = "0.1.5"
+version = "0.1.6"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/MCMC.jl

@@ -0,0 +1,20 @@
+include("Random.jl")
+
+function MCMC(f::Function, count::UInt64, burnInPer::UInt64)
+    x = zeros(count + burnInPer)
+    x[1] = 1.0
+
+    for i = 2:(count + burnInPer)
+        currx = x[i - 1]
+        # TODO: Replace randn and rand with my own implementation
+        propx = currx + randn(1)[1]
+        a = f(propx) / f(currx)
+        if (rand() < a)
+            x[i] = propx
+        else
+            x[i] = currx
+        end
+    end
+
+    return x[(burnInPer+1):end]
+end

src/NumericalAlgorithms.jl

@@ -12,6 +12,7 @@ module NumericalAlgorithms
 include("Constants.jl")
 include("Differentiation.jl")
 include("Integration.jl")
+include("MCMC.jl")
 include("Random.jl")
 include("RootFinding.jl")
 include("StatisticalTests.jl")
@@ -25,5 +26,6 @@ export Dual
 export CalcSingleIntegral, CalcDoubleIntegral, CalcMonteCarloIntegral
 export LFG, RANMAR, vanderCorputSeq, haltonSeq, faureSeq, sobolSeq
 export RunsTest
+export MCMC
 
 end # module

test/runtests.jl

@@ -155,4 +155,12 @@ end
     z2 = RunsTest(rndnums_rng2)
     z2 = round(z2, digits=2)
     @test -2.58 < z1 < 2.58    
+end
+
+@testset "Markov Chain Monte Carlo" begin
+    f(x::Float64)::Float64 = x < 0 ? 0 : exp(-x)
+    a::UInt64 = 10000
+    b::UInt64 = 2500
+    res = MCMC(f, a, b)
+    @test length(res) == a
 end