Métodos de amostragem MCMC - Algoritmo de Metropolis-Hastings

Introdução

Os métodos de Monte Carlo via cadeias de Markov (MCMC) constituem uma classe de algoritmos de amostragem que permite extrair amostras de uma distribuição-alvo complexa p(x). O objetivo do MCMC é obter um conjunto de amostras que represente com precisão essa distribuição-alvo, permitindo estimar médias, variâncias e outras características da distribuição. A essência do MCMC está na construção de uma cadeia de Markov específica cuja distribuição estacionária coincide com a distribuição-alvo.

Esses métodos são amplamente utilizados na inferência bayesiana, em machine learning e em outras áreas que exigem a aproximação de distribuições a posteriori. Diferentemente de métodos determinísticos, como inferência variacional ou aproximação de Laplace, os algoritmos MCMC possuem uma propriedade única: quando configurados corretamente, garantem assintoticamente a obtenção de amostras que seguem exatamente a distribuição-alvo. Entre os diversos algoritmos MCMC, o algoritmo de Metropolis-Hastings (MH) ocupa um lugar de destaque, sendo um método fundamental que serve de base para muitas abordagens modernas.

Neste artigo, analisaremos o algoritmo de Metropolis-Hastings começando por seus fundamentos teóricos e conceitos principais. Em seguida, apresentaremos sua implementação em MQL5 na forma da classe MHSampler e veremos exemplos simples de aplicação a distribuições unidimensionais e multidimensionais, incluindo as variantes Random Walk Metropolis e Metropolis-Hastings Independente.

Vale destacar que, ao utilizar métodos MCMC, é preciso dedicar atenção especial ao diagnóstico e ao ajuste do algoritmo, pois seu funcionamento correto exige uma análise cuidadosa da convergência da cadeia com o auxílio de ferramentas estatísticas e visuais.

Algoritmo de Metropolis

O algoritmo básico de Metropolis (Metropolis, 1953) gera amostras da distribuição-alvo propondo, a cada etapa, uma transição do estado atual x para um novo estado x′ com probabilidade q(x′∣x), em que q é chamada de distribuição auxiliar (distribuição de proposta, proposal distribution). O principal requisito para a distribuição de proposta é que ela garanta que a cadeia não fique presa em uma única região, ou seja, que seja possível transitar de qualquer ponto do suporte da distribuição-alvo para qualquer outro ponto. Para isso, basta que a probabilidade q(x′∣x) seja diferente de zero para todos os estados possíveis x e x′.

O estado proposto x′ é aceito ou rejeitado segundo um critério que garante que, no longo prazo, a cadeia visite os estados proporcionalmente a p(x). Se o estado x′ for aceito, ele se torna o novo estado; caso contrário, a cadeia permanece no estado atual, o que resulta na repetição da amostra.

O algoritmo de Metropolis é um caso particular do algoritmo de Metropolis-Hastings, no qual a distribuição de proposta q(x'|x) é simétrica, ou seja, a probabilidade de propor uma transição de x para x′ é igual à probabilidade da transição inversa: q(x′∣x)=q(x∣x′).

Graças a essa simetria, o critério de aceitação do candidato x′ fica consideravelmente mais simples:

A = min (1, p*(x')/p*(x))

onde p*(x) é a densidade não normalizada da distribuição-alvo, enquanto a própria distribuição é definida como p(x)=p*(x)/Zp.

A constante de normalização Zp pode ser desconhecida, o que torna o algoritmo especialmente conveniente em problemas nos quais o cálculo analítico de Zp é difícil, como ocorre frequentemente na inferência bayesiana.

Se a densidade no novo ponto x′ for maior do que no ponto atual x, isto é, p*(x′) > p*(x), então a probabilidade de aceitação A=1 e a transição para x′ ocorre sempre, com probabilidade 1. Isso permite que a cadeia avance rapidamente em direção aos picos da distribuição, favorecendo a busca por regiões de alta probabilidade.

Se a densidade diminuir, p*(x′) < p*(x), a transição para x′ ocorre com probabilidade A = p*(x′)/p*(x). Esse mecanismo permite que a cadeia se afaste dos máximos da distribuição e explore suas encostas e vales, percorrendo de forma eficiente todo o espaço e evitando que fique permanentemente presa em um único ponto.

É justamente esse equilíbrio entre a aceitação incondicional de movimentos "para cima" e a aceitação probabilística de movimentos "para baixo" que garante que a cadeia preserve as proporções da distribuição-alvo, permanecendo mais tempo nas regiões em que a densidade é maior.

A sequência de amostras obtida dessa forma constitui uma cadeia de Markov, na qual cada estado seguinte depende do anterior. Isso gera autocorrelação entre as amostras, o que pode reduzir sua representatividade, pois, quando a autocorrelação é alta, amostras vizinhas contêm praticamente a mesma informação.

Por exemplo, se você coletar 1000 amostras fortemente correlacionadas, elas podem conter a mesma quantidade de informação que apenas 50 amostras independentes. Nesse caso, as 1000 amostras não representam 1000 observações únicas da distribuição, mas apenas cerca de 50 observações efetivas. Como consequência, as estimativas da média, da variância e de outras características da distribuição obtidas com essas amostras correlacionadas terão maior variância, ou erro, do que teriam se fosse utilizado o mesmo número de amostras independentes. Por isso, para reduzir a correlação, utiliza-se o chamado período de aquecimento (burn-in), que descarta as iterações iniciais, e o afinamento (thinning), que mantém apenas cada i-ésima amostra.

Fig. 1 Algoritmo de Metropolis

Cadeias de Markov

Os métodos de Monte Carlo via cadeias de Markov utilizam cadeias de Markov para gerar amostras de distribuições complexas. Para entender como isso funciona, vamos analisar duas propriedades fundamentais das cadeias de Markov que devem ser satisfeitas para construir um algoritmo MCMC: sua "memória" e sua capacidade de convergir para a distribuição-alvo.

1. Propriedade de Markov

Uma cadeia de Markov é uma sequência de estados x(1), x(2), …, na qual cada estado seguinte x(t+1) depende apenas do estado atual x(t) e não de todo o histórico anterior da cadeia x(1), …, x(t−1):

P(x(t+1) |x(t), …, x(1)) = P(x(t+1) |x(t))

Essa propriedade significa que, para prever o estado futuro, basta conhecer apenas o estado atual.

2. Convergência para a distribuição-alvo

O objetivo dos algoritmos MCMC é construir uma cadeia de Markov que, ao longo do tempo, após o período inicial de aquecimento, ou burn-in, gere amostras da distribuição-alvo p(x). Isso é obtido por meio de uma regra específica de transição que satisfaz a condição de equilíbrio detalhado (Detailed Balance):

p(x)T(x′∣x) = p(x′)T(x∣x′)

onde:

• p(x) e p(x′) são as densidades da distribuição-alvo,
• T(x′∣x) é a probabilidade total de transição de x para x′.

No algoritmo de Metropolis-Hastings, essa importante condição é garantida pela escolha de uma probabilidade de aceitação A adequada, que faz com que, ao longo do tempo, a cadeia "visite" os estados proporcionalmente a p(x), mesmo quando a distribuição é conhecida apenas parcialmente, sem a constante de normalização Zp. Nesse caso, a probabilidade total de transição T(x′∣x) é formada pelo produto da probabilidade de proposta q(x′∣x) pela probabilidade de aceitação A(x′∣x).

Algoritmo de Metropolis-Hastings

O algoritmo de Metropolis-Hastings (Hastings, 1970) generaliza o algoritmo de Metropolis ao permitir o uso de distribuições de proposta assimétricas, isto é, q(x'|x) ≠ q(x|x'). Isso torna o algoritmo mais flexível e aplicável a uma ampla variedade de problemas nos quais distribuições simétricas não são adequadas, por exemplo, quando são usadas propostas lognormais ou exponenciais.

Para compensar essa assimetria, o critério de aceitação A incorpora a correção de Hastings, que corrige o viés introduzido pela distribuição de proposta assimétrica. A probabilidade de aceitação do candidato x′ é definida como:

A = min (1, p*(x')q(x|x') / p*(x)q(x'|x))

onde:

• p*(x) é a densidade não normalizada da distribuição-alvo,
• q(x|x') / q(x'|x) é a correção de Hastings, que leva em conta a assimetria da distribuição de proposta.

Para calcular A, assim como antes, não é necessário conhecer a constante de normalização Zp da distribuição de probabilidade p(x) = p*(x)/Zp, pois ela se cancela.

Influência da distribuição de proposta no desempenho

A escolha de q(x'|x) tem grande impacto sobre a eficiência do algoritmo. Em espaços contínuos, é comum utilizar uma distribuição normal centrada no estado atual x(t), o que resulta no comportamento conhecido como Random Walk Metropolis-Hastings (RWMH). Nesse caso, o parâmetro de variância, ou passo da cadeia, desempenha um papel fundamental, pois determina a escala das transições propostas:

• Variância pequena: a taxa de aceitação (acceptance rate) será alta, mas a cadeia avançará lentamente. Isso resultará em um longo tempo de autocorrelação das amostras devido à exploração lenta do espaço.
• Variância grande: a taxa de aceitação será baixa, pois a maioria dos passos propostos cairá em regiões de baixa probabilidade. Isso leva a rejeições frequentes e ao uso ineficiente dos recursos computacionais.

Na prática, o melhor desempenho é obtido quando a taxa de aceitação fica na faixa de 20 a 40%. Isso exige um ajuste cuidadoso dos parâmetros de q para equilibrar a velocidade de exploração do espaço e a probabilidade de aceitação dos candidatos.

Fig. 2 Algoritmo de Metropolis-Hastings

O algoritmo de Metropolis-Hastings foi implementado em MQL5 com uma abordagem orientada a objetos. A classe principal MHSampler exige a implementação de três classes abstratas: LogPDF, para o logaritmo da densidade-alvo, PropRND, para gerar um novo estado, e LogPropPDF, para o logaritmo da densidade de transição. Graças a essa estrutura, o algoritmo pode ser facilmente configurado para qualquer distribuição-alvo e distribuição de proposta.

Classe LogPDF

A classe LogPDF define a interface para calcular o logaritmo da densidade de probabilidade da distribuição-alvo p*(x). O método LogPdf(const vector &x) recebe o vetor x, que representa um ponto no espaço de estados, e retorna ln p*(x).

O uso do logaritmo da densidade é uma prática padrão em MCMC, pois evita problemas numéricos ao trabalhar com probabilidades muito pequenas, especialmente em problemas multidimensionais. A classe pode ser implementada para qualquer distribuição-alvo, por exemplo, uma distribuição normal ou uma mistura de distribuições. Na inferência bayesiana, essa classe é implementada para calcular o logaritmo da densidade a posteriori não normalizada, expresso como a soma dos logaritmos: ln p*(x) = ln(verossimilhança) + ln(distribuição a priori).

Classe PropRND

A classe PropRND é responsável por gerar estados candidatos aleatórios x′ a partir da distribuição de proposta q(x′∣x). O método PropRnd(const vector &x) recebe o estado atual x e retorna um novo vetor x′ gerado a partir da distribuição q(x′∣x).

Por exemplo, no Random Walk MH, o usuário pode implementar PropRND com uma distribuição gaussiana centrada no estado atual x. Já no Independent MH, pode usar uma distribuição fixa que não depende de x.

Classe LogPropPDF

A classe LogPropPDF é responsável por calcular o logaritmo da densidade da distribuição de proposta q(x′∣x). Essa classe não gera novos estados, apenas avalia a densidade de probabilidade da transição. O método LogPropPdf(const vector &x_to, const vector &x_from) retorna ln q(x′∣x), onde x′ é o novo estado proposto e x é o estado atual.

Essa classe é usada apenas com distribuições de proposta assimétricas. Nesses casos, ln q(x′∣x) é necessário para calcular a correção de Hastings no critério de aceitação. Para distribuições simétricas, essa classe não é utilizada.

Classe MHSampler

O método MHSample executa o algoritmo de Metropolis ou Metropolis-Hastings e recebe os seguintes parâmetros:

• start - estado inicial x0, representado por um vetor,
• samples - matriz para armazenar as amostras obtidas,
• accept_rate - variável para armazenar a taxa de aceitação das propostas,
• log_pdf - ponteiro para um objeto da classe LogPDF, usado para calcular ln p(x),
• prop_rnd - ponteiro para um objeto da classe PropRND, usado para gerar os candidatos,
• log_prop_pdf - ponteiro para um objeto da classe LogPropPDF, ou NULL no caso simétrico,
• params - estrutura MH_Params que contém os seguintes parâmetros:
• nsamples - número de amostras armazenadas,
• burnin - período de aquecimento, isto é, o número de iterações iniciais descartadas para que a cadeia alcance a distribuição estacionária,
• thin - intervalo de afinamento (thinning), no qual apenas uma amostra a cada i é armazenada para reduzir a correlação,
• symmetric — sinalizador que indica se a distribuição de proposta é simétrica, isto é, se deve ser usada a fórmula simplificada (Metropolis) ou a fórmula completa (MH).

Em cada iteração, partindo do estado atual x0, são executadas as seguintes etapas:

1. Geração do candidato: x' ∼ q(x'∣x0) com o objeto prop_rnd.
2. Critério de aceitação: calcula-se o logaritmo da razão de aceitação r:

• Para o caso simétrico (Metropolis):

r = lnp*(x') − lnp*(x0)

• Para o caso assimétrico (Metropolis-Hastings):

r = [ lnp*(x') + lnq(x0∣x') ] − [ lnp*(x0) + lnq(x'∣x0) ]

3. Aceitação/Rejeição: gera-se um número aleatório uniforme u ∼ U(0,1). O passo é aceito se a condição ln u ≤ min (r, 0) for satisfeita. Nesse caso, o estado atual é atualizado para x0 = x'; caso contrário, x0 permanece inalterado.

4. Coleta das amostras: após o período de aquecimento (burnin) e considerando o intervalo de afinamento (thin), o estado atual x0 é gravado na matriz samples.

Ao final do laço, a taxa de aceitação (accept_rate) é sempre calculada. Ela é uma das principais métricas de diagnóstico para avaliar a eficiência da distribuição de proposta escolhida.

#include <Math\Stat\Uniform.mqh>

//+------------------------------------------------------------------+
//| Abstract class for the target log-density P(x)                   |
//+------------------------------------------------------------------+
class LogPDF
  {
public:
   virtual double    LogPdf(const vector &x) = 0;
  };

//+------------------------------------------------------------------+
//| Abstract class for the proposal log-density q(x'|x)              |
//| (For the asymmetric case)                                        |
//+------------------------------------------------------------------+
class LogPropPDF
  {
public:
   // LogPropPdf(x_to, x_from) returns q(x_to | x_from)
   virtual double    LogPropPdf(const vector &x_to, const vector &x_from) = 0;
  };

//+------------------------------------------------------------------+
//| Abstract class for a candidate generator q(x'|x)                 |
//+------------------------------------------------------------------+
class PropRND
  {
public:
   virtual vector    PropRnd(const vector &x) = 0;
  };

//+------------------------------------------------------------------+
//| Parameter structure for MHSampler                                |
//+------------------------------------------------------------------+
struct MH_Params
  {
   int               nsamples;     // number of samples
   int               burnin;       // burn-in period
   int               thin;         // thinning
   bool              symmetric;   // A symmetric proposal distribution q(x', x) = q(x, x')
  };

//+------------------------------------------------------------------+
//| MHSampler class implementing the Metropolis-Hastings algorithm   |
//+------------------------------------------------------------------+
class MHSampler
  {
public:
   //+------------------------------------------------------------------+
   //| Main MH sampling function                                        |
   //+------------------------------------------------------------------+
   bool              MHSample(
      const vector &start,
      matrix &samples,
      double &accept_rate,
      LogPDF *log_pdf,
      PropRND *prop_rnd,
      LogPropPDF *log_prop_pdf, // May be NULL if symmetric = true
      const MH_Params &params
   )
     {
      // --- 1. Validation and initialization ---
      if(params.nsamples <= 0 || params.thin <= 0)
        {
         Print("Error: nsamples and thin should be > 0");
         return false;
        }
      if(!params.symmetric && log_prop_pdf == NULL)
        {
         Print("Error: log_prop_pdf is required for asymmetric proposal distribution.");
         return false;
        }

      int dim = (int)start.Size();
      int total_steps = params.nsamples * params.thin + params.burnin;
      vector x0 = start; // Current value
      double accepted_count = 0;
      int sample_idx = 0;

      // Initialize the sample matrix
      samples.Resize(params.nsamples, dim);

      // --- 2. Main Metropolis-Hastings loop ---
      for(int i = 1 - params.burnin; i <= params.nsamples * params.thin; i++)
        {

         // 2.1. Candidate generation: x' ~ q(x'| x0)
         vector x_new = prop_rnd.PropRnd(x0);

         // 2.2. Computing log probabilities
         double log_pdf_new = log_pdf.LogPdf(x_new);
         double log_pdf_x0 = log_pdf.LogPdf(x0);

         // 2.3. Computing the log acceptance ratio A
         double r;
         if(params.symmetric)
           {
            // Random Walk MH (Symmetric q(x,x') = q(x',x)): r = log(P(x')/P(x0))
            r = log_pdf_new - log_pdf_x0;
           }
         else
           {
            // General MH (Asymmetric q): r = log( P(x')q(x0|x') / P(x0)q(x'|x0) )
            // q(x_new | x0) - Forward transition probability
            double log_prop_x0_to_xnew = log_prop_pdf.LogPropPdf(x_new, x0); // q(x' | x0)
            // q(x0 | x_new) - Reverse transition probability
            double log_prop_xnew_to_x0 = log_prop_pdf.LogPropPdf(x0, x_new); // q(x0 | x')

            r = (log_prop_xnew_to_x0 + log_pdf_new) - (log_prop_x0_to_xnew + log_pdf_x0);
           }

         // 2.4. Acceptance/Rejection
         int err;
         double U = MathLog(MathRandomUniform(0.0, 1.0,err));
         //A step is accepted if log(U) ≤ min(0, r), which is equivalent to
         //the classical condition U ≤ A, where the acceptance probability is: A = min(1, exp(r))
         if(U <= MathMin(0.0, r))
           {
            // Acceptance: x0 = x'
            x0 = x_new;
            accepted_count++;
           }
         // Rejection: x0 remains unchanged

         // 2.5. Sampling (accounting for burn-in and thinning)
         if(i > 0 && i % params.thin == 0)
           {
            if(sample_idx < params.nsamples)
              {
               samples.Row(x0, sample_idx);
               sample_idx++;
              }
           }
        }

      // 3. Calculating the acceptance rate
      accept_rate = accepted_count / total_steps;

      return true;
     }
  };

Caso unidimensional com distribuição de proposta simétrica (Random Walk Metropolis)

Exemplo nº 1

Este exemplo demonstra a implementação do algoritmo de Metropolis com caminhada aleatória (Random Walk Metropolis, RWM), no qual a distribuição de proposta é simétrica.

Como distribuição-alvo, foi escolhida a distribuição normal padrão unidimensional N(0,1). Os candidatos são gerados com uma distribuição de proposta simétrica:

x' = x + u, onde u ∼ Uniform(−δ, δ)

São utilizados os seguintes parâmetros:

• ponto inicial x0 = 1,
• período de aquecimento burnin = 500,
• intervalo de afinamento (thin) = 10,
• δ = 4, que determina a largura da distribuição uniforme U(−δ, +δ).

Se o passo da cadeia δ for grande demais, a taxa de aceitação será baixa e a cadeia ficará frequentemente "presa" devido à rejeição dos candidatos. Se o passo δ for pequeno demais, a taxa de aceitação será alta, mas a cadeia avançará lentamente, aumentando a correlação entre as amostras.

O ajuste do parâmetro δ permite equilibrar a velocidade de exploração do espaço e a probabilidade de aceitação dos candidatos. O usuário precisa escolher um valor adequado de δ para obter uma amostra independente, ou quase independente.

#include <Graphics\Graphic.mqh>
#include <MCMC\MH.mqh>

//+------------------------------------------------------------------+
//| LogPDF implementation for N(0, 1)                                |
//+------------------------------------------------------------------+
class NormalLogPDF : public LogPDF
  {
public:
   virtual double    LogPdf(const vector &x) override
     {
      // Log(N(x | 0, 1)) = -0.5*x^2 - 0.5*log(2*pi)
      return -0.5 * x[0] * x[0] - 0.5 * MathLog(2.0 * M_PI);
     }
  };

//+------------------------------------------------------------------+
//| PropRND implementation (Random Walk)                             |
//+------------------------------------------------------------------+
class RndWalkPropRND : public PropRND
  {
public:
   virtual vector    PropRnd(const vector &x) override
     {
      double delta = 4; // Random Walk range [-4, 4]
      double x0 = x[0];
      int err;
      double u = MathRandomUniform(-delta, delta,err) ; // Uniform(-d, d)
      vector x_new(1);
      x_new[0] = x0 + u; // x' = x + u
      return x_new;
     }
  };

//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
// Parameters
   vector start_point(1);
   start_point[0] = 1.0;

// Object initialization
   NormalLogPDF target_pdf;
   RndWalkPropRND proposal_rnd;

   MH_Params params;
   params.nsamples = 5000;
   params.burnin = 500;
   params.thin = 10;
   params.symmetric = true; // Random Walk — symmetric proposal distribution

   MHSampler sampler;
   matrix samples;
   double accept_rate;

   Print("Launch Random Walk MH for N(0, 1)...");

// Running MH
   if(sampler.MHSample(start_point, samples, accept_rate, &target_pdf, &proposal_rnd,
                       NULL, // log_prop_pdf is not needed because symmetric = true
                       params))
     {
      Print("=========================================");
      PrintFormat("Acceptance probability: %.2f%%", accept_rate * 100.0);

      vector mean = samples.Mean(0);
      PrintFormat("Mean (0 expected): %.4f", mean[0]);

      Print("First 5 samples:");
      for(int i = 0; i < 5; i++)
        {
         PrintFormat("Sample %d: %.4f", i + 1, samples[i][0]);
        }
     }

// Saving samples to CSV
   MatrixToCSV("MH_RW/MH_samples.csv", samples);

   plotTrace(samples);
  }

//+------------------------------------------------------------------+
//|Sample plot                                                       |
//+------------------------------------------------------------------+
void plotTrace(matrix &smp)
  {
   ChartSetInteger(0, CHART_SHOW, false);
   CGraphic graphic;
   ulong width = ChartGetInteger(0, CHART_WIDTH_IN_PIXELS);
   ulong height = ChartGetInteger(0, CHART_HEIGHT_IN_PIXELS);
   graphic.Create(0, "MH_RW_MCMC", 0, 0, 0, int(width), int(height));
   graphic.BackgroundMain("Trace Plot MH_RW");
   graphic.BackgroundMainSize(16);

   vector v = smp.Col(0);
   double x[];
   v.Swap(x);
   graphic.CurveAdd(x, CURVE_LINES, "Sample");
   graphic.CurvePlotAll();
   graphic.Update();
   Sleep(10 * 1000);
   ChartSetInteger(0, CHART_SHOW, true);
   graphic.Destroy();
   ChartRedraw(0);
  }
//+------------------------------------------------------------------+

Depois de obter as amostras, elas são analisadas com Python. As principais ferramentas de diagnóstico são:

• Gráfico de trajetória (Trace Plot): gráfico dos valores x(t) ao longo das iterações, que permite avaliar a convergência da cadeia e seu grau de mistura.
• Histograma: comparação da distribuição empírica das amostras com a densidade teórica da distribuição-alvo N(0,1).
• Função de autocorrelação (ACF): mostra o nível de correlação entre amostras sucessivas. No RWM, a autocorrelação costuma ser alta devido à própria natureza da caminhada aleatória, por isso é necessário ajustar o afinamento (thinning).

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf

# File path
file_path = "C:/Program Files/MetaTrader 5/MQL5/Files/MH_RW/MH_samples.csv"

# Load data
data = pd.read_csv(file_path, header=None)
samples = data.iloc[:, 0].values

# Statistics
print("=========================================")
print("MCMC Samples Analysis (MH_RW, N(0, 1)):")
print(f"Number of samples: {len(samples)}")
print(f"Mean: {np.mean(samples):.4f} (expected 0.0000)")
print(f"Std. deviation: {np.std(samples):.4f} (expected 1.0000)")
print("=========================================")

# Plots
fig, axes = plt.subplots(3, 1, figsize=(8, 8))

# Trace Plot
axes[0].plot(samples)
axes[0].set_title('Trace Plot')
axes[0].set_xlabel('Iteration')
axes[0].set_ylabel('X')

# Histogram and PDF
x_range = np.linspace(samples.min(), samples.max(), 100)
pdf_theoretical = (1 / np.sqrt(2 * np.pi)) * np.exp(-0.5 * x_range**2)
axes[1].hist(samples, bins=50, density=True, alpha=0.6, label='Samples')
axes[1].plot(x_range, pdf_theoretical, 'r-', label='PDF N(0, 1)')
axes[1].set_title('Histogram and PDF')
axes[1].legend()

# ACF
plot_acf(samples, lags=50, ax=axes[2], title='ACF')
axes[2].set_xlabel('Lags')
axes[2].set_ylabel('Correlation')

plt.tight_layout()
plt.show()

Fig. 3 Gráfico de trajetória e PDF de N(0,1)

Caso multidimensional com distribuição de proposta simétrica (RWM)

Exemplo nº 2

Este exemplo demonstra a aplicação do algoritmo RWM em um espaço multidimensional com variáveis correlacionadas.

#include <Math\Stat\Normal.mqh>
#include <Graphics\Graphic.mqh>
#include <MCMC\MH.mqh>

//+------------------------------------------------------------------+
//|Class for the bivariate normal distribution                       |
//+------------------------------------------------------------------+
class BivariateNormalLogPDF : public LogPDF
  {
public:
   virtual double    LogPdf(const vector &x) override
     {
      vector mu(2);
      mu[0] = 0;
      mu[1] = 0; // Mean
      matrix Sigma(2, 2);   // Covariance matrix
      Sigma[0][0] = 1;
      Sigma[0][1] = 0.8;
      Sigma[1][0] = 0.8;
      Sigma[1][1] = 1;
      matrix Sigma_inv(2, 2);
      Sigma_inv = Sigma.Inv();
      double det = Sigma.Det();
      vector diff = x - mu;
      double quad_form = diff @ Sigma_inv @ diff;
      return -0.5 * (quad_form + MathLog(det) + 2 * MathLog(2 * M_PI));
     }
  };

//
//+------------------------------------------------------------------+
//|Symmetric proposal distribution (Random Walk)                     |
//+------------------------------------------------------------------+
class MultiRndWalkPropRND : public PropRND
  {
public:
   virtual vector    PropRnd(const vector &x) override
     {
      double sigma = 1;
      vector delta(x.Size());
      int err;
      for(ulong i = 0; i < x.Size(); i++)
        {
         delta[i] = MathRandomNormal(0, sigma, err);
        }
      return x + delta;
     }
  };

//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
   vector start_point(2);
   start_point[0] = 1.0;
   start_point[1] = 1.0;
   BivariateNormalLogPDF target_pdf;
   MultiRndWalkPropRND proposal_rnd;

   MH_Params params;
   params.nsamples = 5000;
   params.burnin = 500;
   params.thin = 10;
   params.symmetric = true;

   MHSampler sampler;
   matrix samples;
   double accept_rate;

   Print("Launching Random Walk MH for bivariate N([0,0], [[1,0.8],[0.8,1]])...");
   if(sampler.MHSample(start_point, samples, accept_rate, &target_pdf, &proposal_rnd, NULL, params))
     {
      PrintFormat("Acceptance probability: %.2f%%", accept_rate * 100.0);
      vector mean = samples.Mean(0);
      PrintFormat("Mean (expected [0,0]): [%.4f, %.4f]", mean[0], mean[1]);
      Print("First 5 samples:");
      for(int i = 0; i < 5; i++)
        {
         PrintFormat("Sample %d: [%.4f, %.4f]", i + 1, samples[i][0], samples[i][1]);
        }

      MatrixToCSV("MH_RW/MH_bivariate_samples.csv", samples);

      plotBivariateTrace(samples);
     }
  }

//+------------------------------------------------------------------+
//| Plot of bivariate samples                                        |
//+------------------------------------------------------------------+
void plotBivariateTrace(matrix &smp)
  {
   ChartSetInteger(0, CHART_SHOW, false);
   CGraphic graphic;
   ulong width = ChartGetInteger(0, CHART_WIDTH_IN_PIXELS);
   ulong height = ChartGetInteger(0, CHART_HEIGHT_IN_PIXELS);
   graphic.Create(0, "MH_Bivariate_MCMC", 0, 0, 0, int(width), int(height));
   graphic.BackgroundMain("Bivariate Normal Distribution");
   graphic.BackgroundMainSize(16);

   vector x1 = smp.Col(0);
   vector y1 = smp.Col(1);
   double x[],y[];
   x1.Swap(x);
   y1.Swap(y);
   graphic.CurveAdd(x, y, CURVE_POINTS, "Samples");
   graphic.CurvePlotAll();
   graphic.Update();
   Sleep(10 * 1000);
   ChartSetInteger(0, CHART_SHOW, true);
   graphic.Destroy();
   ChartRedraw(0);
  }

Como distribuição-alvo, para simplificar, escolheremos uma distribuição normal bidimensional N(μ, Σ) com os seguintes parâmetros:

• vetor de médias μ=[0,0],
• matriz de covariância Σ= { {1, 0.8}, {0.8, 1} }

Essa matriz reflete uma forte correlação entre as duas variáveis, o que torna a amostragem mais difícil.

Para calcular o logaritmo da densidade da distribuição-alvo, criaremos a classe BivariateNormalLogPDF, que utiliza a fórmula padrão da distribuição normal multivariada:

ln p(x)= −0.5(x − μ)^T Σ^−1(x − μ) − 0.5ln∣Σ∣ − ln(2π)

Como usamos o algoritmo RWM simétrico, a distribuição de proposta também deve ser simétrica. Ela é implementada na classe MultiRndWalkPropRND. O candidato x′ é gerado a partir de uma distribuição normal isotrópica da seguinte forma:

x' = x + δ

onde:

• δ ∼ N(0,σ^2*I),
• σ^2=1,
• I é a matriz identidade.

Essa escolha corresponde ao comportamento de uma caminhada aleatória centrada no estado atual x.

A amostragem é realizada com os seguintes parâmetros:

• nsamples=5000, burnin=500, thin=10.
• Ponto inicial: [1.0,1.0].

Para avaliar a qualidade das amostras, utiliza-se um script em Python que gera:

• Gráficos de trajetória das amostras: gráficos dos valores x1(t) e x2(t) ao longo das iterações para avaliar a convergência e a mistura da cadeia,
• Função de autocorrelação (ACF): mostra a correlação entre amostras sucessivas de cada variável (x1,x2),
• Gráfico de dispersão: nuvem de amostras com elipses de confiança de 95% e 99% sobrepostas, correspondentes à distribuição teórica N(μ, Σ).

Fig. 4 Gráfico de amostras bidimensionais e elipses de confiança

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf
from matplotlib.patches import Ellipse
from scipy.linalg import eigh

# File path
file_path = "C:/Program Files/MetaTrader 5/MQL5/Files/MH_RW/MH_bivariate_samples.csv"

# Load data
data = pd.read_csv(file_path, header=None)
samples_x1 = data.iloc[:, 0].values
samples_x2 = data.iloc[:, 1].values

# Statistics
print("=========================================")
print("MCMC Samples Analysis (MH_RW, Bivariate):")
print(f"Number of samples: {len(samples_x1)}")
print(f"Mean X1: {np.mean(samples_x1):.4f} (expected 0.0000)")
print(f"Mean X2: {np.mean(samples_x2):.4f} (expected 0.0000)")
print(f"Std. deviation X1: {np.std(samples_x1):.4f} (expected 1.0000)")
print(f"Std. deviation X2: {np.std(samples_x2):.4f} (expected 1.0000)")
print(f"Covariance: {np.cov(samples_x1, samples_x2)[0, 1]:.4f} (expected 0.8000)")
print("=========================================")

# Ellipse plotting function
def plot_ellipse(mean, cov, ax, n_std=1.96, **kwargs):
    vals, vecs = eigh(cov)
    order = vals.argsort()[::-1]
    vals, vecs = vals[order], vecs[:, order]
    theta = np.degrees(np.arctan2(*vecs[:, 0][::-1]))
    width, height = 2 * n_std * np.sqrt(vals)
    ellipse = Ellipse(xy=mean, width=width, height=height, angle=theta, **kwargs)
    ax.add_patch(ellipse)

# Trace Plots and ACF
fig1, axes = plt.subplots(2, 2, figsize=(10, 8))
fig1.suptitle('MCMC Samples Analysis (Trace and ACF)', fontsize=16)

# Trace Plot for X1
axes[0, 0].plot(samples_x1)
axes[0, 0].set_title('Trace Plot (X1)')
axes[0, 0].set_xlabel('Iteration')
axes[0, 0].set_ylabel('X1')

# Trace Plot for X2
axes[0, 1].plot(samples_x2)
axes[0, 1].set_title('Trace Plot (X2)')
axes[0, 1].set_xlabel('Iteration')
axes[0, 1].set_ylabel('X2')

# ACF for X1
plot_acf(samples_x1, lags=50, ax=axes[1, 0], title='ACF (X1)')
axes[1, 0].set_xlabel('Lags')
axes[1, 0].set_ylabel('Correlation')

# ACF for X2
plot_acf(samples_x2, lags=50, ax=axes[1, 1], title='ACF (X2)')
axes[1, 1].set_xlabel('Lags')
axes[1, 1].set_ylabel('Correlation')

plt.tight_layout(rect=[0, 0.03, 1, 0.95])

# Scatter Plot with Ellipses
fig2, ax = plt.subplots(figsize=(8, 8))
ax.scatter(samples_x1, samples_x2, s=10, alpha=0.5, label='Samples')
plot_ellipse([0, 0], np.array([[1, 0.8], [0.8, 1]]), ax, n_std=1.96, edgecolor='red', facecolor='none', label='95% Ellipse')
plot_ellipse([0, 0], np.array([[1, 0.8], [0.8, 1]]), ax, n_std=3.034, edgecolor='green', facecolor='none', label='99% Ellipse')
ax.set_title('Bivariate Scatter Plot with Ellipses')
ax.set_xlabel('X1')
ax.set_ylabel('X2')
ax.legend()
ax.grid(True)
ax.set_aspect('equal')

plt.show()

Distribuição de proposta assimétrica

Exemplo nº 3

Este exemplo demonstra a necessidade de usar o algoritmo completo de Metropolis-Hastings. Isso ocorre porque a distribuição-alvo possui um domínio limitado aos valores positivos, o que exige o uso de uma distribuição de proposta q assimétrica.

Como distribuição-alvo, escolheremos uma distribuição gama unidimensional Gamma(α=3,β=2). Essa distribuição é definida apenas para x>0.

#include <Math\Stat\Math.mqh>
#include <Math\Stat\Gamma.mqh>
#include <Math\Stat\Lognormal.mqh>
#include <Graphics\Graphic.mqh>
#include <MCMC\MH.mqh>

// Class for the target gamma distribution
class GammaLogPDF : public LogPDF
  {
public:
   virtual double    LogPdf(const vector &x) override
     {
      double alpha = 3.0; // Shape parameter
      double beta = 2.0;  // Scale parameter
      return alpha * MathLog(beta) - MathGammaLog(alpha) + (alpha - 1) * MathLog(x[0]) - beta * x[0];
     }
  };

//+--------------------------------------------------------------------------+
//| Class for an asymmetric proposal distribution (lognormal)                |
//+--------------------------------------------------------------------------+
class LogNormalPropRND : public PropRND
  {
public:
   virtual vector    PropRnd(const vector &x) override
     {
      double sigma = 2;
      double mu = MathLog(x[0]);
      vector x_new(1);
      int err;
      x_new[0] = MathRandomLognormal(mu, sigma, err);
      return x_new;
     }
  };

//+------------------------------------------------------------------+
//|Proposal log-density class                                        |
//+------------------------------------------------------------------+
class LogNormalLogPropPDF : public LogPropPDF
  {
public:
   virtual double    LogPropPdf(const vector &x, const vector &y) override
     {
      double sigma = 2;
      double mu = MathLog(y[0]);
      int err;
      double log_pdf = MathProbabilityDensityLognormal(x[0], mu, sigma, true, err);
      return log_pdf;
     }
  };

//+------------------------------------------------------------------+
//| Sample plot                                                      |
//+------------------------------------------------------------------+
void plotTrace(matrix &smp)
  {
   ChartSetInteger(0, CHART_SHOW, false);
   CGraphic graphic;
   ulong width = ChartGetInteger(0, CHART_WIDTH_IN_PIXELS);
   ulong height = ChartGetInteger(0, CHART_HEIGHT_IN_PIXELS);
   graphic.Create(0, "MH_Gamma_MCMC", 0, 0, 0, int(width), int(height));
   graphic.BackgroundMain("Gamma Distribution Samples");
   graphic.BackgroundMainSize(16);

   vector x = smp.Col(0);
   double x_arr[], y_arr[];
   x.Swap(x_arr);
   graphic.CurveAdd(x_arr, CURVE_POINTS, "Samples");
   graphic.CurvePlotAll();
   graphic.Update();
   Sleep(10 * 1000);
   ChartSetInteger(0, CHART_SHOW, true);
   graphic.Destroy();
   ChartRedraw(0);
  }

//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
   vector start_point(1);
   start_point[0] = 1.5;
   GammaLogPDF target_pdf;
   LogNormalPropRND proposal_rnd;
   LogNormalLogPropPDF proposal_pdf;

   MH_Params params;
   params.nsamples = 5000;
   params.burnin = 500;
   params.thin = 10;
   params.symmetric = false;

   MHSampler sampler;
   matrix samples;
   double accept_rate;

   Print("Launching MH with asymmetric proposal for Gamma(3, 2)...");
   if(sampler.MHSample(start_point, samples, accept_rate, &target_pdf, &proposal_rnd, &proposal_pdf, params))
     {
      PrintFormat("Acceptance probability: %.2f%%", accept_rate * 100.0);
      PrintFormat("Mean (expected %.4f): %.4f", 3.0/2.0, samples.Mean(0)[0]);
      Print("First 5 samples:");
      for(int i = 0; i < 5; i++)
        {
         PrintFormat("Sample %d: %.4f", i + 1, samples[i][0]);
        }
      MatrixToCSV("MH_RW/MH_gamma_samples.csv", samples);
      plotTrace(samples);
     }
  }

O logaritmo da densidade, implementado na classe GammaLogPDF, herdada da classe abstrata LogPDF, é calculado como:

ln p(x) = αlnβ − lnΓ(α) + (α − 1)lnx − βx

Como distribuição de proposta q(x′∣x), usaremos a distribuição lognormal LogNormal(ln x, σ²). Essa distribuição é naturalmente adequada para gerar candidatos positivos x′ a partir do estado positivo atual x, pois centraliza o logaritmo da nova proposta no logaritmo do estado atual. No entanto, ela é assimétrica, portanto precisamos implementar duas classes:

• LogNormalPropRND para gerar os candidatos x′,
• LogNormalLogPropPDF para calcular o logaritmo da densidade da proposta.

A amostragem é realizada com os seguintes parâmetros:

• nsamples=5000, burnin=500, thin=10, symmetric=false,
• Ponto inicial: x0=1.5

Para avaliar o funcionamento do algoritmo neste exemplo didático, comparamos os resultados empíricos com as características estatísticas esperadas da distribuição gama:

• média esperada: E[X] = α/β = 1.5,
• variância esperada: Var[X] = α/β^2 =0.75

Em tarefas práticas reais, quando as características exatas de p(x) são desconhecidas, a qualidade da cadeia é diagnosticada, como de costume, por meio de gráficos de trajetória, taxa de aceitação e função de autocorrelação (ACF). Os resultados do script MQL5 demonstram uma convergência bem-sucedida: a taxa de aceitação permanece na faixa ideal (≈30%−40%), e a média empírica fica próxima do valor esperado.

Fig. 5 Gráfico de trajetória e PDF

Metropolis-Hastings Independente

Exemplo nº 4

O Metropolis-Hastings Independente (IMH) utiliza uma distribuição de proposta independente q(x'∣x)=q(x'), que não depende do estado atual x. Essa é a principal característica do IMH em comparação com o Random Walk MH.

Graças à independência, q(x'∣x)=q(x') e q(x∣x')=q(x), o logaritmo da razão de aceitação é dado por:

r = ln(p(x')q(x) / p(x)q(x')) = [ lnp(x′) + lnq(x) ] − [ lnp(x) + lnq(x′) ]

O único requisito do algoritmo IMH é que q(x′) cubra o suporte da distribuição-alvo p(x).

Como distribuição-alvo, foi escolhida a distribuição beta Beta(α=2, β=5), definida no intervalo unitário x∈[0,1]. O logaritmo da densidade, implementado na classe BetaLogPDF, é calculado usando o logaritmo da função beta para garantir estabilidade numérica:

ln⁡p*(x) = (α−1)ln⁡x + (β−1)ln⁡(1−x) − ln⁡B(α, β)

Como distribuição de proposta q(x′), escolheremos a distribuição uniforme Uniform(0,1). Essa escolha corresponde perfeitamente ao domínio da distribuição beta-alvo (x∈[0,1]). Ela é implementada por meio de duas classes:

• UniformPropRND: para gerar candidatos x′ ∼ Uniform(0,1).
• UniformLogPropPDF: para calcular lnq(x′).

Uma simplificação importante: como a densidade da distribuição uniforme é q(x)=1 e q(x′)=1 para x∈[0,1], a razão q(x)/q(x′) é igual a 1, e seu logaritmo é igual a 0. Assim, neste caso particular do IMH, a razão de aceitação se reduz à forma de Metropolis:

r = ln p*(x′) − ln p*(x)

A amostragem foi realizada com os seguintes parâmetros:

• ponto inicial: x0 = 0.3
• nsamples=5000, burnin=500, thin=10, symmetric=false, embora, neste caso, isso não faça diferença, pois lnq(x)−lnq(x′)=0.

A principal vantagem do IMH é sua capacidade de passar rapidamente para qualquer ponto do espaço de estados, o que, quando q(x′) é bem escolhida, resulta em baixa autocorrelação em comparação com o RWM.

Para avaliar a qualidade da amostragem neste exemplo didático, comparamos a média empírica obtida das amostras com o valor médio esperado conhecido: E[X]= α + βα ≈ 0.2857.

Os resultados empíricos confirmam que a amostragem da distribuição-alvo Beta(2,5) foi bem-sucedida.

Fig. 6 Gráfico de trajetória e PDF

Escolha do algoritmo MH em MCMC: RWM vs. IMH

A escolha entre Random Walk Metropolis (RWM) e Metropolis-Hastings Independente (IMH) depende do nível de conhecimento prévio sobre a distribuição-alvo p(x).

O RWM é a opção mais versátil e tem maior relevância prática.

O candidato x′ é gerado como um "passeio aleatório" em torno do estado atual x. Essa abordagem é usada quando a forma da distribuição-alvo é pouco conhecida ou quando o espaço tem alta dimensionalidade e estrutura complexa, por exemplo, quando é multimodal. Ela exige um ajuste cuidadoso da escala do passo (δ ou σ²) para alcançar uma taxa de aceitação adequada e evitar uma mistura lenta da cadeia, isto é, alta autocorrelação.

O IMH é mais eficiente para gerar amostras independentes, mas exige conhecimento prévio sobre a distribuição-alvo.

O candidato x′ é gerado a partir de uma distribuição de proposta fixa q(x′), independente de x. Essa abordagem é usada quando é possível escolher q(x′) de modo que ela aproxime bem a distribuição-alvo p(x). Quando q(x′) é bem escolhida, obtém-se uma mistura muito rápida da cadeia, com baixa autocorrelação. No entanto, se q(x′) não cobrir adequadamente as caudas de p(x), a cadeia poderá ficar presa e produzir resultados incorretos.

Considerações finais

O algoritmo de Metropolis-Hastings é um método de Monte Carlo via cadeias de Markov (MCMC) que permite gerar amostras de distribuições complexas para as quais a amostragem direta não é possível. Neste artigo, analisamos tanto os fundamentos teóricos do algoritmo quanto sua implementação em MQL5 na forma da classe MHSampler. Essa classe permite ao usuário adaptar facilmente o algoritmo a diferentes variantes de MCMC, incluindo Random Walk Metropolis e Metropolis-Hastings Independente, bastando implementar a lógica da distribuição-alvo (LogPDF) e da geração de propostas (PropRND, LogPropPDF).

Para demonstrar de forma clara o funcionamento do algoritmo, utilizamos distribuições unidimensionais e bidimensionais como distribuições-alvo. A qualidade das amostras geradas foi avaliada com ferramentas padrão de diagnóstico, como gráficos de trajetória, histogramas e funções de autocorrelação construídos em Python.

A taxa de aceitação das amostras (acceptance rate) serviu como principal indicador da eficiência da cadeia e foi usada para ajustar a escala da distribuição de proposta q(x′∣x), definida, por exemplo, pelo desvio-padrão no caso de uma distribuição normal ou pela largura do intervalo no caso de uma distribuição uniforme. A escolha correta dessa escala é decisiva: uma distribuição estreita demais torna a exploração do espaço de estados mais lenta e aumenta a autocorrelação entre as amostras, enquanto uma distribuição ampla demais reduz a taxa de aceitação e, consequentemente, a eficiência computacional.

Embora o algoritmo seja bastante versátil, sua eficiência diminui em problemas multidimensionais e quando há forte correlação entre as variáveis. Nesses casos, métodos mais especializados, como Monte Carlo híbrido (HMC) ou amostragem por fatias (slice sampling), são preferíveis, pois aproveitam melhor a geometria da distribuição-alvo e exigem menos ajuste manual. Ainda assim, devido à sua simplicidade, o algoritmo de Metropolis-Hastings continua sendo uma ferramenta fundamental e amplamente utilizada na inferência bayesiana.

Programas utilizados no artigo

#	Nome	Tipo	Descrição
1	MH.mqh	Arquivo de inclusão	Algoritmo de Metropolis-Hastings
2	Exmp1_MH_RW.mq5	Script	Exemplo de amostragem de uma distribuição unidimensional com RWM
3	Exmp1_Plot_MH_RW.py	Script	Diagnóstico em Python
4	Exmp2_MH_RW.mq5	Script	Exemplo de amostragem de uma distribuição bidimensional com o algoritmo RWM
5	Exmp2_Plot_MH_RW.py	Script	Diagnóstico em Python
6	Exmp3_MH_RW.mq5	Script	Exemplo de amostragem de uma distribuição unidimensional com distribuição de proposta assimétrica
7	Exmp3_Plot_MH_RW.py	Script	Diagnóstico em Python
8	Exmp4_IMH.mq5	Script	Exemplo de amostragem de uma distribuição unidimensional com o algoritmo IMH
9	Exmp4_Plot_IMH.py	Script	Diagnóstico em Python