monte carlo integration

Author: mrtkp9993

Project.toml

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

README.md

@@ -15,28 +15,28 @@ Install the package with ```add https://github.com/mrtkp9993/NumericalAlgorithms
 Currently implemented:
 
 * Root finding algorithms
-    * Secant method
-    * Broyden's method
+  * Secant method
+  * Broyden's method
 * Differentation
-    * Automatic differentiation via dual numbers
+  * Automatic differentiation via dual numbers
 * Integration
-   * Composite Simpson - One dim.
-   * Double Simpson - Two dim.
-   * [WIP] Monte Carlo Integration
+  * Composite Simpson - One dim.
+  * Double Simpson - Two dim.
+  * Monte Carlo Integration - Arbitrary dimension
 * Random Number Generators (RNGs)
-    * Pseudo-random numbers
-        * Lagged Fibonacci generator
-        * RANMAR
-    * [WIP] Quasi-random numbers
-        * van der Corput sequences
-        * Halton sequences
-        * Faure sequences
-        * [WIP] Sobol sequences
-
+  * Pseudo-random numbers
+    * Lagged Fibonacci generator
+    * RANMAR
+  * Quasi-random numbers
+    * van der Corput sequences
+    * Halton sequences
+    * Faure sequences
+    * Sobol sequences
 * Statistical Tests
-    * Wald–Wolfowitz runs test
+  * Wald–Wolfowitz runs test
+
+## Next steps to do
 
-*Todo*
 * Optimization
 * Differential Equations (ODEs)& Systems of diff. eqs.
 * Partial Differential Equations (PDEs)

src/Integration.jl

@@ -79,4 +79,22 @@ function CalcDoubleIntegral(f::Function, a::Real, b::Real, c::Function, d::Funct
     t = t * hx / 3.0
     println("Approximation to integral: $t")
     return t
+end
+
+function CalcMonteCarloIntegral(f::Function, an, bn, n::UInt64)
+    @assert length(an) == length(bn)
+    rndNums = zeros(n, length(an))
+    mult = 1
+    for i = 1:length(an)
+        ani = an[i]
+        bni = bn[i]
+        rndNums[:, i] = uniformSeq(ani, bni, n)
+        mult *= bni - ani
+    end
+    mult /= n
+    intt = 0
+    for j = 1:n
+        intt += f(rndNums[j,:]...)
+    end
+    mult * intt
 end

src/NumericalAlgorithms.jl

@@ -22,8 +22,8 @@ export C_SQRT1_2, C_SQRTPI, C_SQRT2PI
 export primes_1000
 export FindRoot1D, FindRootND
 export Dual
-export CalcSingleIntegral, CalcDoubleIntegral
-export LFG, RANMAR, vanderCorputSeq, haltonSeq, faureSeq
+export CalcSingleIntegral, CalcDoubleIntegral, CalcMonteCarloIntegral
+export LFG, RANMAR, vanderCorputSeq, haltonSeq, faureSeq, sobolSeq
 export RunsTest
 
 end # module

src/Random.jl

@@ -19,7 +19,7 @@ end
 "Constructor for LFG with default parameter values"
 LFG() = LFG(24, 55, 2^32, [rand(UInt) % 2^32 for i = 1:55], +)
 "Constructor for LFG with user supplied values, p, q, m"
-LFG(p, q, m) = LPG(p, q, m, [rand(UInt) % m for i = 1:q], +)
+LFG(p, q, m) = LFG(p, q, m, [rand(UInt) % m for i = 1:q], +)
 
 
 function Base.rand(rng::LFG)
@@ -111,4 +111,17 @@ function faureSeq(nb, m, b)
         end
     end
     bn
+end
+
+"Sobol sequence"
+function sobolSeq()
+    throw("not implemented yet")
+end
+
+"Generate n random numbers with uniform dist. in [a, b]"
+function uniformSeq(a, b, n)
+    rng1 = LFG(9739, 23209, 2^32 - 1)
+    rndNums = [rand(rng1) for i = 1:n]
+    rndNums = a .+ ((b - a) / (maximum(rndNums) - minimum(rndNums))) .* (rndNums .- minimum(rndNums))
+    rndNums
 end