Выборочные методы MCMC: Алгоритм выборки по уровням (Slice sampling)
245
0
Введение
Методы Монте-Карло по схеме марковских цепей (MCMC) широко применяются для выборки из сложных многомерных распределений, особенно в байесовской статистике. Классические алгоритмы, такие как выборка Гиббса или алгоритм Метрополиса, часто требуют тщательной предварительной настройки. В первом случае, необходимо получить аналитические выражения для всех полных условных распределений, а во втором — кропотливо подбирать масштаб и форму предлагающего распределения. Это затрудняет их быстрое и эффективное использование в повседневной практике.
В данной статье исследуется метод выборки по уровням (slice sampling) — адаптивная разновидность MCMC. Его ключевое преимущество заключается в том, что он автоматически подстраивается под особенности целевого распределения, устраняя необходимость в ручной настройке параметров, таких как ширина шага.
После детального описания алгоритма и его реализации в среде MQL5, мы проверим его на примере байесовской линейной и логистической регрессии, сравнив результаты с классическими частотными методами, дающими точечные оценки параметров.
Методы выборки по уровням
Одна из главных трудностей алгоритмов MCMC, таких, как алгоритм Метрополиса, — выбор подходящего масштаба шага. Слишком маленький шаг замедляет движение цепи, увеличивая автокорреляцию между выборками и требуя множества итераций для получения независимых сэмплов. Слишком большой шаг приводит к частым отклонениям предложенных точек снижая эффективность алгоритма.
Для решения этой проблемы была предложена техника выбора по уровням или выбора срезами (Neal, 2003), которая требует лишь возможность оценивать ненормированную плотность f(x). Данные методы можно разделить на два типа в зависимости от подхода к многомерным распределениям:
- одномерная выборка по уровням с поочерёдным обновлением координат,
- многомерная выборка по уровням, которая предполагает прямое обновление всего вектора переменных одновременно.
Одномерная выборка по уровням
Одномерная выборка по уровням используется или для одномерных целевых распределений, или при поочерёдном обновлении координат многомерного вектора x = (x1,…, xn), подобно выборке Гиббса. Давайте рассмотрим механику этого алгоритма.
Обновление текущей точки x0 на новое значение x1 включает три основных шага:
- (a) Определение среза: выбираем вспомогательную переменную y ∼ Uniform (0, f(x0)) и определяем горизонтальный срез S = {x: f(x) > y} (то есть множество таких точек x для которых значение плотности больше значения вспомогательной переменной y).
- (b) Формирование интервала: находим интервал (L, R) шириной w случайно размещенный вокруг x0. Интервал расширяется с шагом w (процедура «stepping-out») до тех пор, пока оба конца интервала не окажутся вне среза.
- (с) Выбор новой точки: новая точка x1 находится путём равномерного выбора из интервала (L, R). Если точка лежит вне среза, интервал сжимается (процедура «shrinkage»), и выбор повторяется, пока точка x1 не попадёт в срез.
На практике удобнее работать в логарифмической шкале:
g(x) = log (f(x)), z = log(y) = g(x0) − e,
где e ∼ Exp (1) - экспоненциально распределенная случайная величина с математическим ожиданием равным единице. Тогда срез S = {x: g(x) > z}
Процедуры «stepping-out» и «shrinkage» предложенные Рэдфордом Нилом (Neal, 2003) известным специалистом в области машинного обучения и методов MCMC, обеспечивают эффективное и корректное сэмплирование в одномерном случае. Графическая иллюстрация данного алгоритма представлена на рис.1
Рис.1. Обновление одномерной выборки по уровням с использованием процедур «расширения» (stepping-out) и «сжатия» (shrinkage).
Анимацию работы алгоритма в одномерном случае можно посмотреть, запустив следующий скрипт, демонстрирующий процесс выборки из мультимодального распределения (рис.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; // Скорость анимации (меньше - быстрее) // Определение типа функции для logpdf typedef double (*LogPdfFunction)(vector &x); // Проверка, находится ли точка внутри среза bool inside(vector &x, double th, LogPdfFunction logpdf) { return logpdf(x) > th; } // целевая плотность 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) { // Проверка входных параметров if (burnin_ < 0 || thin_ <= 0) { Print("Ошибка: burnin должен быть >= 0, thin > 0"); return; } // Инициализация int dim = (int)initial_.Size(); rnd.Resize(nsamples_, dim); rnd.Fill(0.0); int maxiter = 200; vector x0 = initial_; neval = nsamples_; // Генерация экспоненциальных и равномерных чисел 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. График целевой плотности === 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 = ""; // Основной цикл for (int i = 1 - burnin_; i <= nsamples_ * thin_; i++) { // === 2. стартовая точка 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. Вертикальный срез === 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. Горизонтальный срез 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_); // Сохраняем имя текущей кривой prev_slice_name = current_slice_name; // === 6. Интервал [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; //--- stepping 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("Ошибка: слишком много итераций в 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("Ошибка: слишком много итераций в stepping-out (right)"); return; } neval += iter; } Sleep(time_); // === 8. xp — точка внутри среза === 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("Ошибка: слишком много итераций в 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); } //+------------------------------------------------------------------+
Рис.2. Выборка по уровням для одномерного распределения
На графике показана начальная точка x0 (отмечена голубым кружком) и значение плотности в этой точке f(x0) (красный кружок). Получившийся вертикальный срез отображается красной пунктирной линией. Из этого среза равномерно выбирается точка z, задающая горизонтальный срез (тонкая черная линия). Интервал ширины w, случайно размещенный вокруг x0, расширяется с шагом w до тех пор, пока его концы не выйдут за пределы среза (отмечено жирной голубой линией). Новая точка x1 выбирается равномерно из этого интервала, пока не будет найдена точка внутри среза. Точки, выбранные вне среза, используются для сокращения интервала. На графике как раз отображена новая точка-кандидат x1 (красный квадрат), которую алгоритм отклоняет, так как она вышла за пределы среза.
Многомерная выборка по уровням
Вместо поочерёдного обновления каждой координаты вектора x = (x1,…, xn) можно применить идею выборки срезами непосредственно к многомерному распределению. В этом случае одномерный интервал (L, R) заменяется гиперпрямоугольником H={x: Li < xi < Ri, i=1,…, n} где Li и Ri задают границы вдоль каждой оси.
Процедура нахождения следующего состояния x1=(x1,1,…,x1,n) из текущего состояния x0=(x0,1,…,x0,n) аналогична одномерной:
- (a) Выбрать y равномерно из (0,f(x0)), определяя срез S = {x: y < f(x)}
- (b) Найти гиперпрямоугольник H = (L1,R1) ×⋯× (Ln, Rn), окружающий x0, который желательно содержит большую часть среза.
- (c) Выбрать новую точку x1 равномерно из части среза внутри H. Если точка лежит вне среза, гиперпрямоугольник сжимается (процедура «shrinkage»).
В отличие от одномерного случая, процедура «stepping-out» очень сложна в многомерном пространстве, так как требует проверки 2^n вершин гиперпрямоугольника, что становится вычислительно затратным при большом n. Поэтому часто используется упрощённый подход: гиперпрямоугольник случайно размещается вокруг x0 без выполнения процедуры расширения. Сжатие гиперпрямоугольника выполняется независимо по каждой оси, пока не найдётся точка внутри среза. Этот метод проще в реализации, но менее эффективен, чем одномерная выборка. При поочерёдном обновлении каждая координата сжимается ровно настолько, насколько необходимо, учитывая локальные особенности плотности. В многомерном случае все оси сжимаются одновременно, что может привести к избыточному сжатию в направлениях с медленно меняющейся плотностью.
В коде, представленном ниже, процедура «stepping-out» используется только для одномерного случая. Для многомерного распределения используется упрощённый подход без «stepping-out», полагаясь на случайное размещение гиперпрямоугольника и его последующее сжатие.
//+------------------------------------------------------------------+ //| 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> // Определение типа функции для logpdf typedef double (*LogPdfFunction)(vector &x); // Проверка, находится ли точка внутри среза 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) { // Проверка входных параметров if (burnin_ < 0 || thin_ <= 0) { Print("Ошибка: burnin должен быть >= 0, thin > 0"); return; } // Инициализация int dim = (int)initial_.Size(); rnd.Resize(nsamples_, dim); rnd.Fill(0.0); int maxiter = 200; // максимальное кол-во итераций для каждого шага vector x0 = initial_; neval = nsamples_; // Генерация экспоненциальных случайных чисел double e[]; MathRandomExponential(1.0, nsamples_ * thin_ + burnin_, e); // равномерные св для рандомизации масштаба шага w matrix rw = matrix::Random(nsamples_ * thin_ + burnin_, dim, 0.0, 1.0); //равномерные св для выбора точки внутри интервала matrix rd = matrix::Random(nsamples_ * thin_ + burnin_, dim, 0.0, 1.0); // Основной цикл slice sampling for (int i = 1 - burnin_; i <= nsamples_ * thin_; i++) { double z = logpdf(x0) - e[i + burnin_ - 1]; // горизонтальный срез S={x:logf(x)>z} // Начальный интервал: vector r = width_ * rw.Row(i + burnin_ - 1); vector xl = x0 - r; vector xr = xl + width_; int iter = 0; //--- step out только для одномерного распределения ----------- if (dim == 1) { // step out to the left while (inside(xl, z, logpdf) && iter < maxiter) { xl -= width_; iter++; } if (iter >= maxiter) { Print("Ошибка: слишком много итераций в 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("Ошибка: слишком много итераций в stepping-out"); return; } neval += iter; } //--- Shrinking --- vector xp = rd.Row(i + burnin_ - 1) * (xr - xl) + xl; // Сжимаем интервал (или гиперпрямоугольник), если выбрана точка за пределами плотности 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]; // Если xp[d] > x0[d], сжимаем правую границу else xl[d] = xp[d]; // Иначе если xp[d] <= x0[d], сжимаем левую границу } vector new_rand = vector::Random(dim, 0.0, 1.0); xp = new_rand * (xr - xl) + xl; iter++; } if (iter >= maxiter) { Print("Ошибка: слишком много итераций в shrinking"); return; } neval += iter; x0 = xp; // обновляем текущее состояние if (i > 0 && i % thin_ == 0) { rnd.Row(xp, i / thin_ - 1); } } neval /= (nsamples_ * thin_ + burnin_); }
Пример №1. Обучение байесовской линейной регрессии
В данном примере рассмотрим использование алгоритма slice sampling для получения выборки из апостериорного распределения параметров в байесовской линейной регрессии. Для применения теоремы Байеса, зададим модель данных, определив функцию правдоподобия. Чтобы слишком не усложнять модель, предположим, что данные y распределены нормально вокруг линейной комбинации признаков Xw:
P(y ∣ w,σ2) = N(y | Xw, σ^2 * I)
где,
- I — единичная матрица,
- σ2 — дисперсия шума, которую для простоты будем считать известной.
Неизвестными в данной модели являются параметры w. А в рамках байесовского подхода, все неизвестные величины трактуются как случайные величины. Следовательно, нам нужно задать распределение для этих случайных величин. Такое распределение называется априорным, так как мы задаем его до того, как модель увидит данные. В качестве такого распределения можно взять многомерное нормальное распределение с нулевым средним и матрицей ковариаций Σw.
P(w) = N (w| 0, Σw)
Получив данные, мы обновляем априорные представления о параметрах, вычисляя апостериорное распределение как произведение правдоподобия и априора:
p(w∣ X, y) ∝ p(y∣ Xw,σ2) ⋅ p(w)
В коде за вычисление апостериора отвечает функция LogPosterior, которую мы передаем нашему семплеру для генерации выборок.
Перед началом сэмплирования создадим синтетические данные с известными параметрами w_true = [1.0, 2.0, 2.0], стандартным отклонением шума σ = 0.5 и размером выборки равным 100.
Зададим следующие параметры для нашего семплера:
- nsamples 10000 - количество генерируемых выборок,
- burnin 1000 - период прогрева цепи,
- thin 10 - интервал прореживания,
- width 0.5 - начальная ширина гиперпрямоугольника
//+------------------------------------------------------------------+ //| 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> // Структура для хранения данных struct BayesianLinearRegression { matrix X; // Матрица признаков (n x (d+1)) vector y; // Выходные данные (n) double sigma; // Стандартное отклонение шума vector mu_w; // Априорное среднее для w = [b, w_1, ..., w_d] matrix Sigma_w; // Априорная ковариационная матрица для w matrix Sigma_w_inv; // Обратная ковариационная матрица double log_det_Sigma_w; // Логарифм определителя Sigma_w }; BayesianLinearRegression model; //+------------------------------------------------------------------+ //| Логарифм плотности многомерного нормального распределения | //+------------------------------------------------------------------+ 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)); } //+------------------------------------------------------------------+ //| Целевая плотность - апостериорное распределение | //+------------------------------------------------------------------+ double LogPosterior(vector ¶ms, BayesianLinearRegression &mdl) { int d = (int)mdl.X.Cols(); // Количество столбцов в X (d+1, включая единичный столбец) int n = (int)mdl.X.Rows(); // Вычислить количество наблюдений //---- Логарифм правдоподобия (Gaussian likelihood) ------------------------------------ // 1. вектор предсказаний vector y_pred = mdl.X @ params; // 2. вектор ошибок vector residuals = mdl.y - y_pred; // 3. квадрат нормы вектора ошибок (e^T e) double sum_sq_errors = residuals @ residuals; // 4. логарифм правдоподобия 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; //----------------------------------------------------------------------------------------------- // Логарифм априорного распределения (многомерный нормальный априор) double log_prior = multivariate_normal_logpdf(params, model.mu_w, model.Sigma_w_inv, model.log_det_Sigma_w, d); // апостериор = правдоподобие * априор return log_likelihood + log_prior; } double LogPost(vector ¶ms) { return LogPosterior(params, model); } //+------------------------------------------------------------------+ //| Оценки параметров по МНК | //+------------------------------------------------------------------+ 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; } //+------------------------------------------------------------------+ //| Инициализация БЛР и семплирование с помощью 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(); // Предварительно вычисляем обратную матрицу и логарифм определителя matrix Sigma_w_inv = Sigma_w.Inv(); double det = Sigma_w.Det(); double log_det = MathLog(det); // Инициализация модели 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; // Подсчёт количества оценок logpdf и времени выполнения 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; // Время в миллисекундах Print("Average number of logpdf evaluations per sample: ", neval); Print("Slice sampling execution time: ", time_ms, " ms"); } //+------------------------------------------------------------------+ //| Script program start function | //+------------------------------------------------------------------+ void OnStart() { // Генерация синтетических данных int n = 100; // Количество наблюдений int d = 2; // Количество признаков matrix X = matrix::Random(n, d+1,0,10); // Матрица признаков X X.Col(vector::Ones(n),0); vector y(n); vector true_w(d + 1); // w = [b, w_1, w_2] Истинные параметры true_w[0] = 1.0; for(int i = 1; i < d + 1; i++) { true_w[i] = 2.0; } // Центрирование предикторов 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; } } // формирования вектора y double sigma = 0.5; // Стандартное отклонение шума 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); } //--- Параметры Априора vector mu_w = vector::Zeros(d + 1); // Априорное среднее для w = [b, w_1, w_2] matrix Sigma_w = matrix::Identity(d + 1,d + 1); Sigma_w = 10*Sigma_w; // Диагональная ковариационная матрица (дисперсия = 10) //--- Параметры сэмплирования int nsamples = 10000; int burnin = 1000; int thin = 10; vector initial(d + 1); initial.Fill(1.0); // Начальные значения [b, w_1, w_2] vector width(d + 1); width.Fill(0.5); // Ширина шага для slice sampling //--- Частотный подход: OLS-оценки параметров 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; // Доверительные интервалы 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; // Доверительные интервалы 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); } //+------------------------------------------------------------------+ //|Доверительные интервалы, частотный подход | //+------------------------------------------------------------------+ 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. Вычисление остатков и sigma^2 vector y_pred = X @ w_ols; vector residuals = y - y_pred; double sigma2 = (residuals @ residuals) / (n - d); // Оценка дисперсии sigma_hat = MathSqrt(sigma2); // Оценка стандартного отклонения шума // 2. Ковариационная матрица matrix XtX_inv = (X.Transpose() @ X).Inv(); matrix cov_matrix = sigma2 * XtX_inv; // 3. Стандартные ошибки se = cov_matrix.Diag(); se = MathSqrt(se); // 4. Критическое значение t-распределения double t_crit = 1.96; // Для 95% доверительного интервала // 5. Доверительные интервалы lower = w_ols - t_crit * se; upper = w_ols + t_crit * se; } //+------------------------------------------------------------------+ //| Доверительные интервалы, байесовский подход | //+------------------------------------------------------------------+ 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); // Массив вероятностей для квантилей double probs[] = {alpha / 2, 1 - alpha / 2}; // {0.025, 0.975} для 95%-го интервала double quantiles[]; MathQuantile(temp, probs, quantiles); lower[j] = quantiles[0]; // Квантиль alpha/2 upper[j] = quantiles[1]; // Квантиль 1-alpha/2 } } //+------------------------------------------------------------------+ //| Сохранение матрицы в CSV | //+------------------------------------------------------------------+ bool MatrixToCSV(string file_name, const matrix &m) { if(m.Rows() == 0 || m.Cols() == 0) { Print("Ошибка: Матрица пустая"); return false; } int file_handle = FileOpen(file_name, FILE_WRITE | FILE_TXT | FILE_ANSI); if(file_handle == INVALID_HANDLE) { Print("Ошибка открытия/создания файла: ", 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("Матрица успешно сохранена в файл: ", file_name); return true; } //+------------------------------------------------------------------+
В частотном подходе оценки параметров вычисляются с помощью метода наименьших квадратов (Ordinary Least Squares), а доверительные интервалы рассчитываются на основе t-распределения.
В байесовском подходе точечные оценки получаются как среднее апостериорного распределения (среднее сэмплов), а неопределённость отражают 95%-ные вероятностные интервалы, вычисленные функцией CredibleIntervals для квантилей 0.025 и 0.975.
По завершении MCMC-сэмплирования, сравним оценки для частотного и байесовского подходов. Полученные результаты демонстрируют, что мы успешно сэмплировали из апостериорного распределения параметров, и полученные оценки тесно согласуются с частотными оценками, что ожидаемо при использовании слабоинформативного априора (обычно это распределения с большой дисперсией, как в нашем случае).
Средние значения, полученные с помощью сэмплера, практически идентичны оценкам OLS. Это доказательство того, что алгоритм сошелся и правильно исследовал область плотности апостериорного распределения. Хотя все оценки немного смещены от истинных значений [1.0, 2.0, 2.0], это нормально для такой небольшой выборки объемом N = 100 и с учетом присутствия шума. Если вы запустите скрипт, ваши оценки будут немного отличаться из-за новой генерации случайных чисел.
| Параметр | Истинное значение | OLS | SS |
|---|---|---|---|
| w_0 (смещение, b) | 1 | 0.91847 | 0.91852 |
| w_1 | 2 | 1.97694 | 1.97688 |
| w_2 | 2 | 2.02333 | 2.02320 |
Байесовские вероятностные интервалы (Credible intervals) также очень близки к частотным доверительным интервалам (Confidence intervals).
| Параметр | 95% Confidence Intervals(frequency approach) | 95% Credible Intervals (bayesian approach) |
|---|---|---|
| w_0 (смещение, 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] |
Эффективность алгоритма slicesample оценивается по среднему количеству вызовов целевой функции logf(x). При данных настройках параметр neval в функции slicesample составил около 5.365.
Для регулирования этого показателя необходимо изменить величину параметра width. Следует учитывать, что если ширина width слишком мала, алгоритм будет выполнять избыточное количество вычислений функции для определения размера среза. Если же ширина слишком велика, алгоритм будет вынужден часто уменьшать интервал до подходящего значения, что также приведёт к избыточному количеству вычислений функции.
Помимо оценки neval, для полного анализа качества сэмплирования необходимо исследовать график трассы выборок и автокорреляционную функцию (ACF), а также построить гистограммы для каждого параметра модели.
import pandas as pd import numpy as np import matplotlib.pyplot as plt from statsmodels.graphics.tsaplots import plot_acf # Путь к файлу slice_file = "C:/Program Files/MetaTrader 5/MQL5/Files/SS/LR_samples.csv" # Загрузка данных slice_logpdf_data = pd.read_csv(slice_file, header=None) samples = slice_logpdf_data.values # Матрица nsamples x 3 dim = samples.shape[1] # Количество измерений (3: b, w_1, w_2) # Ожидаемые значения параметров expected_values = [1.0] + [2.0] * (dim - 1) # [1.0, 2.0, 2.0] # Статистика 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}") # Ковариация между b и w_1 (для примера) print(f"Covariance (b, w_1): {np.cov(samples[:, 0], samples[:, 1])[0, 1]:.4f}") print("=========================================") print_stats(samples, dim, expected_values) # Построение трасс и 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) # Добавляем горизонтальную линию для ожидаемого значения 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() # Построение гистограмм апостериорных распределений 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') # Добавляем вертикальную линию для ожидаемого значения 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()
График апостериорной плотности для параметра смещения w0 представлен на рис.3
Рис.3. Апостериорная плотность параметра смещения модели линейной регрессии
Пример №2. Обучение байесовской логистической регрессии
В этом примере мы применяем алгоритм slice sampling для моделирования бинарных исходов y ∈ {0,1} в задаче байесовской логистической регрессии. Мы сгенерируем синтетические данные с истинными параметрами w_true = [0.0, 1.0, −1.0], после чего попытаемся восстановить эти истинные значения с помощью частотных и байесовских методов.
В первую очередь определимся с моделью данных, то есть с правдоподобием. Правдоподобие принимает форму распределения Бернулли, параметризованного сигмоидной функцией от линейного предиктора Xw:
В качестве априора для параметров w возьмем многомерное нормальное распределение:
P(w) = N (w | 0, 10 * I)
Как и в предыдущем примере, выберем слабоинформативный априор с довольно широкой дисперсией. Напомним, что в байесовской статистике априор называется слабоинформативным, если он оказывает минимальное влияние на апостериорное распределение, позволяя правдоподобию играть основную роль в формировании выводов.
Что касается классического частотного подхода, то здесь точечные оценки можно найти с помощью алгоритма IRLS (Iteratively Reweighted Least Squares). IRLS — это итеративный метод численного решения задачи оптимизации, который хорошо подходит для модели логистической регрессии. Он минимизирует логарифмическую функцию потерь, эквивалентную максимизации правдоподобия. Обновление параметров w вычисляется по формуле:
где,
- X – матрица признаков (n×(d+1)), включая столбец единиц для смещения,
- W — диагональная матрица весов размером n×n, где элементы матрицы wii = pi (1 − pi), а pi = σ(Xiw) — вероятности, вычисленные через сигмоидную функцию,
- z = Xw + (y − p) / p*(1 − p) – так называемая рабочая переменная, где y — вектор наблюдаемых бинарных откликов, а p — вектор предсказанных вероятностей
- X'WX — информационная матрица Фишера, используемая для вычисления ковариационной матрицы и стандартных ошибок параметров.
Параметры семплера:
- nsamples 10000 - количество генерируемых сэмплов,
- burnin 1000 - период прогрева цепи,
- thin 10 - интервал прореживания,
- width 1 - начальная ширина гиперпрямоугольника.
//+------------------------------------------------------------------+ //| 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> // Структура для хранения данных struct BayesianLogisticRegression { matrix X; // Матрица признаков (n x (d+1)) vector y; // Бинарные отклики (n) vector mu_w; // Априорное среднее для w = [b, w_1, w_2] matrix Sigma_w; // Априорная ковариационная матрица для w matrix Sigma_w_inv; // Обратная ковариационная матрица double log_det_Sigma_w; // Логарифм определителя Sigma_w }; BayesianLogisticRegression model; //+------------------------------------------------------------------+ //| Сигмоидная функция | //+------------------------------------------------------------------+ double sigmoid(double z) { return 1.0 / (1.0 + MathExp(-z)); } //+------------------------------------------------------------------+ //| Логарифм плотности многомерного нормального распределения | //+------------------------------------------------------------------+ 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)); } //+------------------------------------------------------------------+ //| Целевая плотность - апостериорное распределение | //+------------------------------------------------------------------+ double LogPosterior(vector ¶ms, BayesianLogisticRegression &mdl) { int d = (int)model.X.Cols(); // Количество параметров (включая смещение) // --- 1. Логарифм правдоподобия --- int n = (int)model.X.Rows(); // Вычисление линейного предиктора: eta = X * w vector eta = model.X @ params; // Вычисление вектора вероятностей p = sigmoid(eta) vector p(n); eta.Activation(p, AF_SIGMOID); // Вычисление вектора log(p) vector log_p(n); log_p = MathLog(p); // Вычисление вектора log(1-p) vector one_minus_p_log(n); one_minus_p_log = MathLog(1.0 - p); // Вычисление логарифма правдоподобия: 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. Логарифм априорного распределения (многомерный нормальный априор) --- double log_prior = multivariate_normal_logpdf(params, mdl.mu_w, mdl.Sigma_w_inv, mdl.log_det_Sigma_w, d); // Апостериор return log_likelihood + log_prior; } // Обертка для соответствия сигнатуре LogPdfFunction double LogPost(vector ¶ms) { return LogPosterior(params, model); } //+---------------------------------------------------------------------+ //|Вычисление оценок методом IRLS(Iteratively Reweighted Least Squares) | //+---------------------------------------------------------------------+ void IRLS(matrix &X, vector &y, vector &w_irls, vector &se_irls) { int n = (int)X.Rows(); // Количество наблюдений int d = (int)X.Cols(); // Количество параметров (включая смещение) // --- Параметры алгоритма --- double convergence_threshold = 1e-6; // Порог для проверки сходимости весов int max_iterations = 20; // Максимальное количество итераций int iterations_used = 0; // Фактически использованное количество итераций w_irls.Resize(d); w_irls.Fill(0.0); // Начальное приближение весов w vector w_old = w_irls; // Вектор для хранения весов предыдущей итерации // Формула обновления параметров: 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); // Вектор вероятностей p_i = sigmoid(eta_i) vector z(n); // Вектор рабочей переменной matrix W_current(n, n); // Диагональная матрица весов W W_current.Fill(0.0); // 1. Вычисление линейного предиктора: eta = X * w vector eta = X @ w_irls; // 2. Вычисление вероятностей: p = sigmoid(eta) eta.Activation(p, AF_SIGMOID); // 3. Вычисление диагонального вектора весов: w_diag = p * (1 - p) vector w_diag = p * (1.0 - p); // 4. Создание диагональной матрицы весов W_current W_current.Diag(w_diag); // 5. Вычисление рабочей переменной z vector diff = y - p; z = eta + diff / w_diag; // Сохраняем текущую W для использования в расчете SE W_final = W_current; // 6. Обновление весов (IRLS Step) matrix XtW = X.Transpose() @ W_current; // Вычисление XtWX (информационная матрица Фишера) matrix XtWX = XtW @ X; // Решение для w: w = (XtWX)^-1 @ XtW @ z w_old = w_irls; // Сохраняем w для проверки сходимости w_irls = XtWX.Inv() @ XtW @ z; // 7. Проверка сходимости // Вычисляем L2-норму разницы весов 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."); } // --- Расчет Стандартных Ошибок (SE) --- // Вычисление XtWX (информационная матрица Фишера) с финальными весами W_final matrix XtWX_final = X.Transpose() @ W_final @ X; // Ковариационная матрица Cov(w) = (X^T W_final X)^-1 matrix cov_matrix = XtWX_final.Inv(); // стандартные ошибки se_irls = cov_matrix.Diag(); se_irls = MathSqrt(se_irls); } //+---------------------------------------------------------------------+ //| Инициализация логистчиеской регрессии и семплирование с помощью 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(); // Предварительное вычисление обратной матрицы и логарифма определителя matrix Sigma_w_inv = Sigma_w.Inv(); double det = Sigma_w.Det(); double log_det = MathLog(det); // Инициализация модели 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; // Подсчёт количества оценок logpdf и времени выполнения 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; // Время в миллисекундах Print("Average number of logpdf evaluations per sample: ", neval); Print("Slice sampling execution time: ", time_ms, " ms"); } //+------------------------------------------------------------------+ //| Доверительные интервалы, байесовский подход | //+------------------------------------------------------------------+ 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() { // Генерация синтетических данных int n = 100; // Количество наблюдений int d = 2; // Количество признаков (без параметра смещения) matrix X(n, d + 1); // Матрица признаков X с параметром смещения vector y(n); int err; vector true_w(d + 1); // w = [b, w_1, w_2] true_w[0] = 0.0; // Истинное смещение true_w[1] = 1.0; // Истинный вес w_1 true_w[2] = -1.0; // Истинный вес w_2 // Заполнение матрицы X и вектора y for(int i = 0; i < n; i++) { X[i][0] = 1.0; // Первый столбец — единицы for(int j = 1; j < d + 1; j++) { X[i][j] = MathRandomNormal(0.0, 1.0, err); // Признаки из 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); // Параметры априора vector mu_w(d + 1); mu_w.Fill(0.0); // Априорное среднее для 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; // Диагональная ковариационная матрица (дисперсия = 10) } // Параметры сэмплирования int nsamples = 10000; int burnin = 1000; int thin = 10; vector initial(d + 1); initial.Fill(0.0); // Начальные значения [b, w_1, w_2] vector width(d + 1); width.Fill(1); // Ширина шага для slice sampling matrix samples; // Частотный подход: IRLS-оценки параметров 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 ])"); // Доверительные интервалы double z_crit = 1.96; // Для 95% доверительного интервала 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], "]"); } // Сохранение w_irls в 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 ])"); // Доверительные интервалы 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); } //+------------------------------------------------------------------+
Оба подхода дают оценки, близкие к истинным значениям, хотя и с некоторыми отклонениями, из-за конечного объема данных. На рис.4 представлен график апостериорной плотности для параметра w1. Истинное значение параметра(w1=1) находится в пределах 95% байесовского доверительного интервала.
Рис. 4. Апостериорное распределение параметра w1 байесовской логистической регрессии
Для визуального подтверждения корректности результатов построим границу решения, которая в логистической регрессии определяется уравнением x^Tw = 0. Для двух признаков (x1,x2) и смещения граница является прямой.
Построим три границы:
- Истинная граница, основанная на истинных значениях параметров w_true = [0.0,1.0,−1.0].
- IRLS-граница основанная на точечных оценках.
- Байесовская граница, основанная на средних значениях апостериорного распределения.
Рис. 5. Байесовская и IRLS граница решения (совпадают)
Границы решения для IRLS и байесовского подхода практически совпадают, но отклоняются от истинной границы из-за ограниченного объёма данных (n=100) и стохастической природы генерации классов через сигмоидную функцию. Это приводит к небольшому смещению оценок параметров. Тем не менее, алгоритм выборки по уровням успешно справился с задачей, корректно идентифицируя параметры, обеспечивающие разделение классов.
Заключение
В этой статье мы исследовали метод выборки по уровням (slice sampling) — адаптивную разновидность MCMC, которая устраняет необходимость тщательной настройки гиперпараметров, например таких, как масштаб шага в алгоритме Метрополиса.
Эффективность метода проверялась с помощью байесовской линейной и логистической регрессии. Байесовские оценки (средние апостериорного распределения и вероятностные интервалы) оказались близки к результатам классических частотных методов — метода наименьших квадратов (OLS) и итеративно взвешенных наименьших квадратов (IRLS). Графики трасс выборок, автокорреляционных функций (ACF) и гистограмм апостериорной плотности подтверждают качество сэмплирования.
Таким образом, реализация семплера slicesample в MQL5, показала себя универсальным и надёжным «чёрным ящиком». Благодаря такому методу, байесовский вывод стал значительно проще и доступнее. От пользователя, по сути дела, требуется только умение закодировать целевую апостериорную плотность. Это позволяет сконцентрироваться непосредственно на задаче статистического вывода, минимизируя усилия на настройку алгоритма MCMC.
Программы используемые в статье:
| # | Имя | Тип | Описание |
|---|---|---|---|
| 1 | SS.mqh | Включаемый файл | Алгоритм выборки по уровням |
| 2 | LR.mq5 | Скрипт | Пример сэмплирования апостериорного распределения для модели байесовской линейной регрессии |
| 3 | LR_plot.py | Скрипт | Диагностика в Python |
| 4 | LogisticR.mq5 | Скрипт | Пример сэмплирования апостериорного распределения для модели байесовской логистической регрессии |
| 5 | LogisticR_plot.py | Скрипт | Диагностика в Python |
| 6 | PlotMM.mq5 | Скрипт | Анимация работы алгоритма для мультимодального одномерного распределения |
Предупреждение: все права на данные материалы принадлежат MetaQuotes Ltd. Полная или частичная перепечатка запрещена.
Данная статья написана пользователем сайта и отражает его личную точку зрения. Компания MetaQuotes Ltd не несет ответственности за достоверность представленной информации, а также за возможные последствия использования описанных решений, стратегий или рекомендаций.
От новичка до эксперта: Раскрываем скрытые уровни коррекции Фибоначчи
Нейросети в трейдинге: Агрегация движения по времени (TMA)
Нейросети в трейдинге: Агрегация движения по времени (Основные компоненты)
Возможности Мастера MQL5, которые вам нужно знать (Часть 54): Обучение с подкреплением с гибридным SAC и тензорами
- Бесплатные приложения для трейдинга
- 8 000+ сигналов для копирования
- Экономические новости для анализа финансовых рынков
Вы принимаете политику сайта и условия использования