Métodos de amostragem MCMC: algoritmo de amostragem por fatias (Slice Sampling)

Introdução

Os métodos de Monte Carlo via cadeias de Markov (MCMC) são amplamente utilizados para gerar amostras de distribuições multidimensionais complexas, sobretudo em estatística bayesiana. Algoritmos clássicos, como a amostragem de Gibbs ou o algoritmo de Metropolis, muitas vezes exigem um ajuste prévio cuidadoso. No primeiro caso, é necessário obter expressões analíticas para todas as distribuições condicionais completas; no segundo, é preciso ajustar criteriosamente a escala e a forma da distribuição de proposta. Isso dificulta sua aplicação rápida e eficiente na prática cotidiana.

Neste artigo, estudamos o método de amostragem por fatias (slice sampling), uma variante adaptativa dos métodos MCMC. Sua principal vantagem é a capacidade de se ajustar automaticamente às características da distribuição-alvo, eliminando a necessidade de ajustar manualmente parâmetros como a largura inicial do intervalo.

Após descrever em detalhes o algoritmo e sua implementação no ambiente MQL5, vamos testá-lo em modelos bayesianos de regressão linear e logística, comparando os resultados com métodos frequentistas clássicos, que fornecem estimativas pontuais dos parâmetros.

Métodos de amostragem por fatias

Uma das principais dificuldades dos algoritmos MCMC, como o algoritmo de Metropolis, é escolher uma escala de passo adequada. Um passo muito pequeno faz a cadeia se mover mais lentamente, aumenta a autocorrelação entre as amostras e exige um grande número de iterações para obter amostras independentes. Um passo muito grande, por sua vez, leva à rejeição frequente dos pontos propostos, reduzindo a eficiência do algoritmo.

Para resolver esse problema, foi proposta a técnica de amostragem por fatias (Neal, 2003), que requer apenas a capacidade de avaliar a densidade não normalizada f(x). Esses métodos podem ser divididos em dois tipos, de acordo com a abordagem adotada para distribuições multidimensionais:

• amostragem unidimensional por fatias, com atualização alternada das coordenadas,
• amostragem multidimensional por fatias, com atualização direta de todo o vetor de variáveis.

Amostragem unidimensional por fatias

A amostragem unidimensional por fatias é utilizada tanto para distribuições-alvo unidimensionais quanto para a atualização alternada das coordenadas de um vetor multidimensional x = (x1,…, xn), de forma semelhante à amostragem de Gibbs. Vejamos como esse algoritmo funciona.

A atualização do ponto atual x0 para um novo valor x1 envolve três etapas principais:

• (a) Definição da fatia: escolhemos uma variável auxiliar y ∼ Uniform(0, f(x0)) e definimos a fatia horizontal S = {x: f(x) > y}, isto é, o conjunto de pontos x para os quais o valor da densidade é maior que o valor da variável auxiliar y.
• (b) Construção do intervalo: determinamos um intervalo (L, R) de largura w, posicionado aleatoriamente em torno de x0. O intervalo é expandido em incrementos de w pelo procedimento "stepping-out", até que suas duas extremidades fiquem fora da fatia.
• (c) Seleção do novo ponto: o novo ponto x1 é escolhido uniformemente no intervalo (L, R). Se o ponto estiver fora da fatia, o intervalo é reduzido pelo procedimento "shrinkage", e a seleção é repetida até que x1 pertença à fatia.

Na prática, é mais conveniente trabalhar em escala logarítmica:

g(x) = log(f(x)), z = log(y) = g(x0) − e,

onde e ∼ Exp(1) é uma variável aleatória com distribuição exponencial e valor esperado igual a um. Nesse caso, a fatia é S = {x: g(x) > z}.

Os procedimentos "stepping-out" e "shrinkage", propostos por Radford Neal (Neal, 2003), renomado especialista em machine learning e métodos MCMC, garantem uma amostragem eficiente e correta no caso unidimensional. Uma ilustração gráfica desse algoritmo é apresentada na Fig. 1.

Fig. 1. Atualização unidimensional por amostragem por fatias com o uso dos procedimentos de expansão (stepping-out) e contração (shrinkage).

Para visualizar o funcionamento do algoritmo no caso unidimensional, pode-se executar o script a seguir, que demonstra a amostragem de uma distribuição multimodal (Fig. 2).

//+------------------------------------------------------------------+
//|                                                       PlotMM.mq5 |
//|                                                           Eugene |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Eugene"
#property link      "https://www.mql5.com"
#property version   "1.00"
#property script_show_inputs

#include <Math\Stat\Exponential.mqh>
#include <Math\Stat\Normal.mqh> 
#include <Math\Stat\Uniform.mqh> 
#include <Math\Stat\Math.mqh>
#include <Graphics\Graphic.mqh>

input int time_ = 700;     // Animation speed (lower values mean faster)

// Defining the function type for logpdf
typedef double (*LogPdfFunction)(vector &x);

// Checking whether a point lies within the slice
bool inside(vector &x, double th, LogPdfFunction logpdf) {
   return logpdf(x) > th;
}

// target density
double pdf(double x) {
   return MathExp(-x * x / 2.0) * (1.0 + MathPow(MathSin(3.0 * x), 2)) * (1.0 + MathPow(MathCos(5.0 * x), 2));
}

// log(pdf)
double logpdf(vector &x) {
   return MathLog(pdf(x[0]));
}

//+--------------------------------------------------------+
//| Slice sampling                                         |
//+--------------------------------------------------------+
void slicesample(vector &initial_, int nsamples_, vector &width_, int burnin_, int thin_, matrix &rnd, double &neval, LogPdfFunction logpdf) {
   // Input parameter validation
   if (burnin_ < 0 || thin_ <= 0) {
      Print("Error: burnin should be >= 0, thin > 0");
      return;
   }
  
   // Initialization
   int dim = (int)initial_.Size();
   rnd.Resize(nsamples_, dim);
   rnd.Fill(0.0);
   int maxiter = 200;
   vector x0 = initial_;
   neval = nsamples_;

   // Generate exponential and uniform random numbers
   double e[];
   MathRandomExponential(1.0, nsamples_ * thin_ + burnin_, e);
   matrix rw = matrix::Random(nsamples_ * thin_ + burnin_, dim, 0.0, 1.0);
   matrix rd = matrix::Random(nsamples_ * thin_ + burnin_, dim, 0.0, 1.0);

   //---------------------------
   ChartSetInteger(0, CHART_SHOW, false);
   CGraphic graphic;
   ulong w = ChartGetInteger(0, CHART_WIDTH_IN_PIXELS);
   ulong h = ChartGetInteger(0, CHART_HEIGHT_IN_PIXELS);
   graphic.Create(0, "SliceSampling", 0, 0, 0, (int)w, (int)h);
   graphic.BackgroundMain("Slice Sampling");
   graphic.BackgroundMainSize(20);
   graphic.XAxis().Name("x"); 
   graphic.XAxis().NameSize(20); 
   graphic.YAxis().Name("f(x)");
   graphic.YAxis().NameSize(20); 

   // === 1. Target density plot ===
   int steps = 200;
   double x_density[], y_density[];
   MathSequenceByCount(-4,4,200,x_density);
   ArrayResize(y_density, steps);
   for (int i = 0; i < steps; i++) {
      y_density[i] = pdf(x_density[i]);  
   }

   CCurve *density = graphic.CurveAdd(x_density, y_density, clrRed, CURVE_LINES);
   density.LinesWidth(2);
   graphic.CurvePlotAll(); 
   graphic.Update();

   static string prev_slice_name = "";
   static string prev_interval_name = "";
   static string prev_point_name = "";

   // Main loop
   for (int i = 1 - burnin_; i <= nsamples_ * thin_; i++) {   
      // === 2. starting point x0  ===
      double px_axis[] = {x0[0]}, py_axis[] = {0.0};
      CCurve *point_x0_axis = graphic.CurveAdd(px_axis, py_axis, clrRed, CURVE_POINTS);
      point_x0_axis.PointsType(POINT_CIRCLE);
      point_x0_axis.PointsSize(10);
      graphic.CurvePlotAll(); 
      graphic.Update();
      Sleep(time_);

      // === 3. f(x0) ===
      double logf_x0 = logpdf(x0);
      double f_x0    = MathExp(logf_x0);

      double px[] = {x0[0]}, py[] = {f_x0};
      CCurve *point_x0 = graphic.CurveAdd(px, py, clrBlue, CURVE_POINTS);
      point_x0.PointsType(POINT_CIRCLE);
      point_x0.PointsSize(10);
      graphic.CurvePlotAll(); 
      graphic.Update();
      Sleep(time_);

      // === 4. Vertical slice ===
      double x_line[] = {x0[0], x0[0]}, y_line[] = {0.0, f_x0};
      CCurve *vert_line = graphic.CurveAdd(x_line, y_line, clrBlue, CURVE_LINES);
      vert_line.LinesStyle(STYLE_DOT);
      vert_line.LinesWidth(2);
      graphic.CurvePlotAll();
      graphic.Update();
      Sleep(time_);
   
      // === 5. Horizontal slice z = logf(x0) - e ===
      double z = logpdf(x0) - e[i + burnin_ - 1];
      double xz[] = {-4, 4}, yz[] = {MathExp(z), MathExp(z)};  
      string current_slice_name = "slice_" + (string)i;

      if (prev_slice_name != "") {
         graphic.CurveRemoveByName(prev_slice_name);
      }

      CCurve *slice_z = graphic.CurveAdd(xz, yz, clrBlack, CURVE_LINES, current_slice_name);
      slice_z.LinesWidth(2);
      graphic.CurvePlotAll();
      graphic.Update();
      Sleep(time_);

      // Saving the current curve name
      prev_slice_name = current_slice_name;

      // === 6. Interval [xl, xr] ===
      vector r = width_ * rw.Row(i + burnin_ - 1);
      vector xl = x0 - r;
      vector xr = xl + width_;
      int iter = 0;

      double x_slice[] = {xl[0], xr[0]}, y_slice[] = {MathExp(z), MathExp(z)};
      string current_interval_name = "interval_" + (string)i;

      if (prev_interval_name != "") {
         graphic.CurveRemoveByName(prev_interval_name);
      }
      CCurve *horiz_slice = graphic.CurveAdd(x_slice, y_slice, clrOrange, CURVE_LINES,current_interval_name);
      horiz_slice.LinesWidth(4);
      graphic.CurvePlotAll(); 
      graphic.Update();
      Sleep(time_);  
      
      prev_interval_name = current_interval_name;

      //--- step out 
      if (dim == 1) {
         while (inside(xl, z, logpdf) && iter < maxiter) {
            xl -= width_;
            iter++;
            double xs[] = {xl[0], xr[0]}, ys[] = {MathExp(z), MathExp(z)};
            horiz_slice.Update(xs, ys);
            graphic.CurvePlotAll();
            graphic.Update();
            Sleep(time_);
         }
         if (iter >= maxiter) {
            Print("Error: too many iterations in stepping-out (left)");
            return;
         }
         neval += iter;

         iter = 0;
         while (inside(xr, z, logpdf) && iter < maxiter) {
            xr += width_;
            iter++;
            double xs[] = {xl[0], xr[0]}, ys[] = {MathExp(z), MathExp(z)};
            horiz_slice.Update(xs, ys);
            graphic.CurvePlotAll();
            graphic.Update();
            Sleep(time_);
         }
         if (iter >= maxiter) {
            Print("Error: too many iterations in stepping-out (right)");
            return;
         }
         neval += iter;
      }   
      Sleep(time_); 

      // === 8. xp — a point inside the slice ===
      vector xp = rd.Row(i + burnin_ - 1) * (xr - xl) + xl;
      double p_xp[] = {xp[0]}, p_yp[] = {MathExp(z)};
      string current_point_name = "point" + (string)i;

      if (prev_point_name != "") {
         graphic.CurveRemoveByName(prev_point_name);
      }
      CCurve *point_xp = graphic.CurveAdd(p_xp, p_yp, clrMagenta, CURVE_POINTS,current_point_name);
      point_xp.PointsType(POINT_SQUARE);
      point_xp.PointsSize(4);
      point_xp.PointsFill(true);
      graphic.CurvePlotAll();
      graphic.Update();
      Sleep(time_); 
      
      prev_point_name = current_point_name;

      // === 9. Shrinking ===
      iter = 0;
      while (!inside(xp, z, logpdf) && iter < maxiter) {
         for (int d = 0; d < dim; d++) {        
            if (xp[d] > x0[d]) xr[d] = xp[d];
            else xl[d] = xp[d];
            double xs[] = {xl[0], xr[0]}, ys[] = {MathExp(z), MathExp(z)};
            horiz_slice.Update(xs, ys);          
         }
         vector new_rand = vector::Random(dim, 0.0, 1.0);
         xp = new_rand * (xr - xl) + xl;
         double pxp[] = {xp[0]}, pyp[] = {MathExp(z)};
         point_xp.Update(pxp, pyp);
         graphic.CurvePlotAll(); 
         graphic.Update();
         Sleep(time_);
         iter++;
      }
      if (iter >= maxiter) {
         Print("Error: too many iterations in shrinking");
         return;
      }
      neval += iter;
          
      x0 = xp;
      if (i > 0 && i % thin_ == 0) {
         rnd.Row(xp, i / thin_ - 1);
      }
   }

   neval /= (nsamples_ * thin_ + burnin_);

   Sleep(3000);
   graphic.Destroy();
   ChartSetInteger(0, CHART_SHOW, true);
   ChartRedraw(0);
}

//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
   vector start = {1.5};
   int n = 10;
   vector w = {1};
   int burn = 0;
   int thin = 1;
   matrix samples;
   double neval;

   slicesample(start, n, w, burn, thin, samples,neval, logpdf);  
  }
//+------------------------------------------------------------------+

Fig. 2. Slice sampling multimodal

O gráfico mostra o ponto inicial x0, indicado por um círculo azul-claro, e o valor da densidade nesse ponto f(x0), indicado por um círculo vermelho. A fatia vertical resultante é representada por uma linha tracejada vermelha. Ao longo dessa fatia vertical, escolhe-se uniformemente um valor z, que define a fatia horizontal, representada por uma linha preta fina. O intervalo de largura w, posicionado aleatoriamente em torno de x0, é expandido em incrementos de w até que suas extremidades ultrapassem os limites da fatia, como indicado pela linha azul-clara espessa. O novo ponto x1 é escolhido uniformemente nesse intervalo, repetindo-se a seleção até que seja encontrado um ponto dentro da fatia. Os pontos selecionados fora da fatia são usados para reduzir o intervalo. O gráfico mostra justamente um novo ponto candidato x1, indicado por um quadrado vermelho, que é rejeitado pelo algoritmo por estar fora da fatia.

Amostragem multidimensional por fatias

Em vez de atualizar alternadamente cada coordenada do vetor x = (x1,…, xn), a ideia da amostragem por fatias pode ser aplicada diretamente a uma distribuição multidimensional. Nesse caso, o intervalo unidimensional (L, R) é substituído por um hiper-retângulo H = {x: Li < xi < Ri, i = 1,…, n}, em que Li e Ri definem os limites ao longo de cada eixo.

O procedimento para obter o próximo estado x1 = (x1,1,…, x1,n) a partir do estado atual x0 = (x0,1,…, x0,n) é semelhante ao caso unidimensional:

• (a) Escolher y uniformemente em (0, f(x0)), definindo a fatia S = {x: y < f(x)}.
• (b) Determinar um hiper-retângulo H = (L1, R1) × ⋯ × (Ln, Rn) em torno de x0, de preferência contendo a maior parte da fatia.
• (c) Escolher uniformemente um novo ponto x1 na parte da fatia contida em H. Se o ponto estiver fora da fatia, o hiper-retângulo é reduzido pelo procedimento "shrinkage".

Ao contrário do caso unidimensional, o procedimento "stepping-out" é muito complexo em um espaço multidimensional, pois exige verificar os 2^n vértices do hiper-retângulo, o que se torna computacionalmente dispendioso quando n é grande. Por isso, costuma-se utilizar uma abordagem simplificada: o hiper-retângulo é posicionado aleatoriamente em torno de x0 sem executar o procedimento de expansão. A redução do hiper-retângulo é feita de forma independente em cada eixo até que seja encontrado um ponto dentro da fatia. Esse método é mais simples de implementar, porém menos eficiente do que a amostragem unidimensional. Na atualização alternada, cada coordenada é reduzida apenas o necessário, levando em conta as características locais da densidade. No caso multidimensional, o hiper-retângulo é reduzido simultaneamente em todas as dimensões, o que pode causar uma redução excessiva nas direções em que a densidade varia lentamente.

No código apresentado a seguir, o procedimento "stepping-out" é utilizado apenas no caso unidimensional. Para distribuições multidimensionais, empregamos a abordagem simplificada, sem "stepping-out", baseada no posicionamento aleatório do hiper-retângulo e em sua posterior redução.

//+------------------------------------------------------------------+
//|                                                           SS.mqh |
//|                                                           Eugene |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Eugene"
#property link      "https://www.mql5.com"
#property version   "1.00"

#include <Math\Stat\Exponential.mqh>
#include <Math\Stat\Math.mqh>

// Define the function type for logpdf
typedef double (*LogPdfFunction)(vector &x);

// Check whether a point lies within the slice
bool inside(vector &x, double th, LogPdfFunction logpdf) {
   return logpdf(x) > th;
}

//+--------------------------------------------------------+
//| Slice sampling                                         |
//+--------------------------------------------------------+
void slicesample(vector &initial_, int nsamples_, vector &width_, int burnin_, int thin_, matrix &rnd, double &neval, LogPdfFunction logpdf) {
   // Checking input parameters
   if (burnin_ < 0 || thin_ <= 0) {
      Print("Error: burnin should be >= 0, thin > 0");
      return;
   }

   // Initialization
   int dim = (int)initial_.Size();
   rnd.Resize(nsamples_, dim);
   rnd.Fill(0.0);
   int maxiter = 200; // maximum number of iterations per step
   vector x0 = initial_;
   neval = nsamples_; 

   // Generate exponential random numbers
   double e[];
   MathRandomExponential(1.0, nsamples_ * thin_ + burnin_, e);
   // uniform random variables for randomizing the step width w  
   matrix rw = matrix::Random(nsamples_ * thin_ + burnin_, dim, 0.0, 1.0);
   //uniform random variables for selecting a point within an interval 
   matrix rd = matrix::Random(nsamples_ * thin_ + burnin_, dim, 0.0, 1.0); 

   // Main slice sampling loop
   for (int i = 1 - burnin_; i <= nsamples_ * thin_; i++) {
      double z = logpdf(x0) - e[i + burnin_ - 1]; // Horizontal slice S = {x: log f(x) &gt; z}
 
       // Initial interval:
      vector r = width_ * rw.Row(i + burnin_ - 1); 
      vector xl = x0 - r;
      vector xr = xl + width_;
      int iter = 0;

//--- "step out" applies only to one-dimensional distributions -----------
      if (dim == 1) {
         // step out to the left
         while (inside(xl, z, logpdf) && iter < maxiter) {
            xl -= width_;
            iter++;
         }
         if (iter >= maxiter) {
            Print("Error: too many iterations in stepping-out");
            return;
         }
         neval += iter;

         iter = 0;
         // step out to the right
         while (inside(xr, z, logpdf) && iter < maxiter) {
            xr += width_;
            iter++;
         }
         if (iter >= maxiter) {
            Print("Error: too many iterations in stepping-out");
            return;
         }
         neval += iter;
      }     
 //--- Shrinking ---
      vector xp = rd.Row(i + burnin_ - 1) * (xr - xl) + xl;
// Shrink the interval (or hyperrectangle) if the selected point lies outside the slice
      iter = 0;
      while (!inside(xp, z, logpdf) && iter < maxiter) {
         for (int d = 0; d < dim; d++) {        
           if (xp[d] > x0[d]) xr[d] = xp[d]; // If xp[d] &gt; x0[d], we shrink the right boundary
           else xl[d] = xp[d]; // Otherwise, if xp[d] &lt;= x0[d], we shrink the left boundary
         }
         vector new_rand = vector::Random(dim, 0.0, 1.0);
         xp = new_rand * (xr - xl) + xl;
         iter++;
      }
      if (iter >= maxiter) {
         Print("Error: too many iterations in shrinking");
         return;
      }
      neval += iter;

      x0 = xp; // update the current state
      if (i > 0 && i % thin_ == 0) {
         rnd.Row(xp, i / thin_ - 1);
      }
   }
   neval /= (nsamples_ * thin_ + burnin_);
}

Exemplo nº 1. Treinamento de uma regressão linear bayesiana

Neste exemplo, veremos como utilizar o algoritmo slice sampling para gerar amostras da distribuição posterior dos parâmetros de uma regressão linear bayesiana. Para aplicar o teorema de Bayes, definimos inicialmente o modelo dos dados por meio da função de verossimilhança. Para não tornar o modelo excessivamente complexo, vamos supor que os dados y sigam uma distribuição normal em torno da combinação linear dos atributos Xw:

P(y ∣ w,σ2) = N(y | Xw, σ^2 * I)

onde:

• I é a matriz identidade,
• σ2 é a variância do ruído, que, por simplicidade, consideraremos conhecida.

Os parâmetros desconhecidos desse modelo são w. Na abordagem bayesiana, todas as grandezas desconhecidas são tratadas como variáveis aleatórias. Portanto, precisamos definir uma distribuição para essas variáveis. Essa distribuição é chamada de distribuição a priori, pois é especificada antes que o modelo observe os dados. Como distribuição a priori, podemos utilizar uma distribuição normal multivariada com média zero e matriz de covariância Σw.

P(w) = N(w | 0, Σw)

Após observar os dados, atualizamos nossas crenças a priori sobre os parâmetros calculando a distribuição posterior como o produto da verossimilhança pela distribuição a priori:

p(w ∣ X, y) ∝ p(y ∣ Xw, σ2) ⋅ p(w)

No código, a função LogPosterior calcula a distribuição posterior, que passamos ao nosso sampler para gerar as amostras.

Antes de iniciar a amostragem, criaremos dados sintéticos com parâmetros conhecidos w_true = [1.0, 2.0, 2.0], desvio-padrão do ruído σ = 0.5 e uma amostra com 100 observações.

Definiremos os seguintes parâmetros para o sampler:

• nsamples 10000 - número de amostras geradas,
• burnin 1000 - período de aquecimento da cadeia,
• thin 10 - intervalo de thinning,
• width 0.5 - largura inicial do hiper-retângulo.

//+------------------------------------------------------------------+
//|                                                           LR.mq5 |
//|                                                           Eugene |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Eugene"
#property link      "https://www.mql5.com"
#property version   "1.00"

#include <Math\Stat\Exponential.mqh>
#include <Math\Stat\Normal.mqh>
#include <MCMC\SS.mqh>
#include <Math\Stat\Math.mqh>

// Data Storage Structure 
struct BayesianLinearRegression
  {
   matrix            X;           // Feature matrix (n × (d+1))
   vector            y;           // Response data / output data (n)
   double            sigma;       // Noise standard deviation
   vector            mu_w;        // Prior mean for w = [b, w_1, ..., w_d]
   matrix            Sigma_w;     // Prior covariance matrix for w
   matrix            Sigma_w_inv; // Inverse covariance matrix 
   double            log_det_Sigma_w; // Log determinant of Sigma_w 
  };

BayesianLinearRegression model;

//+------------------------------------------------------------------+
//| Log density of a multivariate normal distribution                |
//+------------------------------------------------------------------+
double multivariate_normal_logpdf(vector &x, vector &mu, matrix &Sigma_inv, double log_det, int dim)
  {
   vector diff = x - mu;
   double quad_form = diff @ Sigma_inv @ diff;
   return -0.5 * (quad_form + log_det + dim * MathLog(2 * M_PI));
  }

//+------------------------------------------------------------------+
//| Target density — posterior distribution                          |
//+------------------------------------------------------------------+
double LogPosterior(vector &params, BayesianLinearRegression &mdl)
  {
   int d = (int)mdl.X.Cols(); // Number of columns in X (d+1, including the column of ones)
   int n = (int)mdl.X.Rows(); // Calculate the number of observations
   
  //---- Log-likelihood (Gaussian likelihood) ------------------------------------
// 1.  Prediction vector
   vector y_pred = mdl.X @ params;
// 2.  Error vector
   vector residuals = mdl.y - y_pred;
// 3. Squared norm of the error vector (e^T e)
   double sum_sq_errors = residuals @ residuals; 
// 4. Log-likelihood
   double log_likelihood_quad_term = -0.5 * sum_sq_errors / MathPow(mdl.sigma, 2);
   double log_likelihood_const_term = n * (-MathLog(mdl.sigma) - 0.5 * MathLog(2 * M_PI));
   double log_likelihood = log_likelihood_quad_term + log_likelihood_const_term;
//-----------------------------------------------------------------------------------------------
// Log prior (multivariate normal prior)
   double log_prior = multivariate_normal_logpdf(params, model.mu_w, model.Sigma_w_inv,
                      model.log_det_Sigma_w, d);

// posterior = likelihood * prior
   return log_likelihood + log_prior;
  }

double LogPost(vector &params)
  {
   return LogPosterior(params, model);
  }

//+------------------------------------------------------------------+
//| OLS parameter estimates                                          |
//+------------------------------------------------------------------+
void OLS(matrix &X, vector &y, vector &w_ols)
  {
// OLS: w_ols = (X^T X)^(-1) X^T y
   matrix XtX_inv = (X.Transpose() @ X).Inv();
   w_ols = XtX_inv @ X.Transpose() @ y;
  }

//+------------------------------------------------------------------+
//| Initializing BLR and sampling with SS                            |
//+------------------------------------------------------------------+
void BayesianLinearRegressionSample(matrix &X, vector &y,
                                    double sigma, vector &mu_w, matrix &Sigma_w,
                                    int nsamples, vector &initial, vector &width,
                                    int burnin, int thin, matrix &samples)
  {
   int d = (int)X.Cols(); 

// Precompute the inverse matrix and the log determinant 
   matrix Sigma_w_inv = Sigma_w.Inv();
   double det = Sigma_w.Det();
   double log_det = MathLog(det);

// Model initialization
   model.X = X;
   model.y = y;
   model.sigma = sigma;
   model.mu_w = mu_w;
   model.Sigma_w = Sigma_w;
   model.Sigma_w_inv = Sigma_w_inv; 
   model.log_det_Sigma_w = log_det; 

   // Count the number of logpdf evaluations and measuring execution time
   double neval = 0;
   ulong start_time = GetMicrosecondCount(); 
   slicesample(initial, nsamples, width, burnin, thin, samples, neval, LogPost);
   ulong end_time = GetMicrosecondCount(); 
   double time_ms = (end_time - start_time) / 1000.0; // Time in milliseconds
   Print("Average number of logpdf evaluations per sample: ", neval);
   Print("Slice sampling execution time: ", time_ms, " ms");
  }

//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
// Synthetic data generation
   int n = 100;  // Number of observations
   int d = 2;    // Number of features 
   matrix X = matrix::Random(n, d+1,0,10); // Feature matrix X
   X.Col(vector::Ones(n),0);  
   vector y(n);

   vector true_w(d + 1); // w = [b, w_1, w_2] True parameters
   true_w[0] = 1.0;      
   for(int i = 1; i < d + 1; i++)
     {
      true_w[i] = 2.0;   
     }

// Center the predictors 
   for(int j = 1; j < d + 1; j++)
     {
      vector col = X.Col(j); 
      double mean = col.Mean(); 
      for(int i = 0; i < n; i++)
        {
         X[i,j] -= mean; 
        }
     }

// formation of the y vector
   double sigma = 0.5; // Standard deviation of noise
   int err;
   for(int i = 0; i < n; i++)
     {
      double y_pred =  X.Row(i) @ true_w ;
      y[i] = y_pred + MathRandomNormal(0.0, sigma, err);
     }

//--- Prior parameters
   vector mu_w = vector::Zeros(d + 1); // Prior mean for w = [b, w_1, w_2]
   matrix Sigma_w = matrix::Identity(d + 1,d + 1);
   Sigma_w = 10*Sigma_w; // Diagonal covariance matrix (variance = 10)
   
//--- Sampling parameters
   int nsamples = 10000;
   int burnin = 1000;
   int thin = 10;
   vector initial(d + 1);
   initial.Fill(1.0); // Initial values [b, w_1, w_2]
   vector width(d + 1);
   width.Fill(0.5);   // Step width for slice sampling

//--- Frequentist approach: OLS parameter estimates
   Print("------------   Frequency approach -----------------");
   vector w_ols;
   OLS(X, y, w_ols);
   Print("OLS estimator w: ", w_ols, " (Expected values [1.0 , 2.0, 2.0])");

   vector lower_ci_freq, upper_ci_freq, se_ols;
   double sigma_hat;
   // Confidence intervals
   ConfidenceIntervals(X, y, w_ols, 0.05, lower_ci_freq, upper_ci_freq, se_ols, sigma_hat);
   Print("95% Confidence intervals:");
   for(int i = 0; i < (int)lower_ci_freq.Size(); i++)
     {
      Print("w_", i, ": [", lower_ci_freq[i], ", ", upper_ci_freq[i], "]");
     }
 
   Print("------------   Bayesian approach -----------------");
   matrix samples;
   BayesianLinearRegressionSample(X, y, sigma, mu_w, Sigma_w,
                                  nsamples, initial, width, burnin, thin, samples);
   vector mean_w(d + 1);
   for(int i = 0; i < d + 1; i++)
     {
      mean_w[i] = samples.Col(i).Mean();
     }
   Print("Average Posterior w: ", mean_w, " (Expected values [1.0 , 2.0, 2.0])");
                                  
   vector lower_ci_bayes, upper_ci_bayes;
   // Confidence intervals
   CredibleIntervals(samples, 0.05, lower_ci_bayes, upper_ci_bayes);
   Print("95% Credible intervals:");
   for(int i = 0; i < (int)lower_ci_bayes.Size(); i++)
     {
      Print("w_", i, ": [", lower_ci_bayes[i], ", ", upper_ci_bayes[i], "]");
     }

   MatrixToCSV("SS/LR_samples.csv", samples);
  }

//+------------------------------------------------------------------+
//|Confidence intervals, frequentist approach                        |
//+------------------------------------------------------------------+
void ConfidenceIntervals(matrix &X, vector &y, vector &w_ols, double alpha, vector &lower, vector &upper, vector &se, double &sigma_hat)
  {
   int n = (int)X.Rows();
   int d = (int)X.Cols();

// 1. Compute residuals and sigma^2
   vector y_pred = X @ w_ols;
   vector residuals = y - y_pred;
   double sigma2 = (residuals @ residuals) / (n - d); // Variance estimate
   sigma_hat = MathSqrt(sigma2); // Noise standard deviation estimate

// 2. Covariance matrix
   matrix XtX_inv = (X.Transpose() @ X).Inv();
   matrix cov_matrix = sigma2 * XtX_inv;

// 3. Standard errors
   se = cov_matrix.Diag();
   se = MathSqrt(se);

// 4. Critical value of the t-distribution
   double t_crit = 1.96; // For a 95% confidence interval 

// 5. Confidence intervals
   lower = w_ols - t_crit * se;
   upper = w_ols + t_crit * se;
  }

//+------------------------------------------------------------------+
//| Credible intervals, Bayesian approach                            |
//+------------------------------------------------------------------+
void CredibleIntervals(matrix &samples, double alpha, vector &lower, vector &upper)
  {
   int d = (int)samples.Cols(); 
   lower.Resize(d);
   upper.Resize(d);

   for(int j = 0; j < d; j++)
     {      
      vector col = samples.Col(j);
      int n = (int)col.Size(); 

      double temp[];
      col.Swap(temp);

      // Array of probabilities for quantiles
      double probs[] = {alpha / 2, 1 - alpha / 2}; // {0.025, 0.975} for a 95% credible interval
      double quantiles[];

      MathQuantile(temp, probs, quantiles);

      lower[j] = quantiles[0]; // alpha/2 quantile
      upper[j] = quantiles[1]; // 1-alpha/2 quantile
     }
  }

//+------------------------------------------------------------------+
//|  Save the matrix to CSV                                          |
//+------------------------------------------------------------------+
bool MatrixToCSV(string file_name, const matrix &m)
  {
   if(m.Rows() == 0 || m.Cols() == 0)
     {
      Print("Error: Matrix is empty");
      return false;
     }
   int file_handle = FileOpen(file_name, FILE_WRITE | FILE_TXT | FILE_ANSI);
   if(file_handle == INVALID_HANDLE)
     {
      Print("File opening/creating error: ", GetLastError());
      return false;
     }
   int rows = (int)m.Rows();
   int cols = (int)m.Cols();
   for(int i = 0; i < rows; i++)
     {
      string line = "";
      for(int j = 0; j < cols; j++)
        {
         line += DoubleToString(m[i][j], 10);
         if(j < cols - 1)
           {
            line += ",";
           }
        }
      FileWriteString(file_handle, line + "\n");
     }
   FileClose(file_handle);
   Print("Matrix successfully saved to file: ", file_name);
   return true;
  }
//+------------------------------------------------------------------+

Na abordagem frequentista, as estimativas dos parâmetros são calculadas pelo método dos mínimos quadrados ordinários (Ordinary Least Squares), enquanto os intervalos de confiança são obtidos com base na distribuição t.

Na abordagem bayesiana, as estimativas pontuais são obtidas pela média da distribuição posterior, isto é, pela média das amostras, enquanto a incerteza é representada por intervalos de credibilidade de 95%, calculados pela função CredibleIntervals a partir dos quantis 0.025 e 0.975.

Ao final da amostragem MCMC, compararemos as estimativas obtidas pelas abordagens frequentista e bayesiana. Os resultados mostram que conseguimos gerar corretamente amostras da distribuição posterior dos parâmetros e que os valores obtidos são muito próximos das estimativas frequentistas. Esse comportamento é esperado quando se utiliza uma distribuição a priori fracamente informativa, normalmente caracterizada por uma variância elevada, como ocorre neste caso.

Os valores médios obtidos com o sampler são praticamente idênticos às estimativas OLS. Isso indica que o algoritmo convergiu e explorou corretamente a região relevante da densidade da distribuição posterior. Embora todas as estimativas estejam ligeiramente deslocadas em relação aos valores verdadeiros [1.0, 2.0, 2.0], isso é esperado para uma amostra relativamente pequena, com N = 100, especialmente na presença de ruído. Ao executar o script, suas estimativas poderão diferir um pouco devido à nova geração de números aleatórios.

Parâmetro	Valor verdadeiro	OLS	SS
w_0 (intercepto, b)	1	0.91847	0.91852
w_1	2	1.97694	1.97688
w_2	2	2.02333	2.02320

Os intervalos de credibilidade bayesianos (Credible Intervals) também são muito próximos dos intervalos de confiança frequentistas (Confidence Intervals).

Parâmetro	95% Confidence Intervals (abordagem frequentista)	95% Credible Intervals (abordagem bayesiana)
w_0 (intercepto, b)	[0.8220, 1.0149]	[0.8227, 1.0175]
w_1	[1.9415, 2.0123]	[1.9409, 2.0128]
w_1	[1.9892, 2.0574]	[1.9883, 2.0585]

A eficiência do algoritmo slicesample é avaliada pelo número médio de chamadas à função-alvo log f(x). Com essas configurações, o parâmetro neval da função slicesample ficou em torno de 5.365.

Para ajustar esse indicador, é necessário alterar o valor do parâmetro width. Se width for muito pequeno, o algoritmo fará um número excessivo de avaliações da função para determinar a extensão da fatia. Por outro lado, se width for muito grande, o algoritmo precisará reduzir frequentemente o intervalo até chegar a um tamanho adequado, o que também aumentará desnecessariamente o número de avaliações da função.

Além da análise de neval, uma avaliação completa da qualidade da amostragem exige examinar o gráfico de traço das amostras e a função de autocorrelação (ACF), bem como construir histogramas para cada parâmetro do modelo.

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

# File path
slice_file = "C:/Program Files/MetaTrader 5/MQL5/Files/SS/LR_samples.csv"

# Load data
slice_logpdf_data = pd.read_csv(slice_file, header=None)
samples = slice_logpdf_data.values  # nsamples x 3 matrix
dim = samples.shape[1]  # Number of dimensions (3: b, w_1, w_2)

# Expected parameter values
expected_values = [1.0] + [2.0] * (dim - 1)  # [1.0, 2.0, 2.0]

# Statistics
def print_stats(samples, dim, expected_values):
    print("=========================================")
    print(f"MCMC Samples Analysis (Bayesian Linear Regression, {dim}D):")
    print(f"Number of samples: {samples.shape[0]}")
    for i in range(dim):
        param_name = "b" if i == 0 else f"w_{i}"
        print(f"Mean {param_name}: {np.mean(samples[:, i]):.4f} (expected {expected_values[i]:.4f})")
        print(f"Std. deviation {param_name}: {np.std(samples[:, i]):.4f}")
    # Covariance between b and w_1 (for example)
    print(f"Covariance (b, w_1): {np.cov(samples[:, 0], samples[:, 1])[0, 1]:.4f}")
    print("=========================================")

print_stats(samples, dim, expected_values)

# Plot traces and ACF
fig, axes = plt.subplots(dim, 2, figsize=(12, 4 * dim))
fig.suptitle('Trace Plots and ACF for Bayesian Linear Regression Parameters', fontsize=16)

for i in range(dim):
    param_name = "b" if i == 0 else f"w_{i}"
    # Trace Plot
    axes[i, 0].plot(samples[:, i])
    axes[i, 0].set_title(f'Trace Plot ({param_name})')
    axes[i, 0].set_xlabel('Iteration')
    axes[i, 0].set_ylabel(param_name)
    axes[i, 0].grid(True)
    # Add a horizontal line for the expected value
    axes[i, 0].axhline(y=expected_values[i], color='r', linestyle='--', label=f'Expected {expected_values[i]}')
    axes[i, 0].legend()

    # ACF Plot
    plot_acf(samples[:, i], lags=50, ax=axes[i, 1], title=f'ACF ({param_name})')
    axes[i, 1].set_xlabel('Lags')
    axes[i, 1].set_ylabel('Correlation')

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

# Plot histograms of posterior distributions
fig, axes = plt.subplots(1, dim, figsize=(12, 3))
fig.suptitle('Posterior Distributions for Bayesian Linear Regression Parameters', fontsize=16)

for i in range(dim):
    param_name = "b" if i == 0 else f"w_{i}"
    axes[i].hist(samples[:, i], bins=30, density=True, alpha=0.7, color='skyblue')
    axes[i].set_title(f'Posterior ({param_name})')
    axes[i].set_xlabel(param_name)
    axes[i].set_ylabel('Density')
    # Add a vertical line for the expected value
    axes[i].axvline(x=expected_values[i], color='r', linestyle='--', label=f'Expected {expected_values[i]}')
    axes[i].legend()

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

O gráfico da densidade posterior do parâmetro de intercepto w0 é apresentado na Fig. 3.

Fig. 3. Densidade posterior do parâmetro de intercepto do modelo de regressão linear

Exemplo nº 2. Treinamento de uma regressão logística bayesiana

Neste exemplo, aplicamos o algoritmo slice sampling à modelagem de resultados binários y ∈ {0,1} em um problema de regressão logística bayesiana. Geraremos dados sintéticos com os parâmetros verdadeiros w_true = [0.0, 1.0, −1.0] e, em seguida, tentaremos recuperar esses valores por meio de métodos frequentistas e bayesianos.

Primeiro, definimos o modelo dos dados, isto é, a função de verossimilhança. A verossimilhança assume a forma de uma distribuição de Bernoulli parametrizada pela função sigmoide aplicada ao preditor linear Xw:

Como distribuição a priori para os parâmetros w, utilizaremos uma distribuição normal multivariada:

P(w) = N (w | 0, 10 * I)

Assim como no exemplo anterior, escolheremos uma distribuição a priori fracamente informativa, com variância relativamente alta. Vale lembrar que, em estatística bayesiana, uma distribuição a priori é considerada fracamente informativa quando exerce pouca influência sobre a distribuição posterior, permitindo que a verossimilhança tenha o papel predominante na inferência.

Quanto à abordagem frequentista clássica, as estimativas pontuais podem ser obtidas pelo algoritmo IRLS (Iteratively Reweighted Least Squares). O IRLS é um método iterativo para resolver numericamente problemas de otimização, bastante adequado para modelos de regressão logística. Ele minimiza a função de perda logarítmica, o que equivale a maximizar a verossimilhança. A atualização dos parâmetros w é calculada pela fórmula:

onde:

• X é a matriz de atributos (n×(d+1)), incluindo uma coluna de uns para o intercepto,
• W é uma matriz diagonal de pesos de dimensão n×n, cujos elementos são wii = pi (1 − pi), sendo pi = σ(Xiw) as probabilidades calculadas pela função sigmoide,
• z = Xw + (y − p) / p*(1 − p) é a chamada variável de trabalho, em que y é o vetor de respostas binárias observadas e p é o vetor de probabilidades previstas,
• X'WX é a matriz de informação de Fisher, utilizada para calcular a matriz de covariância e os erros-padrão dos parâmetros.

Os intervalos de confiança são construídos com base na matriz de informação de Fisher. Vamos compará-los com as estimativas bayesianas obtidas por meio do nosso sampler.
Parâmetros do sampler:

• nsamples 10000 - número de amostras geradas,
• burnin 1000 - período de aquecimento da cadeia,
• thin 10 - intervalo de thinning,
• width 1 - largura inicial do hiper-retângulo.

//+------------------------------------------------------------------+
//|                                                    LogisticR.mq5 |
//|                                                           Eugene |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Eugene"
#property link      "https://www.mql5.com"
#property version   "1.00"

#include <Math\Stat\Math.mqh>
#include <Math\Stat\Normal.mqh>
#include <Math\Stat\Uniform.mqh>
#include <MCMC\SS.mqh>

// Data storage structure 
struct BayesianLogisticRegression
  {
   matrix            X;           // Feature matrix (n × (d+1))
   vector            y;           // Binary responses (n)
   vector            mu_w;        // Prior mean for w = [b, w₁, w₂]
   matrix            Sigma_w;     // The prior covariance matrix for w
   matrix            Sigma_w_inv; // Inverse covariance matrix 
   double            log_det_Sigma_w; // Log determinant of Sigma_w
  };

BayesianLogisticRegression model;

//+------------------------------------------------------------------+
//| Sigmoid function                                                 |
//+------------------------------------------------------------------+
double sigmoid(double z)
  {
   return 1.0 / (1.0 + MathExp(-z));
  }

//+------------------------------------------------------------------+
//| Log-density of a multivariate normal distribution                |
//+------------------------------------------------------------------+
double multivariate_normal_logpdf(vector &x, vector &mu, matrix &Sigma_inv, double log_det, int dim)
  {
   vector diff = x - mu;
   double quad_form = diff @ Sigma_inv @ diff;
   return -0.5 * (quad_form + log_det + dim * MathLog(2 * M_PI));
  }

//+------------------------------------------------------------------+
//| Target density - posterior distribution                          |
//+------------------------------------------------------------------+
double LogPosterior(vector &params, BayesianLogisticRegression &mdl)
{
    int d = (int)model.X.Cols(); // Number of parameters (including intercept)
   
    // --- 1. Log-likelihood ---
    int n = (int)model.X.Rows();
    
    // Calculate the linear predictor: eta = X * w
    vector eta = model.X @ params; 
    
    // Calculate the probability vector p = sigmoid(eta)
    vector p(n);
    eta.Activation(p, AF_SIGMOID);
   
    // Calculate the vector log(p)
    vector log_p(n);
    log_p =  MathLog(p);
    
    // Calculate the vector log(1-p)
    vector one_minus_p_log(n);
    one_minus_p_log = MathLog(1.0 - p);
    
    // Calculate the log-likelihood: L = sum( y * log(p) + (1-y) * log(1-p) )
    vector term_1 = mdl.y * log_p;
    vector term_2 = (1.0 - mdl.y) * one_minus_p_log;
    double log_likelihood = (term_1 + term_2).Sum(); 

    // --- 2. Log-prior (multivariate normal prior) ---
    double log_prior = multivariate_normal_logpdf(params, mdl.mu_w, mdl.Sigma_w_inv,
                                                 mdl.log_det_Sigma_w, d);
    // Posterior 
    return log_likelihood + log_prior;
}

// A wrapper to match the LogPdfFunction signature
double LogPost(vector &params)
  {
   return LogPosterior(params, model);
  }
  
//+---------------------------------------------------------------------+
//|Estimation using IRLS (Iteratively Reweighted Least Squares)         |
//+---------------------------------------------------------------------+
void IRLS(matrix &X, vector &y, vector &w_irls, vector &se_irls)
{
    int n = (int)X.Rows(); // Number of observations
    int d = (int)X.Cols(); // Number of parameters (including intercept)

    // --- Algorithm parameters ---
    double convergence_threshold = 1e-6; // Threshold for checking the convergence of weights 
    int max_iterations = 20;             // Maximum number of iterations
    int iterations_used = 0;             // Actual number of iterations used

    w_irls.Resize(d);
    w_irls.Fill(0.0);      // Initial estimate of weights w 
    vector w_old = w_irls; // A vector for storing the weights from the previous iteration

    // Parameter update formula: w_new = (X^T W X)^(-1) X^T W z
    matrix W_final(n, n);
    W_final.Fill(0.0);

    for(int iter = 0; iter < max_iterations; iter++)
    {
        iterations_used = iter + 1;
        vector p(n);    // Probability vector p_i = sigmoid(eta_i)
        vector z(n);    // Working variable vector 
        matrix W_current(n, n); // Diagonal weight matrix W
        W_current.Fill(0.0);

        // 1. Calculate the linear predictor: eta = X * w
        vector eta = X @ w_irls;

        // 2. Calculate probabilities: p = sigmoid(eta)
        eta.Activation(p, AF_SIGMOID);

        // 3. Calculate the diagonal weight vector: w_diag = p * (1 - p)
        vector w_diag = p * (1.0 - p);

        // 4. Create the diagonal weight matrix W_current
        W_current.Diag(w_diag); 

        // 5. Calculate the working variable z
        vector diff = y - p;
        z = eta + diff / w_diag; 

        // Save the current W for use in the SE calculation
        W_final = W_current; 

        // 6. Weight update (IRLS step)
        matrix XtW = X.Transpose() @ W_current;

        // Calculate XtWX (Fisher information matrix)
        matrix XtWX = XtW @ X;

        // Solve for w: w = (XtWX)^-1 @ XtW @ z
        w_old = w_irls; // Save w to check for convergence
        w_irls = XtWX.Inv() @ XtW @ z;

        // 7. Convergence check
        // Compute the L2 norm of the difference between the weights
        double diff_norm = (w_irls - w_old).Norm(VECTOR_NORM_P, 2);
        double w_norm = w_old.Norm(VECTOR_NORM_P, 2);
        
        double relative_change = (w_norm > 0.0) ? diff_norm / w_norm : diff_norm;

        if (relative_change < convergence_threshold)
        {
          Print("IRLS: Convergence achieved on iteration ", iterations_used, ". Relative change: ", DoubleToString(relative_change, 8));
          break;
        }
        
    }
    
    if (iterations_used < max_iterations)
    {
      Print("IRLS: Algorithm completed successfully (convergence). Iterations used: ", iterations_used, " out of ", max_iterations);
    }
    else
    {
      Print("IRLS: Algorithm terminated due to maximum iteration limit (", max_iterations, "). Convergence check: not achieved.");
    }

    // --- Calculation of Standard Errors (SE) ---   
    // Calculation of XtWX (Fisher information matrix) using the final weights W_final
    matrix XtWX_final = X.Transpose() @ W_final @ X; 

    // Covariance matrix Cov(w) = (X^T W_final X)^-1
    matrix cov_matrix = XtWX_final.Inv(); 

   // standard errors 
    se_irls = cov_matrix.Diag();
    se_irls = MathSqrt(se_irls);
}

//+---------------------------------------------------------------------+
//| Initialize logistic regression and sampling with SS                 |
//+---------------------------------------------------------------------+
void BayesianLogisticRegressionSample(matrix &X, vector &y,
                                      vector &mu_w, matrix &Sigma_w,
                                      int nsamples, vector &initial, vector &width,
                                      int burnin, int thin, matrix &samples)
  {
   int d = (int)X.Cols();

// Precompute the inverse matrix and the log determinant
   matrix Sigma_w_inv = Sigma_w.Inv();
   double det = Sigma_w.Det();
   double log_det = MathLog(det);

// Model initialization
   model.X = X;
   model.y = y;
   model.mu_w = mu_w;
   model.Sigma_w = Sigma_w;
   model.Sigma_w_inv = Sigma_w_inv;
   model.log_det_Sigma_w = log_det;

   // Count the number of logpdf evaluations and execution time
   double neval = 0;
   ulong start_time = GetMicrosecondCount(); 
   slicesample(initial, nsamples, width, burnin, thin, samples, neval, LogPost);
   ulong end_time = GetMicrosecondCount(); 
   double time_ms = (end_time - start_time) / 1000.0; // Time in milliseconds
   Print("Average number of logpdf evaluations per sample: ", neval);
   Print("Slice sampling execution time: ", time_ms, " ms");
  }

//+------------------------------------------------------------------+
//| Credible intervals, Bayesian approach                            |
//+------------------------------------------------------------------+
void CredibleIntervals(matrix &samples, double alpha, vector &lower, vector &upper)
  {
   int d = (int)samples.Cols();
   lower.Resize(d);
   upper.Resize(d);

   for(int j = 0; j < d; j++)
     {
      vector col = samples.Col(j);
      int n = (int)col.Size();
      double temp[];
      ArrayResize(temp, n);
      for(int i = 0; i < n; i++)
        {
         temp[i] = col[i];
        }
      double probs[] = {alpha / 2, 1 - alpha / 2};
      double quantiles[];
      
      MathQuantile(temp, probs, quantiles);

      lower[j] = quantiles[0];
      upper[j] = quantiles[1];
     }
  }
//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
// Synthetic data generation
   int n = 100;  // Number of observations
   int d = 2;    // Number of features (excluding the intercept parameter)
   matrix X(n, d + 1); // Feature matrix X including the intercept parameter
   vector y(n);
   int err;
   vector true_w(d + 1); // w = [b, w_1, w_2]
   true_w[0] = 0.0;      // True intercept
   true_w[1] = 1.0;      // True weight w_1
   true_w[2] = -1.0;     // True weight w_2

// Populate the matrix X and vector y
   for(int i = 0; i < n; i++)
     {
      X[i][0] = 1.0; // The first column is a column of ones
      for(int j = 1; j < d + 1; j++)
        {
         X[i][j] = MathRandomNormal(0.0, 1.0, err); // Features drawn from N(0, 1)
        }
      double z = true_w @ X.Row(i); 
      double p = sigmoid(z);        
      y[i] = MathRandomUniform(0.0, 1.0, err) < p ? 1.0 : 0.0; 
     }

   MatrixToCSV("SS/X_logistic.csv", X);
   matrix y_matrix(n, 1);
   for(int i = 0; i < n; i++)
     {
      y_matrix[i][0] = y[i];
     }
   MatrixToCSV("SS/y_logistic.csv", y_matrix);

// Prior parameters
   vector mu_w(d + 1);
   mu_w.Fill(0.0); // Prior mean for w = [b, w_1, w_2]
   matrix Sigma_w(d + 1, d + 1);
   for(int i = 0; i < d + 1; i++)
     {
      Sigma_w[i][i] = 10.0; // Diagonal covariance matrix (variance = 10)
     }

// Sampling parameters
   int nsamples = 10000;
   int burnin = 1000;
   int thin = 10;
   vector initial(d + 1);
   initial.Fill(0.0); // Initial values [b, w_1, w_2]
   vector width(d + 1);
   width.Fill(1);   // Step width for slice sampling
   matrix samples;

// Frequentist approach: IRLS parameter estimates
   Print("------------   Frequency approach -----------------");
   vector w_irls, se_irls;
   IRLS(X, y, w_irls, se_irls);
   Print("IRLS Estimator w: ", w_irls, " (Expected values [ 0.0, 1.0, -1.0 ])");

// Confidence intervals
   double z_crit = 1.96; // For a 95% confidence interval
   vector lower_ci_freq = w_irls - z_crit * se_irls;
   vector upper_ci_freq = w_irls + z_crit * se_irls;
   Print("95% Confidence intervals:");
   for(int i = 0; i < (int)lower_ci_freq.Size(); i++)
     {
      Print("w_", i, ": [", lower_ci_freq[i], ", ", upper_ci_freq[i], "]");
     }
     
   // Save w_irls to CSV
   matrix w_irls_matrix(1, d + 1);
   for(int i = 0; i < d + 1; i++)
      w_irls_matrix[0][i] = w_irls[i];

   MatrixToCSV("SS/w_irls_logistic.csv", w_irls_matrix);  

   Print("------------   Bayesian approach -----------------");
   BayesianLogisticRegressionSample(X, y, mu_w, Sigma_w, nsamples, initial, width, burnin, thin, samples);
  
   vector mean_w(d + 1);
   for(int i = 0; i < d + 1; i++)
     {
      mean_w[i] = samples.Col(i).Mean();
     }
   Print("Average w: ", mean_w, " (Expected values[ 0.0, 1.0, -1.0 ])");
   
// Confidence intervals
   vector lower_ci_bayes, upper_ci_bayes;
   CredibleIntervals(samples, 0.05, lower_ci_bayes, upper_ci_bayes);
   Print("95% Credible intervals:");
   for(int i = 0; i < (int)lower_ci_bayes.Size(); i++)
     {
      Print("w_", i, ": [", lower_ci_bayes[i], ", ", upper_ci_bayes[i], "]");
     }

   MatrixToCSV("SS/LR_samples_logistic.csv", samples);
  }
//+------------------------------------------------------------------+

As duas abordagens produzem estimativas próximas dos valores verdadeiros, embora apresentem alguns desvios devido à quantidade limitada de dados. A Fig. 4 mostra o gráfico da densidade posterior do parâmetro w1. O valor verdadeiro do parâmetro (w1 = 1) está dentro do intervalo de credibilidade bayesiano de 95%.

Fig. 4. Distribuição posterior do parâmetro w1 da regressão logística bayesiana

Para confirmar visualmente a consistência dos resultados, construiremos a fronteira de decisão, que na regressão logística é definida pela equação x^Tw = 0. Para dois atributos (x1, x2) e um intercepto, essa fronteira é uma reta.

Construiremos três fronteiras:

• Fronteira verdadeira, baseada nos valores verdadeiros dos parâmetros w_true = [0.0, 1.0, −1.0].
• Fronteira IRLS, baseada nas estimativas pontuais.
• Fronteira bayesiana, baseada nos valores médios da distribuição posterior.

Fig. 5. Fronteiras de decisão bayesiana e IRLS (coincidentes)

As fronteiras de decisão obtidas pelo IRLS e pela abordagem bayesiana praticamente coincidem, mas se afastam da fronteira verdadeira devido ao volume limitado de dados (n = 100) e à natureza estocástica da geração das classes por meio da função sigmoide. Isso provoca um pequeno viés nas estimativas dos parâmetros. Ainda assim, o algoritmo de amostragem por fatias foi bem-sucedido, identificando corretamente os parâmetros que proporcionam a separação entre as classes.

Conclusão

Neste artigo, estudamos o método de amostragem por fatias (slice sampling), uma variante adaptativa dos métodos MCMC que elimina a necessidade de um ajuste cuidadoso de hiperparâmetros, como a escala do passo no algoritmo de Metropolis.

A eficiência do método foi avaliada por meio de modelos bayesianos de regressão linear e logística. As estimativas bayesianas, incluindo as médias das distribuições posteriores e os intervalos de credibilidade, ficaram próximas dos resultados obtidos pelos métodos frequentistas clássicos, ou seja, mínimos quadrados ordinários (OLS) e mínimos quadrados iterativamente reponderados (IRLS). Os gráficos de traço das amostras, as funções de autocorrelação (ACF) e os histogramas das densidades posteriores confirmam a qualidade da amostragem.

Assim, a implementação do sampler slicesample em MQL5 mostrou-se uma solução universal e confiável do tipo "caixa-preta". Com esse método, a inferência bayesiana se torna consideravelmente mais simples e acessível. Na prática, o usuário precisa apenas saber codificar a densidade posterior-alvo. Isso permite concentrar os esforços diretamente no problema de inferência estatística, minimizando o esforço necessário para ajustar o algoritmo MCMC.

Programas utilizados no artigo:

#	Nome	Tipo	Descrição
1	SS.mqh	Arquivo de inclusão	Algoritmo de amostragem por fatias
2	LR.mq5	Script	Exemplo de amostragem da distribuição posterior de um modelo de regressão linear bayesiana
3	LR_plot.py	Script	Diagnóstico em Python
4	LogisticR.mq5	Script	Exemplo de amostragem da distribuição posterior de um modelo de regressão logística bayesiana
5	LogisticR_plot.py	Script	Diagnóstico em Python
6	PlotMM.mq5	Script	Animação do funcionamento do algoritmo para uma distribuição unidimensional multimodal