Monte Carlo Integration

Share:

Description

A function to perform very simple Monte Carlo integration. This algorithm evaluates the integrand at a series of random parameter values drawn from a uniform distribution defined between the limits of integration. These limits must be finite in this implementation.

Usage

1
MCIntegration(func, lower, upper, initial.step = 10000, step.size = 10000, max.rel.var = .Machine$double.eps^0.25, max.runs = +Inf, ...)

Arguments

func

The integrand. The first argument of the integrand function must be a vector argument with an element for each parameter to be integrated over.

lower

A vector of values with each element containing the lower limit of integration for each parameter to be integrated. Limits must be finite.

upper

A vector of values with each element containing the upper limit of integration for each parameter to be integrated. Limits must be finite.

initial.step

An integer containing the initial number of integrand evaluations to perform before testing the series for a stopping condition. This value must be larger than one.

step.size

An integer containing the number of integrand evaluations to perform in every testing step. After the initial number of evaluations are made (set using the initial.step parameter) a test is made for the algorithm stopping condition every step.size iterations. This value must be larger than one.

max.rel.var

The maximum relative variance allowed for the integral estimate before the stopping condition is met. This value must be larger than zero.

max.runs

The maximum number of integrand evaluations to perform before the algorithm is stopped. This value must be larger than one.

...

Additional arguments to be passed to the integrand.

Details

The estimate of the integral Q is the product of the volume of integration and the sample mean of the series of integrand evaluations. The absolute variance of this estimate, Var(Q), is defined as

Var(Q) = V^2 Var(f) / N

where Var(f) is the sample variance of evaluations of the integrand. N is the sample size and V is the volume of integration. The relative variance is the absolute variance divided by the absolute value of the integral estimate, Q.

Value

A list containing the following elements:

val

The final estimate of the integral (Q).

abs.var

The absolute variance of the integral estimate. The relative variance is equal to abs.var / abs(val).

runs

The number of integrand evaluations performed before the stopping condition is met.

Note

This algorithm can take a long time to run, particularly when there is no upper limit to the number of integrand evaluations set. For univariate integration, adaptive quadrature methods such as those implemented in the integrate function may well perform better.

Author(s)

Joseph Chipperfield <joechip90@googlemail.com>

References

Kalos, M. H. and Whitlock, P. A. (2008) Monte Carlo methods (2nd edition), Wiley VCH.

See Also

LatticeTransitionProbs

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
# The limits of the integration we wish to perform
# for example between 3 and 6
lower.lim <- 3
upper.lim <- 6
# Parameters controlling the Monte Carlo integration algorithm
initial.step <- 10000	# Evaluate the integrand function 10000 times before testing algorithm stopping condition
step.size <- 10000	# Test stopping condition every 10000 integrand evaluations
# Stop when relative variance is less than .Machine$double.eps^0.25
max.rel.var <- .Machine$double.eps^0.25
# ... or stop when 100000 integrand evaluations have been performed
max.runs <- 100000
# Prepare a simple integrand function -
# we use here a function to square the
# input
sqrd <- function(x)
{
	x^2
}
# Perform the Monte Carlo integration
est <- MCIntegration(sqrd, lower.lim, upper.lim, initial.step, step.size, max.rel.var, max.runs)
# Function with analytical result of intration (indefinate integral)
sqrd.analytic.int <- function(x)
{
	x^3 / 3.0
}
# Calculate the analytic result
real.val <- sqrd.analytic.int(6) - sqrd.analytic.int(3)
# Calculate the difference between the analytic result and the estimate
est$val - real.val