Keywords - Function groups - @ A B C D E F G H I J K L M N O P Q R S T U V W X Y Z

 Quantlet: mcint Description: plain Monte Carlo integration, computes an integral and its standard error (deviation) based on randomly distributed points with uniform distribution

 Usage: res = mcint(F,xl,xu,N) Input: F string, name of the function of n variables to be integrated. It should have one input - N x n matrix (randomly generated points) and one output - N x 1 vector of function values at the input points. xl n x 1 vector, vector of lower limits of integration xu n x 1 vector, vector of upper limits of integration N scalar, number of randomly distributed points, N >= 16*n Output: res.int scalar, the estimated integral of F over the integration region from xl to xu res.stderr scalar, the standard deviation of the estimated integral int of F over the integration region from xl to xu res.outtext string vector, text with the results (int and stderr)

Note:
This numerical method is designed only for a finite integration region (e.g. \int_a^b, a = - infinity or b = infinity, is not possible).

Example:
library("nummath")
; we integrate the constant function f = 1 over the region <0,2>
; using N = 10000 randomly distributed points
proc(p) = ftion(x)
p = 1 + 0.*x
endp
randomize(0)
I = mcint("ftion", 0, 2, 10000)
I


Result:
Contents of I.int
[1,]   2

Contents of I.stderr
[1,]   0

Contents of I.outtext
[1,] "  P L A I N   M O N T E   C A R L O  "
[2,] "_____________________________________"
[3,] "  Result   =       2.00000"
[4,] "  Error    =    0.0000e+00"

Example:
library("nummath")
; we integrate the function f = x^2 over the region <0,2>
; using N = 500000 randomly distributed points
proc(p) = ftion(x)
p = x^2
endp
randomize(0)
I = mcint("ftion", 0, 2, 500000)
I


Result:
Contents of I.int
[1,]   2.6667

Contents of I.stderr
[1,]  0.0023851

Contents of I.outtext
[1,] "  P L A I N   M O N T E   C A R L O  "
[2,] "_____________________________________"
[3,] "  Result   =       2.66667"
[4,] "  Error    =    2.3851e-03"

Example:
library("nummath")
; we integrate the function f = x^2 + 3*(y-1)^4 over
; the region <0,1> x <0,1> using N = 10000 randomly distributed points
proc(p) = ftion(x)
p = x[,1]^2 + 3*(x[,2]-1)^4
endp
randomize(0)
I = mcint("ftion", #(0,0), #(1,1), 10000)
I


Result:
Contents of I.int
[1,]  0.93561

Contents of I.stderr
[1,]  0.008528

Contents of I.outtext
[1,] "  P L A I N   M O N T E   C A R L O  "
[2,] "_____________________________________"
[3,] "  Result   =       0.93561"
[4,] "  Error    =    8.5280e-03"

Example:
library("nummath")
; we integrate the function f = x^2 * y^3 - exp(x - y) over
; the region <0,2> x <-1,3> using N = 250000 randomly distributed points
proc(p) = ftion(x)
p = x[,1]^2.*x[,2]^3 - exp(x[,1] - x[,2])
endp
randomize(0)
I = mcint("ftion", #(0,-1), #(2,3), 250000)
I


Result:
Contents of I.int
[1,]   35.962

Contents of I.stderr
[1,]  0.086717

Contents of I.outtext
[1,] "  P L A I N   M O N T E   C A R L O  "
[2,] "_____________________________________"
[3,] "  Result   =      35.96214"
[4,] "  Error    =    8.6717e-02"


Author: M. Svojik, L. Cizkova, 20030521 license MD*Tech
(C) MD*TECH Method and Data Technologies, 05.02.2006