Monte Carlo Integration

The Function class

The function to be integrated has its own datatype, the GSL::Monte::Function class.

Monte Carlo plans, alrgorithms

PLAIN Monte Carlo




Accessing internal state of the Monte Carlo classes

Miser Parameters (GSL-1.13 or later)

Accessors of Miser Params

This parameter introduces a random fractional variation of size dither into each bisection, which can be used to break the symmetry of integrands which are concentrated near the exact center of the hypercubic integration region. The default value of dither is zero, so no variation is introduced. If needed, a typical value of dither is 0.1.

Vegas Parameters (GSL-1.13 or later)

Accessors of Vegas Params


#!/usr/bin/env ruby
include GSL::Monte
include Math

proc_f = { |k, dim, params|
  pi = Math::PI
  a = 1.0/(pi*pi*pi)
  a/(1.0 - cos(k[0])*cos(k[1])*cos(k[2]))

def display_results(title, result, error)
  exact = 1.3932039296856768591842462603255

  diff = result - exact
  printf("%s ==================\n", title);
  printf("result = % .6f\n", result);
  printf("sigma  = % .6f\n", error);
  printf("exact  = % .6f\n", exact);
  printf("error  = % .6f = %.1g sigma\n", diff, diff.abs/error)

dim = 3
xl = Vector.alloc(0, 0, 0)
xu = Vector.alloc(PI, PI, PI)
G = Monte::Function.alloc(proc_f, dim)
calls = 500000
r = GSL::Rng.alloc(Rng::DEFAULT)

plain = Monte::Plain.alloc(dim)
result, error = G.integrate(xl, xu, dim, calls, r, plain)
display_results("plain", result, error)

miser = Monte::Miser.alloc(dim)
result, error = G.integrate(xl, xu, dim, calls, r, miser)
display_results("miser", result, error)

vegas = Monte::Vegas.alloc(dim)
result, error = G.integrate(xl, xu, dim, 10000, r, vegas)
display_results("vegas warm-up", result, error)
  result, error = G.integrate(xl, xu, dim, calls/5, r, vegas)
  printf("result = % .6f sigma = % .6f chisq/dof = %.1f\n",
          result, error, vegas.chisq)
end while (vegas.chisq-1.0).abs > 0.5
display_results("vegas final", result, error)

prev next

Reference index top