In mathematics, Monte Carlo integration is a technique for numerical integration using random numbers. It is a particular Monte Carlo method that numerically computes a definite integral. While other algorithms usually evaluate the integrand at a regular grid, Monte Carlo randomly choose points at which the integrand is evaluated.

mc.swift

import Foundation

/// Approximates integral using Monte-Carlo method.
///
/// - parameter bounds: N-dimensional array of boundary limits.
/// - parameter f: integrand function.
/// - parameter N: number of points randomly picked to approximate.
///
/// - Returns: Approximated value of integral.
public func mcr(in bounds: NBounds, f: FuncND, N: Double) -> Double {
    let ps = (0..<Int(N)).map { _ -> NPoint in .random(in: bounds) }
    return bounds.map(-).reduce(1, *) * ps.map(f).reduce(0, +) / N
}

/// Approximates integral using Monte-Carlo method.
///
/// - parameter bounds: N-dimensional array of boundary limits.
/// - parameter f: integrand function.
/// - parameter N: number of points randomly picked to approximate.
///
/// - Returns: Approximated value of integral and estimated error.
public func mc(in bounds: NBounds, f: FuncND, N: Double) -> MCResult {
    let length = bounds.map(-).reduce(1, *)
    let ps = (0..<Int(N)).map { _ -> NPoint in .random(in: bounds) }
    
    let fAverage = ps.map(f).reduce(0, +) / N
    let result = length * fAverage
    
    let eAverage = ps.map { pow(f($0), 2) }.reduce(0, +) / N
    let error = length * sqrt((eAverage - pow(fAverage, 2)) / N)
    
    return (result, error)
}

mcLowDimInterfaces.swift

/// Approximates integral using Monte-Carlo method.
///
/// - parameter bounds: 1-dimensional array of boundary limits.
/// - parameter f: integrand function.
/// - parameter N: number of points randomly picked to approximate.
///
/// - Returns: Approximated value of integral.
public func mc(in bounds: Bounds, f: Func1D, N: Double) -> MCResult {
    return mc(in: [bounds], f: { f($0[0]) }, N: N)
}

/// Approximates integral using Monte-Carlo method.
///
/// - parameter bounds: 2-dimensional array of boundary limits.
/// - parameter f: integrand function.
/// - parameter N: number of points randomly picked to approximate.
///
/// - Returns: Approximated value of integral.
public func mc(in rect: Rect, f: Func2D, N: Double) -> MCResult {
    return mc(in: [(rect.a, rect.b), (rect.c, rect.d)], f: { f($0[0], $0[1]) }, N: N)
}

MonteCarloIntegration.swift

import Foundation

extension Double { static var e: Double { return 2.72 } }

let function: Func1D = sqrt
let bounds: Bounds = (0, 4)

// Using compact 1D interfaces 
// 1D
mc(in: bounds, f: function, N: 1000)
mc(in: (0, 1), f: { 4 / (1 + pow($0, 2)) }, N: 100)

// 2D
mc(in: (1, 2, 0, 2), f: { sin($0) + cos($1) }, N: 1000)


// Using N-dimensional interface
// 1D case
mc(in: [(-.pi / 4, 0)], f: { pow(tan($0[0]), 2) }, N: 1000)

// 2D case
mc(in: [(1.0 / 6, 1.0 / 3), (log(3), log(4))], 
   f: { 12 * $0[1] * pow(.e, 6 * $0[0] * $0[1]) }, 
   N: 1000
)

Types.swift

/// 1-dimensional bounds
public typealias Bounds = (a: Double, b: Double)
/// 2-dimensional bounds
public typealias Rect = (a: Double, b: Double, c: Double, d: Double)

/// MCPair - result and error approximation
public typealias MCResult = (result: Double, error: Double)

/// 1-dimensonal integrand function
public typealias Func1D = (Double) -> Double
/// 2-dimensonal integrand function
public typealias Func2D = (Double, Double) -> Double

/// N-dimensional bounds
public typealias NBounds = [Bounds]
/// N-dimensional point
public typealias NPoint = [Double]
/// N-dimensonal integrand function
public typealias FuncND = (NPoint) -> Double

public extension NPoint {
    /// Random N-dimensional point generator.
    ///
    /// - parameter bounds: array of bounds for each dimension.
    /// - Returns: Random point in **bounds**.
    static func random(in bounds: NBounds) -> NPoint {
        return (0..<bounds.count).map {  Double.random(in: bounds[$0].a...bounds[$0].b) }
    }
}