# Monte Carlo Integration

### 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 |

`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 ( |

`abs.var` |
The absolute variance of the integral estimate. The relative variance is equal to |

`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
``` |