24 Integração de Monte Carlo

Esse capítulo descreve rotinas para integração multidimensional de Monte Carlo. Esse capítulo inclui o método tradicional de Monte Carlo e algoritmos adaptativos tais como VEGAS e MISER que usam amostragem por importância e técnicas de amostragem estratificada. Cada algoritmo calcula uma estimativa de uma integral definida multidimensional da forma,

As funções são definidas em arquivos de cabeçalho separados para cada rotina, ‘gsl_monte_plain.h’, ‘gsl_monte_miser.h’ e ‘gsl_monte_vegas.h’.

24.1 Interface

Todas as rotinas de integração de Monte Carlo usam a mesma forma geral de interface. Existe um alocador para alocar memória para variáveis de controle e espaço de trabalho, uma rotina para inicializar essas variáveis de controle, o integrador propriamente dito, e uma função para liberar o espaço quando tudo acabar.

Cada função de integração requer um gerador de números aleatórios para ser fornecido, e retorna uma estimativa da integral e seu desvio padrão. A precisão do resultado é determinada pelo número de chamadas de função especificada pelo usuário. se um nível conhecido de precisão é requerido isso pode ser encontrado chamado o integrador muitas vezes e avaliando os resultados individuais até que a precisão desejada seja obtida.

Pontos de amostra aleatória usados dentro das rotinas de Monte Carlo são sempre escolhidos estritamente dentro da região de integraço, De forma que singularidades de extremidade são automaticamente evitadas.

A função a ser integrada tem seu próprio tipo de dado, definido no arquivo de cabeçalho ‘gsl_monte.h’.

Data Type: gsl_monte_function

Esse tipo de dado define uma função geral com parâmetros para integração de Monte Carlo.

double (* f) (double * x, size_t dim, void * params)

essa função deve retornar o valor f(x,params) para o argumento x e parâmetros params, onde x é um vetor estático de tamanho dim fornecendo as coordenadas do ponto onde a função é para ser avaliada.

size_t dim

o número de dimensões para x.

void * params

um apontador para os parâmetros da função.

Aqui está um exemplo para uma função quadrática em duas dimensões,

struct my_f_params { double a; double b; double c; };

double
my_f (double x[], size_t dim, void * p) {
   struct my_f_params * fp = (struct my_f_params *)p;

   if (dim != 2)
      {
        fprintf (stderr, "error: dim != 2");
        abort ();
      }

   return  fp->a * x[0] * x[0] 
             + fp->b * x[0] * x[1] 
               + fp->c * x[1] * x[1];
}

gsl_monte_function F;
struct my_f_params params = { 3.0, 2.0, 1.0 };

F.f = &my_f;
F.dim = 2;
F.params = &params;

A função f(x) pode ser avaliada usando a seguinte macro,

#define GSL_MONTE_FN_EVAL(F,x) 
    (*((F)->f))(x,(F)->dim,(F)->params)

24.2 Monte Carlo Tradicional

O algoritmo tradicional de Monte Carlo faz amostragem de pontos aleatóriamente a partir da região de integração para estimar a integral e seu erro. Usando esse algoritmo a estimativa da integral E(f; N) para N pontos aleatóriamente distribuídos x_i é dada por,

As funções descritas nessa seção são declaradas no arquivo de cabeçalho ‘gsl_monte_plain.h’.

Function: gsl_monte_plain_state * gsl_monte_plain_alloc (size_t dim)

Essa função aloca e inicializa um espaço de trabalho para uma integração de Monte Carlo em dim dimensões.

Function: int gsl_monte_plain_init (gsl_monte_plain_state* s)

Essa função inicializa um estado de integração anteriormente alocado. Isso permite que um espaço de trabalho existente seja reutilizado para diferentes integrações.

Function: int gsl_monte_plain_integrate (gsl_monte_function * f, const double xl[], const double xu[], size_t dim, size_t calls, gsl_rng * r, gsl_monte_plain_state * s, double * result, double * abserr)

Essas rotinas usam o algoritmo tradicional de Monte Carlo para integrar a função f sobre a região hipercúbica dim-dimensional definida pelos limites inferior e superior nos vetores estáticos xl e xu, cada um de tamanho dim. A integração usa um número fixo de chamadas a função calls, e obtém pontos de amostragem aleatórios usando o gerador de números aleatórios r. Um espaço de trabalho previamene alocado s deve ser aplicado. O resultado da integração é retornado em result, com um erro absoluto estimado abserr.

Function: void gsl_monte_plain_free (gsl_monte_plain_state * s)

Essa função libera a memória associada ao estado do integrador s.

24.3 MISER

O algorittmo MISER de Press e Farrar é baseado sobre amostragem recursiva estratificada. Essa técnica pretende reduzir o erro geral de integração concentrando os pontos de integração em regiões de grande variância.

A idéia da amostragem estratificada inicia-se com a observação de que para duas regiões disjuntas a e b com estimativas de integrais feitas através do método Monte Carlo E_a(f) e E_b(f) e variâncias \sigma_a^2(f) e \sigma_b^2(f), a variância \Var(f) da estimativa combinada E(f) = (1/2) (E_a(f) + E_b(f)) é dada por,

O algoritmo de MISER continua através da bissecção da região de integração em torno de um eixo coordenado para fornecer duas sub-regiões a cada degrau. A direção é escolhida examinando todas as d possíveis bissecções e selecionando aquela que irá minimizar a variância combinada das duas sub-regiões. A variância nas sub-regiões é estimada por amostragem com uma fração do número total de pontos disponíveis para o passo atual. O mesmo procedimento é então repetido recursivamente para cada um dos dois meio espaços a partir da melhor bissecção. Os restantes pontos a amostra são alocados às sub-regiões usando a fórmula para N_a e N_b. Essa alocação recursiva de pontos de integração continua convergindo para um nível definido pelo usuário onde cada sub-região é integrada usando uma estimativa pelo método original de Monte Carlo. Esses valores individuais e suas estimativas de erro são então combinados para cima de forma a fornecer um resultado geral e uma estimativa de seu erro.

As funções descritas nessa seção são declaradas no arquivo de cabeçalho ‘gsl_monte_miser.h’.

Function: gsl_monte_miser_state * gsl_monte_miser_alloc (size_t dim)

Essa função aloca e inicializa um espaço de trabalh para integração de Monte Carlo em dim dimensões. O espaço de trabalho é usado para manter o estado de integração.

Function: int gsl_monte_miser_init (gsl_monte_miser_state* s)

Essa função inicializa um previamente alocado estado de integração. Isso permite que um espaço de trabalho existente seja reutilizado por diferentes integrações.

Function: int gsl_monte_miser_integrate (gsl_monte_function * f, const double xl[], const double xu[], size_t dim, size_t calls, gsl_rng * r, gsl_monte_miser_state * s, double * result, double * abserr)

Essa rotina usa o algoritmo MISER Monte Carlo para integrar a função f sobre a região hipercúbica dim-dimensional definida pelos limites inferior e superior nos vetores estáticos xl e xu, cada um de tamanho dim. A integração usa um número fixo de chamadas a funções calls, e obtém pontos de amostragem aleatórios usando o gerador de números aleatórios r. Um espaço de trabalho alocado previamente s deve ser fornecido. O resultado da integração é retornado em result, com um erro absoluto estimado abserr.

Function: void gsl_monte_miser_free (gsl_monte_miser_state * s)

Essa função libera a memória associada ao estado de integrador s.

O algoritmo MISER tem muitos parâmetros configuráveis que podem ser modificados usando as seguintes duas funções.(47)

Function: void gsl_monte_miser_params_get (const gsl_monte_miser_state * s, gsl_monte_miser_params * params)

Essa função copia os parâmetros do estado do integrador para a estrutura params fornecida pelo usuário.

Function: void gsl_monte_miser_params_set (gsl_monte_miser_state * s, const gsl_monte_miser_params * params)

Essa função ajusta os parâmetros do integrador baseado em valores fornecidos na estrutura params.

Tipicamente os valores dos parâmetros são primeiramente lidos usando gsl_monte_miser_params_get, as mudanças necessárias são feitas para os campos da estrutura params, e os valores são copiados de volta ao estado do integrador usando gsl_monte_miser_params_set. As fuções usam a estrutura gsl_monte_miser_params que conté os seguintes campos:

Variable: double estimate_frac

Esse parâmetro especifica a fração do número atualmente disponível de chamadas a funções que são slocados para estimar a variância em cada passo recursivo. O valor padrão é 0.1.

Variable: size_t min_calls

esse parâmetro especifica o número mínimo de chamadas a funções requeridas para cada estimativa de variância. Se o número de chamadas a funções alocadas para a estimativa usando estimate_frac encontra-se abaixo de min_calls então min_calls é usada em seu lugar. Isso garante que cada estimativa mantenha uma nível rasoável de precisão. O valor padrão de min_calls é 16 * dim.

Variable: size_t min_calls_per_bisection

Esse parâmetro especifica o número mínimo de chamadas a funções requeridos para proceder a um passo de bissecção. Quando um passo recursivo tem poucas chamadas disponíveis que min_calls_per_bisection executa uma estimativa original de Monte Carlo da atual sub-região e termina seu ramo de recursão. O valor padrão desse parâmetro é 32 * min_calls.

Variable: double alpha

Esse parâmetro controla como as variâncias estimadas para as duas sub-regiões de uma bissecção são copiadas quando se faz alocação de pontos. Com amostragem recursiva a variância média deve ajustar-se proporcionalmente melhor que 1/N, uma vez que os valores a partir de sub-regiões irão ser obtidos usando um procedimento que explicitamente minimiza sua variância. Para acomodar esse comportamento o algoritmo MISER permite que a variância total dependa um parâmetro de ajuste proporcional \alpha,

Var(f) = σa Naα + σb Nbα .

Os autores do artigo original descrevendo MISER recomendam o valor \alpha = 2 como uma boa escolha, obtido a partir de experimentos numéricos, e esse valor é usado como o valor padrão nessa implementação.

Variable: double dither

Esse parâmetro introduz uma variação fracionária alearória de tamanho dither em cada bissecção, que pode ser usada para parar a simetria de integrandos que estejam concentrados próximo do centro exato da região hipercúbica de integração. O valor padrão do tamanho dither é zero, de forma que nenhuma variação é introduzida. Se necessário, um valor típico de dither is 0.1.

24.4 VEGAS

O algoritmo VEGAS de Lepage é baseado em amostragem por importância. Esses pontos de amostra a partir da distribuição de probabilidade descrita pela função |f|, de forma que os pontos são concentrados nas regiões que fazem a maior contribuição para a integral.

Em geral, se a integral de Monte Carlo de f recebe amostra com pontos distribuídos conforme a distribuição de probabilidade descrita pela função g, obtemos uma estimativa E_g(f; N),

O algoritmo VEGAS aproxima-se da distribuição exata fazendo um número de passagens sobre a região de integração enquanto faz o histrograma da função f. Cada histograma é usado para definir uma distribuição de amostragem para a passagem seguinte. Assintóticqamente esse procedimento converge para a distribuição desejada. Com o objetivo de evitar que o número de bins de histograma aumente como K^d a distribuição de probabilidade é aproximada por uma função separada: g(x_1, x_2, ...) = g_1(x_1) g_2(x_2) ... de forma que o número de bins requerido é somente Kd. Isso é equivalente a localizar as seleções da função a partir das projeções do integrando sobre os eixos coordenados. A eficiência de VEGAS depende da validade dessa premissa. é mais eficiente quando a seleção do integrando são bem localizados. Se um integrando pode ser re-escrito de forma que seja aproximadamente separável isso pode melhorar a eficiência da integração com VEGAS.

O algoritmo VEGAS incorpora um número de recursos adicioonais, e combina ambas amostragem por estratificação e amostragem por imprtância. A região de integração é dividida em um número de “caixas”, com cada caixa pegando um número fixo de pontos (a meta é 2). Cada caixa pode então ter um número de bins fracional, mas se a razão de bins-por-caixa for menor que dois, Vegas alterna para um tipo de redução de variância (em lugar da amostragem por importância).

Function: gsl_monte_vegas_state * gsl_monte_vegas_alloc (size_t dim)

Essa função aloca e inicializa um espaço de trabalho para a integração de Monte Carlo em dim dimensões. O espaço de trabalho é usado para manter o estado de integração.

Function: int gsl_monte_vegas_init (gsl_monte_vegas_state* s)

Essa função inicializa um estado de integração previamente alocado. Isso permite que um espaço de trabalho existente seja reutilizado por diferentes integrações.

Function: int gsl_monte_vegas_integrate (gsl_monte_function * f, double xl[], double xu[], size_t dim, size_t calls, gsl_rng * r, gsl_monte_vegas_state * s, double * result, double * abserr)

Essa rotina usa o algoritmo VEGAS Monte Carlo para integrar a função f sobre a região hipercúbica dim-dimensional definida pelos limites inferior e superior nos vetores estáticos xl e xu, cada um de tamanho dim. A integração utiliza um número fixo de chamadas a funções calls, e obtém pontos aleatórios de amostragem usando o gerador de números aleatórios r. Um espaço de trabalho alocado previamente s deve ser fornecido. O resultado da integração é retornado em result, com um erro absoluto estimado abserr. O resultado e seu erro estimado são baseados em um peso médio de amostras independentes. O chi-quadrado por grau de liberdade para o peso médio é retornado via componente da estrutura de estado, s->chisq, e deve ser consistente com 1 para peso médio ser confiável.

Function: void gsl_monte_vegas_free (gsl_monte_vegas_state * s)

essa função libera a memória associada ao estado do integrador s.

O algoritmo VEGAS calcula um número de estimativas independentes da integral internamente, conforme o parâmetro iterations descrito abaixo, e retorna seu peso médio. Amostragem aleatória do integrando pode ocasionalmente produzir uma estimativa onde o erro é zero, particularmente se a função for constante em algumas regiões. Uma estimative com erro zero faz com que o peso médio pare abaixo do valor e deve ser tratada separadamente. Na implementações fortran originais de VEGAS o erro estimado é feito diferente de zero por substituição por um pequeno valor (typicamente 1e-30). A implementação na GSL difere dessas implementações fortran e evita o uso de uma constrante arbitrária—isso ou atribui o valor um peso que é a média ponderada do procedimento estimado ou descarta-a conforme o seguinte procedimento,

estimativa atual tem erro zero, média ponderada tem ero finito

A atual estimativa é atribuída um peso que é peso médio das estimativas do procedimento.

esstimativa atual tem erro finito, estimativa anterior tem erro zero

As estimativas anteriores são descartadas e o procedimento do peso médio começa com a estimativa atual.

estimativa atual tem erro zero, estimativa anterior também tinha erro zero

As estimativas teem sua média calculada usando a média artimétca, mas nenhum erro é calculado.

A convergência do algoritmo pode ser testada usando o valor médio de chi-quadrado dos resultados, que estão disponíveis a partir da seguinte função:

Function: double gsl_monte_vegas_chisq (const gsl_monte_vegas_state * s)

Essa função retorna o chi-quadrado por grau de liberdade para a estimativa de peso da integral. O valor retornado deve estar perto de 1. Um valor que difere significativamente de 1 indica que os valores de diferentes iterações são inconsistentes. Ness caso o erro de peso irá ser estimado para baixo, e iterações adicionais do algoritmo são necessárias para obter resultados confiáveis.

Function: void gsl_monte_vegas_runval (const gsl_monte_vegas_state * s, double * result, double * sigma)

Essa função retorna os valores de linha (não incluídos na média) da integral result e seu erro sigma a partir da mais recente iteração do algoritmo.

O algoritmo VEGAS é altamente configurável. Muitos parâmetros podem ser mudados usando as seguintes duas funções.

Function: void gsl_monte_vegas_params_get (const gsl_monte_vegas_state * s, gsl_monte_vegas_params * params)

Essa função copia os parâmetros de estado do integrador para a estrutura params fornecida pelo usuário.

Function: void gsl_monte_vegas_params_set (gsl_monte_vegas_state * s, const gsl_monte_vegas_params * params)

Essa função ajusta os parâmetros do integrador baseada nos valores fornecidos na estrutura params.

Tipicamente os valores dos parâmetros são primeiro lidos usando gsl_monte_vegas_params_get, as mudanças necessárias são feitas nos campos da estrutura params, e os valores são copiados de volta para o estado do integrador usando gsl_monte_vegas_params_set. As funções usam a estrutura gsl_monte_vegas_params que contém os seguintes campos:

Variable: double alpha

O parâmetro alpha controla a rigidez do algoritmo de de rebinning (48). É tipicamente ajustado entre um e dois. Um valor zero evita rebinning da grade. O valor padrão é 1.5.

Variable: size_t iterations

O número de iterações para executar para cada chamada para a rotina. O valor padrão é 5 iterações.

Variable: int stage

Ajustando esse parâmetro determinamos o estágio stage do cálculo. Normalmente, stage = 0 que começa com uma nova grade uniforme e peso médio vazio. Chamando VEGAS com stage = 1 retemos a grade da execução anterior mas descartamos o peso médio, de forma que se pode ajustar “ajustar” a grade usando um número relativamente pequeno de pontos e então fazer uma grande execução com stage = 1 sobre a grade otimizada. Ajustando stage = 2 mantém-se a grade e o peso médio da execução anterior, mas pode-se aumentar (ou decrescer) o número de bins de histograma na grade dependendo do número de chamadas disponíveis. Escolhendo stage = 3 entramos no laço principal, de forma que nada é modificado, e é equivalente a executar iterações adicionais usando os dados de uma chamada anterior.

Variable: int mode

As escolhas possíveis são GSL_VEGAS_MODE_IMPORTANCE, GSL_VEGAS_MODE_STRATIFIED, GSL_VEGAS_MODE_IMPORTANCE_ONLY. Determina se VEGAS irá usar amostragem por improtância ou amostragem estratificada, ou se é possível selecionar sobre si mesmo. Em pequenas dimensões VEGAS usa amostragem estritamente estratificada (mais precisamente, amostragem estratificada é escolhida se houverem menos que 2 bins por caixa).

Variable: int verbose

Variable: FILE * ostream

Esses parâmetros ajustam o nível de informação mostrado por VEGAS. Toda informação é escrita para o fluxo ostream. O ajuste padrão de verbose é -1, que desliga toda saída. Um valor para verbose sendo 0 mostra informação resumida sobre o peso médio e resultado final, enquanto um valor de 1 também mostra as coordenadas da grade. Um valor de 2 mostra informações do procedimento de rebinning para cada iteração.

Os campos acima e o valor chisq podem também serem acessados diretamente em gsl_monte_vegas_state mas a utilização dessa forma está desatualizada.

24.5 Exemplos

O programa exemplo abaixo usa as rotinas de Monte Carlo para estimar o valor das seguintes integrais tridimensionais partindo da teoria de caminhadas aleatórias,

Por simplicidade iremos calcular a integral sobre a região de (0,0,0) a (\pi,\pi,\pi) e multiplicar por 8 para obter o resultado completo. A integral está variando lentamente no meio da região mas tem singularidades integráveis nos cantos (0,0,0), (0,\pi,\pi), (\pi,0,\pi) e (\pi,\pi,0). As rotinas de Monte Carlo somente selecionam pontos que estão estritamente dentro da região de integração e de forma que nenhuma medida especial é preciso para evitar essas singularidades.

#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_miser.h>
#include <gsl/gsl_monte_vegas.h>

/* Computation of the integral,

      I = int (dx dy dz)/(2pi)^3  1/(1-cos(x)cos(y)cos(z))

   over (-pi,-pi,-pi) to (+pi, +pi, +pi).  The exact answer
   is Gamma(1/4)^4/(4 pi^3).  This example is taken from
   C.Itzykson, J.M.Drouffe, "Statistical Field Theory -
   Volume 1", Section 1.1, p21, which cites the original
   paper M.L.Glasser, I.J.Zucker, Proc.Natl.Acad.Sci.USA 74
   1800 (1977) */

/* For simplicity we compute the integral over the region 
   (0,0,0) -> (pi,pi,pi) and multiply by 8 */

double exact = 1.3932039296856768591842462603255;

double
g (double *k, size_t dim, void *params)
{
  double A = 1.0 / (M_PI * M_PI * M_PI);
  return A / (1.0 - cos (k[0]) * cos (k[1]) * cos (k[2]));
}

void
display_results (char *title, double result, double error)
{
  printf ("%s ==================\n", title);
  printf ("result = % .6f\n", result);
  printf ("sigma  = % .6f\n", error);
  printf ("exact  = % .6f\n", exact);
  printf ("error  = % .6f = %.2g sigma\n", result - exact,
          fabs (result - exact) / error);
}

int
main (void)
{
  double res, err;

  double xl[3] = { 0, 0, 0 };
  double xu[3] = { M_PI, M_PI, M_PI };

  const gsl_rng_type *T;
  gsl_rng *r;

  gsl_monte_function G = { &g, 3, 0 };

  size_t calls = 500000;

  gsl_rng_env_setup ();

  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

  {
    gsl_monte_plain_state *s = gsl_monte_plain_alloc (3);
    gsl_monte_plain_integrate (&G, xl, xu, 3, calls, r, s, 
                               &res, &err);
    gsl_monte_plain_free (s);

    display_results ("plain", res, err);
  }

  {
    gsl_monte_miser_state *s = gsl_monte_miser_alloc (3);
    gsl_monte_miser_integrate (&G, xl, xu, 3, calls, r, s,
                               &res, &err);
    gsl_monte_miser_free (s);

    display_results ("miser", res, err);
  }

  {
    gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (3);

    gsl_monte_vegas_integrate (&G, xl, xu, 3, 10000, r, s,
                               &res, &err);
    display_results ("vegas warm-up", res, err);

    printf ("converging...\n");

    do
      {
        gsl_monte_vegas_integrate (&G, xl, xu, 3, calls/5, r, s,
                                   &res, &err);
        printf ("result = % .6f sigma = % .6f "
                "chisq/dof = %.1f\n", res, err, gsl_monte_vegas_chisq (s));
      }
    while (fabs (gsl_monte_vegas_chisq (s) - 1.0) > 0.5);

    display_results ("vegas final", res, err);

    gsl_monte_vegas_free (s);
  }

  gsl_rng_free (r);

  return 0;
}

Com 500,000 chamadas a funções o algoritmo de Monte Carlo original encontra um erro fracional de 1%. O erro estimado sigma é grosseiramente consistente com o erro atual–o resultado computado difere do verdadeiro resultado em torno do valor de 1.4 de desvios padrões,

plain ==================
result =  1.412209
sigma  =  0.013436
exact  =  1.393204
error  =  0.019005 = 1.4 sigma

O algoritmo MISER reduz o erro por uma fator de quatro, e também corretamente estima o erro,

miser ==================
result =  1.391322
sigma  =  0.003461
exact  =  1.393204
error  = -0.001882 = 0.54 sigma

No caso do algoritmo VEGAS o programa usa uma estimativa inicial de execução de 10,000 chamadas a funções para preparar, ou “esquentar”, a grade. Isso é seguido por uma execução principal com cindo iterações de 100,000 chamadas a funções. O chi-quadrado por grau de liberdade para as cinco iterações são verificadas para a consistência com 1, e a execução é repetida se os resultados não tiveem convergindo. Nesse caso a estimativa está consistente sobre o pimeiro passo.

vegas warm-up ==================
result =  1.392673
sigma  =  0.003410
exact  =  1.393204
error  = -0.000531 = 0.16 sigma
converging...
result =  1.393281 sigma =  0.000362 chisq/dof = 1.5
vegas final ==================
result =  1.393281
sigma  =  0.000362
exact  =  1.393204
error  =  0.000077 = 0.21 sigma

Se o valor de chisq diferir significativamente de 1 pode indicar resultados inconsistententes, com um correspondentemente erro estimado baixo. A estimativa final de VEGAS (usando um número similar de chamadas a funções) é significativamente mais preciso que os outros dois algoritmos.

24.6 Referências e Leituras Adicionais

O algoritmo de MISER é descrito no seguinte artigo de autoria de Press e de Farrar,

• W.H. Press, G.R. Farrar, Recursive Stratified Sampling for Multidimensional Integração de Monte Carlo, Computers in Physics, v4 (1990), pp190–195.

O algoritmo VEGAS é descrito nos seguintes artigos,

• G.P. Lepage, A New Algorithm for Adaptive Multidimensional Integration, Journal of Computational Physics 27, 192–203, (1978)
• G.P. Lepage, VEGAS: An Adaptive Multi-dimensional Integration Program, Cornell preprint CLNS 80-447, March 1980

Esse documento foi gerado em 23 de Julho de 2013 usando texi2html 5.0.