MCMC Sampling Methods: The Slice Sampling Algorithm

Introduction

Markov chain Monte Carlo (MCMC) methods are widely used for sampling from complex multidimensional distributions, particularly in Bayesian statistics. Classical algorithms, such as Gibbs sampling or the Metropolis algorithm, often require careful fine-tuning. In the first case, it is necessary to derive analytical expressions for all full conditional distributions, and in the second, to painstakingly select the scale and shape of the proposal distribution. This makes it difficult to use them quickly and effectively in everyday practice.

This article examines slice sampling — an adaptive variant of MCMC. Its key advantage is that it automatically adapts to the characteristics of the target distribution, eliminating the need to manually adjust parameters such as step width.

After providing a detailed description of the algorithm and its implementation in the MQL5 environment, we will test it using the examples of Bayesian linear and logistic regression, comparing the results with classical frequentist methods that provide point estimates of the parameters.

Slice Sampling Methods

One of the main challenges of MCMC algorithms, such as the Metropolis algorithm, is choosing an appropriate step size. A step size that is too small slows down the chain's movement, increasing the autocorrelation between samples and requiring many iterations to obtain independent samples. A step size that is too large leads to frequent rejection of proposed points, reducing the algorithm's efficiency.

To address this problem, a technique known as slice sampling (Neal, 2003) was proposed, which requires only the ability to evaluate the unnormalized density f(x). These methods can be divided into two types depending on the approach to multidimensional distributions:

• one-dimensional slice sampling with coordinate-wise updating,
• multidimensional slice sampling, which directly updates the entire vector of variables at once.

One-Dimensional Slice Sampling

One-dimensional slice sampling is used either for one-dimensional target distributions or for coordinate-wise updating of a multidimensional vector x = (x1, …, xn), similarly to Gibbs sampling. Let's take a look at how this algorithm works.

Updating the current point x0 to the new value x1 involves three main steps:

• (a) Defining the slice: we choose an auxiliary variable y ∼ Uniform(0, f(x0)) and define the horizontal slice S = {x: f(x) > y} (that is, the set of points x for which the density value is greater than the value of the auxiliary variable y).
• (b) Forming the interval: find an interval (L, R) of width w randomly placed around x0. The interval is expanded in increments of w (the “stepping-out” procedure) until both ends of the interval lie outside the slice.
• (c) Selecting a new point: the new point x1 is chosen by sampling uniformly from the interval (L, R). If the point lies outside the slice, the interval is shrunk (the “shrinkage” procedure), and the selection is repeated until the point x1 falls within the slice.

In practice, it is more convenient to work on a logarithmic scale:

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

where e ∼ Exp (1) is an exponentially distributed random variable with a mean of one. Then the slice S = {x: g(x) > z}

The “stepping-out” and “shrinkage” procedures proposed by Radford Neal (Neal, 2003), a renowned expert in machine learning and MCMC methods, ensure efficient and correct sampling in the one-dimensional case. A graphical illustration of this algorithm is shown in Fig. 1.

Fig. 1. One-dimensional slice sampling update using the “stepping-out” and “shrinkage” procedures.

You can view an animation of the algorithm in action in the one-dimensional case by running the following script, which demonstrates the process of sampling from a multimodal distribution (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 for a One-Dimensional Distribution

The graph shows the initial point x0 (marked with a blue circle) and the density value at that point, f(x0) (red circle). The resulting vertical slice is shown as a red dashed line. A point z is selected uniformly from this slice to define the horizontal slice (the thin black line). An interval of width w, randomly placed around x0, is expanded in increments of w until its endpoints lie outside the slice (marked by the bold blue line). A new point x1 is selected uniformly from this interval until a point within the slice is found. Points selected outside the slice are used to shrink the interval. The graph shows a new candidate point, x1 (the red square), which the algorithm rejects because it lies outside the slice.

Multidimensional Slice Sampling

Instead of coordinate-wise updating of the vector x = (x1, …, xn), the idea of slice sampling can be applied directly to a multidimensional distribution. In this case, the one-dimensional interval (L, R) is replaced by the hyperrectangle H = {x: Li < xi < Ri, i = 1, …, n}, where Li and Ri define the boundaries along each axis.

The procedure for finding the next state x1=(x1,1,…,x1,n) from the current state x0=(x0,1,…,x0,n) is similar to the one-dimensional case:

• (a) Choose y uniformly from (0, f(x0)), defining the slice S = {x: y < f(x)}
• (b) Find a hyperrectangle H = (L1, R1) ×⋯× (Ln, Rn) surrounding x0 that, preferably, contains most of the slice.
• (c) Select a new point x1 uniformly from the part of the slice inside H. If a point lies outside the slice, the hyperrectangle is shrunk (the "shrinkage" procedure).

Unlike the one-dimensional case, the “stepping-out” procedure is very complex in multidimensional space, since it requires checking 2^n vertices of the hyperrectangle, which becomes computationally expensive when n is large. Therefore, a simplified approach is often used: the hyperrectangle is randomly placed around x0 without performing the expansion procedure. The hyperrectangle is shrunk independently along each axis until a point is found inside the slice. This method is easier to implement, but less efficient than one-dimensional slice sampling. With coordinate-wise updating, each coordinate is shrunk exactly as much as necessary, taking into account local features of the density. In the multidimensional case, all axes are shrunk simultaneously, which can lead to excessive shrinking in directions where the density changes slowly.

In the code shown below, the "stepping-out" procedure is used only for the one-dimensional case. For a multivariate distribution, a simplified approach without "stepping-out" is used, relying on the random placement of a hyperrectangle and its subsequent shrinking.

//+------------------------------------------------------------------+
//|                                                           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_);
}

Example No. 1. Training a Bayesian Linear Regression Model

In this example, we will examine the use of the slice sampling algorithm to draw samples from the posterior distribution of the parameters in Bayesian linear regression. To apply Bayes' theorem, let's define a data model by specifying the likelihood function. To avoid overcomplicating the model, let’s assume that the data y are normally distributed around a linear combination of features Xw:

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

where:

• I is the identity matrix,
• σ² — the noise variance, which we will assume to be known for simplicity.

The unknowns in this model are the parameters w. Under the Bayesian approach, all unknown quantities are treated as random variables. Therefore, we need to specify a distribution for these random variables. This distribution is called a prior distribution, since we specify it before the model sees the data. For this distribution, we can use a multivariate normal distribution with zero mean and covariance matrix Σw.

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

Once we receive the data, we update our prior beliefs about the parameters by calculating the posterior distribution as the product of the likelihood and the prior:

p(w | X, y) ∝ p(y | Xw, σ²) · p(w)

In the code, the LogPosterior function is responsible for computing the posterior, and we pass it to our sampler to generate samples.

Before we begin sampling, let’s generate synthetic data with known parameters w_true = [1.0, 2.0, 2.0], a noise standard deviation of σ = 0.5, and a sample size of 100.

Let's set the following parameters for our sampler:

• nsamples 10000 — number of samples to generate,
• burnin 1000 — burn-in period,
• thin 10 — thinning interval,
• width 0.5 — initial width of the hyperrectangle

//+------------------------------------------------------------------+
//|                                                           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;
  }
//+------------------------------------------------------------------+

In the frequentist approach, parameter estimates are computed using ordinary least squares (OLS), and confidence intervals are calculated based on the t-distribution.

In the Bayesian approach, point estimates are obtained as the mean of the posterior distribution (the sample mean), and uncertainty is represented by 95% credible intervals computed by the CredibleIntervals function for the 0.025 and 0.975 quantiles.

Once the MCMC sampling is complete, let us compare the estimates for the frequentist and Bayesian approaches. The results show that we have successfully sampled from the posterior distribution of the parameters, and the resulting estimates are in close agreement with the frequentist estimates, which is expected when using a weakly informative prior (typically distributions with large variance, as in our case).

The means obtained using the sampler are virtually identical to the OLS estimates. This is evidence that the algorithm has converged and correctly explored the region of the posterior density. Although all estimates differ slightly from the true values [1.0, 2.0, 2.0], this is normal for such a small sample of size N = 100, especially given the presence of noise. If you run the script, your estimates will differ slightly because a new set of random numbers will be generated.

Parameter	true value	OLS	SS
w_0 (intercept / bias, b)	1	0.91847	0.91852
w_1	2	1.97694	1.97688
w_2	2	2.02333	2.02320

Bayesian credible intervals are also very similar to frequentist confidence intervals.

Parameter	95% Confidence Intervals (frequentist approach)	95% Credible Intervals (Bayesian approach)
w_0 (intercept / bias, 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]

The efficiency of the `slicesample` algorithm is measured by the average number of calls to the target function `log f(x)`. With these settings, the `neval` parameter in the `slicesample` function was approximately 5.365.

To adjust this metric, you need to change the value of the `width` parameter. Keep in mind that if the `width` parameter is too small, the algorithm will perform an excessive number of function evaluations to determine the size of the slice. If, on the other hand, the width is too large, the algorithm will have to frequently shrink the interval to a suitable value, which will also result in an excessive number of function evaluations.

In addition to assessing `neval`, a comprehensive analysis of sampling quality requires examining the trace plots and the autocorrelation function (ACF), as well as plotting histograms for each model parameter.

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()

The posterior density plot for the intercept parameter w0 is shown in Fig. 3.

Fig. 3. Posterior density of the intercept parameter in the linear regression model

Example No. 2. Training a Bayesian Logistic Regression Model

In this example, we apply the slice sampling algorithm to model binary outcomes y ∈ {0,1} in a Bayesian logistic regression problem. We will generate synthetic data using the true parameters w_true = [0.0, 1.0, −1.0], and then attempt to recover these true values using frequentist and Bayesian methods.

First, let's define the data model — that is, the likelihood function. The likelihood takes the form of a Bernoulli distribution parameterized by the sigmoid function of the linear predictor Xw:

As a prior for the parameters w, let’s use a multivariate normal distribution:

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

As in the previous example, we will choose a weakly informative prior with fairly high variance. As a reminder, in Bayesian statistics, a prior is considered weakly informative if it has a minimal effect on the posterior distribution, allowing the likelihood to play the primary role in drawing conclusions.

As for the classical frequentist approach, point estimates can be obtained using the IRLS (Iteratively Reweighted Least Squares) algorithm. IRLS is an iterative method for numerically solving the optimization problem and is well suited to logistic regression models. It minimizes the log-loss function, which is equivalent to maximizing the likelihood. The update of the parameters w is computed using the following formula:

where:

• X – a feature matrix (n×(d+1)), including a column of ones for the intercept,
• W is an n×n diagonal weight matrix, with diagonal elements wii = pi(1 − pi), where pi = σ(Xiw) are the probabilities computed using the sigmoid function,
• z = Xw + (y − p) / (p*(1 − p)) — the so-called working variable, where y is the vector of observed binary responses and p is the vector of predicted probabilities
• X'WX — the Fisher information matrix, used to compute the covariance matrix and the standard errors of the parameters.

Confidence intervals are constructed based on the Fisher information matrix. We will compare them with the Bayesian estimates we obtain using our sampler.
Sampler parameters:

• nsamples 10000 — the number of samples to generate,
• burnin 1000 — the burn-in period,
• thin 10 — the thinning interval,
• width 1 — the initial width of the hyperrectangle.

//+------------------------------------------------------------------+
//|                                                    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);
  }
//+------------------------------------------------------------------+

Both approaches yield estimates that are close to the true values, although with some deviations due to the finite sample size. Figure 4 shows a plot of the posterior density for the parameter w1. The true value of the parameter (w1 = 1) lies within the 95% Bayesian credible interval.

Fig. 4. The posterior distribution of the parameter w1 in Bayesian logistic regression

To visually confirm that the results are correct, let’s plot the decision boundary, which in logistic regression is defined by the equation x^Tw = 0. For two features (x1, x2) and an intercept, the decision boundary is a straight line.

Let’s plot three decision boundaries:

• The true decision boundary, based on the true values of the parameters w_true = [0.0, 1.0, −1.0].
• The IRLS decision boundary, based on point estimates.
• The Bayesian decision boundary, based on the posterior means.

Fig. 5. Bayesian and IRLS decision boundaries (overlap)

The decision boundaries for the IRLS and Bayesian approaches are virtually identical, but they deviate from the true boundary due to the limited amount of data (n=100) and the stochastic nature of class generation via the sigmoid function. This results in a slight bias in the parameter estimates. Nevertheless, the slice sampling algorithm successfully accomplished the task, correctly identifying the parameters that ensure class separation.

Conclusion

In this article, we explored slice sampling — an adaptive variant of MCMC that eliminates the need for careful tuning of hyperparameters, such as the step size in the Metropolis algorithm.

The effectiveness of the method was tested using Bayesian linear and logistic regression. Bayesian estimates (posterior means and credible intervals) were found to be close to the results obtained using classical frequentist methods — ordinary least squares (OLS) and iteratively reweighted least squares (IRLS). Trace plots, autocorrelation functions (ACF), and posterior-density histograms confirm the quality of the sampling.

Thus, the implementation of the `slicesample` sampler in MQL5 has proven to be a versatile and reliable “black box.” With this method, Bayesian inference has become much simpler and more accessible. Essentially, the user only needs to be able to encode the target posterior density. This allows us to focus directly on the task of statistical inference, minimizing the effort required to tune the MCMC algorithm.

Programs used in this article:

#	Name	Type	Description
1	SS.mqh	Include file	Slice Sampling Algorithm
2	LR.mq5	Script	Example of sampling from the posterior distribution for a Bayesian linear regression model
3	LR_plot.py	Script	Diagnostics in Python
4	LogisticR.mq5	Script	Example of sampling from the posterior distribution for a Bayesian logistic regression model
5	LogisticR_plot.py	Script	Diagnostics in Python
6	PlotMM.mq5	Script	Animation of the algorithm for a multimodal one-dimensional distribution