Exponenta Banner
exponenta event banner

Байесовская линейная регрессия с использованием гамильтонова Монте-Карло

Открыть сценарий

В этом примере показано, как выполнить байесовский вывод для модели линейной регрессии с использованием гамильтонова образца Монте-Карло (HMC).

В байесовском выводе параметров целью является анализ статистических моделей с включением предшествующего знания параметров модели. Апостериорное распределение свободных параметров$\theta$ объединяет функцию правдоподобия$P(y | \theta)$ с предыдущим распределением, $P(\theta)$используя теорему Байеса:

$$ P(\theta | y) = \frac{P(y | \theta) P(\theta)}{P(y)}.$$

Обычно лучшим способом суммирования заднего распределения является получение образцов из этого распределения с использованием методов Монте-Карло. Используя эти выборки, можно оценить предельные апостериорные распределения и производную статистику, такую как заднее среднее, среднее и стандартное отклонение. HMC - это градиентный образец Маркова цепи Монте-Карло, который может быть более эффективным, чем стандартные образцы, особенно для среднеразмерных и высокоразмерных проблем.

Модель линейной регрессии

Проанализируйте модель линейной регрессии с пересечением, $\alpha$линейными коэффициентами ($\beta$вектором столбца) и дисперсией шума$\sigma^2$ распределения данных в качестве свободных параметров. Предположим, что каждая точка данных имеет независимое гауссово распределение:

$$ y_i|\theta \sim \mathcal{N}(\mu_i(\theta), \sigma^2). $$

Моделирование среднего$\mu_i$ гауссова распределения как функции предикторов$x_i$ и параметров модели как

$$\mu_i = \alpha + x_i^T \beta.$$

В байесовском анализе необходимо также назначить всем свободным параметрам предыдущие распределения. Назначьте независимые пауссовы приоры на пересечении и линейные коэффициенты:

$\alpha \sim \mathcal{N}(\alpha_0, \sigma_\alpha^2)$,

$\beta_i \sim \mathcal{N}(\beta_{0}, \sigma_{\beta}^2)$.

Для использования HMC все переменные выборки должны быть неограниченными, что означает, что задняя плотность и ее градиент должны быть хорошо определены для всех реальных значений параметров. Если параметр ограничен интервалом, его необходимо преобразовать в неограниченный. Чтобы сохранить вероятность, необходимо умножить предыдущее распределение на соответствующий коэффициент Якобиана. Также этот фактор учитывается при вычислении градиента заднего.

Дисперсия шума является (в квадрате) параметром масштаба, который может быть только положительным. Тогда проще и естественнее рассматривать его логарифм как свободный параметр, который неограничен. Назначьте нормаль перед логарифмом дисперсии шума:

$\log \sigma^2 \sim \mathcal{N}(\kappa, \omega^2)$.

Запишите логарифм задней плотности свободных параметров$\theta = (\alpha; \beta; \log \sigma^2)$ как

$$\log P(\theta | y) = \mathrm{const.} + \log P(y | \theta) + \log
P(\theta).$$

Игнорируйте постоянный член и вызывайте сумму последних двух членов. $g(\theta)$Чтобы использовать HMC, создайте дескриптор функции, который вычисляет$g(\theta)$ и его градиент${\partial g}/{\partial \theta}$ для любого значения. $\theta$Функции, используемые для вычисления$g(\theta)$, находятся в конце сценария.

Создать набор данных

Определение истинных значений параметров для пересечения, линейных коэффициентов Betaи стандартное отклонение шума. Знание истинных значений параметров позволяет сравнивать с выходными данными дискретизатора HMC. На реакцию влияет только первый предиктор.

NumPredictors = 2;

trueIntercept = 2;
trueBeta = [3;0];
trueNoiseSigma = 1;

Используйте эти значения параметров для создания нормально распределенного набора данных выборки при случайных значениях двух предикторов.

NumData = 100;
rng('default') %For reproducibility
X = rand(NumData,NumPredictors);
mu = X*trueBeta + trueIntercept;
y = normrnd(mu,trueNoiseSigma);

Определение задней плотности вероятности

Выберите средства и стандартные отклонения гауссовых приоров.

InterceptPriorMean = 0;
InterceptPriorSigma = 10;
BetaPriorMean = 0;
BetaPriorSigma = 10;
LogNoiseVarianceMean = 0;
LogNoiseVarianceSigma = 2;

Сохранение функции logPosterior на пути MATLAB ®, который возвращает логарифм произведения предшествующего и правдоподобия, и градиент этого логарифма. logPosterior определяется в конце этого примера. Затем вызовите функцию с аргументами, чтобы определить logpdf входной аргумент для hmcSampler функция.

logpdf = @(Parameters)logPosterior(Parameters,X,y, ...
    InterceptPriorMean,InterceptPriorSigma, ...
    BetaPriorMean,BetaPriorSigma, ...
    LogNoiseVarianceMean,LogNoiseVarianceSigma);

Создание образца HMC

Определите начальную точку для начала выборки, а затем вызовите hmcSampler для создания гамильтонова пробоотборника в виде HamiltonianSampler объект. Отображение свойств пробоотборника.

Intercept = randn;
Beta = randn(NumPredictors,1);
LogNoiseVariance = randn;
startpoint = [Intercept;Beta;LogNoiseVariance];
smp = hmcSampler(logpdf,startpoint,'NumSteps',50);

smp

smp = 

  HamiltonianSampler with properties:

                  StepSize: 0.1000
                  NumSteps: 50
                MassVector: [4x1 double]
              JitterMethod: 'jitter-both'
      StepSizeTuningMethod: 'dual-averaging'
    MassVectorTuningMethod: 'iterative-sampling'
                    LogPDF: [function_handle]
             VariableNames: {4x1 cell}
                StartPoint: [4x1 double]

Оценка точки MAP

Оценить точку MAP (maximum-a-posteriori) задней плотности. Можно начать выборку из любой точки, но часто эффективнее оценить точку MAP, а затем использовать ее в качестве начальной точки для настройки образца и образцов чертежа. Оцените и отобразите точку MAP. Вы можете показать больше информации во время оптимизации, установив 'VerbosityLevel' значение 1.

[MAPpars,fitInfo] = estimateMAP(smp,'VerbosityLevel',0);
MAPIntercept = MAPpars(1)
MAPBeta = MAPpars(2:end-1)
MAPLogNoiseVariance = MAPpars(end)

MAPIntercept =

    2.3857


MAPBeta =

    2.5495
   -0.4508


MAPLogNoiseVariance =

   -0.1007

Чтобы проверить сходимость оптимизации к локальному оптимуму, постройте график fitInfo.Objective поле. Это поле содержит значения отрицательной логарифмической плотности при каждой итерации оптимизации функции. Конечные значения все похожи, поэтому оптимизация сошлась.

plot(fitInfo.Iteration,fitInfo.Objective,'ro-');
xlabel('Iteration');
ylabel('Negative log density');

Настройка сэмплера

Важно выбрать хорошие значения для параметров пробоотборника, чтобы получить эффективную выборку. Лучший способ найти хорошие значения - это автоматически настроить MassVector, StepSize, и NumSteps параметры с использованием tuneSampler способ. Используйте метод для:

1. Настройка MassVector пробоотборника.

2. Мелодия StepSize и NumSteps для фиксированной длины моделирования для достижения определенного коэффициента приемки. Целевой коэффициент принятия по умолчанию 0,65 в большинстве случаев является хорошим.

Начните настройку в расчетной точке MAP для более эффективной настройки.

[smp,tuneinfo] = tuneSampler(smp,'Start',MAPpars);

Постройте график изменения размера шага во время настройки, чтобы гарантировать сходимость настройки размера шага. Просмотрите достигнутый коэффициент приемки.

figure;
plot(tuneinfo.StepSizeTuningInfo.StepSizeProfile);
xlabel('Iteration');
ylabel('Step size');

accratio = tuneinfo.StepSizeTuningInfo.AcceptanceRatio

accratio =

    0.6400

Рисование образцов

Извлеките образцы из задней плотности, используя несколько независимых цепей. Выберите различные начальные точки для цепей, случайным образом распределенные вокруг расчетной точки MAP. Укажите количество выгорающих образцов, которые будут отброшены с начала цепочки Маркова, и количество образцов, которые будут сгенерированы после выгорания.

Установите 'VerbosityLevel' значение для печати подробных данных во время выборки для первой цепочки.

NumChains = 4;
chains = cell(NumChains,1);
Burnin = 500;
NumSamples = 1000;
for c = 1:NumChains
    if (c == 1)
        level = 1;
    else
        level = 0;
    end
    chains{c} = drawSamples(smp,'Start',MAPpars + randn(size(MAPpars)), ...
        'Burnin',Burnin,'NumSamples',NumSamples, ...
        'VerbosityLevel',level,'NumPrint',300);
end


|==================================================================================|
|   ITER   |    LOG PDF    |  STEP SIZE  |  NUM STEPS  |  ACC RATIO  |  DIVERGENT  |
|==================================================================================|
|      300 | -1.484164e+02 |   2.532e-01 |          12 |   9.500e-01 |           0 |
|      600 | -1.492436e+02 |   2.128e-02 |           4 |   9.450e-01 |           0 |
|      900 | -1.509753e+02 |   2.171e-01 |           5 |   9.444e-01 |           0 |
|     1200 | -1.493455e+02 |   1.128e-01 |          16 |   9.358e-01 |           0 |
|     1500 | -1.489602e+02 |   2.532e-01 |          12 |   9.373e-01 |           0 |

Проверка диагностики сходимости

Используйте diagnostics способ вычисления стандартной диагностики MCMC. Для каждого параметра выборки метод использует все цепочки для вычисления этих статистических данных:

  • Задняя средняя оценка (Mean)

  • Оценка стандартной ошибки Монте-Карло (MCSE), которое является стандартным отклонением задней средней оценки

  • Оценка заднего стандартного отклонения (SD)

  • Оценки 5-го и 95-го квантилей краевого заднего распределения (Q5 и Q95)

  • Эффективный размер выборки для задней средней оценки (ESS)

  • Статистика сходимости Гельмана-Рубина (RHat). Как правило, значения RHat меньше, чем 1.1 интерпретируются как признак сходимости цепи к требуемому распределению. Если RHat для любой переменной больше, чем 1.1, затем попробуйте нарисовать больше образцов с помощью drawSamples способ.

Просмотрите диагностическую таблицу и истинные значения параметров выборки, определенных в начале примера. Поскольку предыдущее распределение является неинформативным для этого набора данных, истинные значения находятся между (или вблизи) 5-м и 95-м квантилями.

diags = diagnostics(smp,chains)
truePars = [trueIntercept;trueBeta;log(trueNoiseSigma^2)]

diags =

  4x8 table

     Name       Mean         MCSE         SD          Q5         Q95       ESS       RHat 
    ______    _________    _________    _______    ________    _______    ______    ______

    {'x1'}       2.3805    0.0053354    0.28028      1.8966     2.8335    2759.7    1.0005
    {'x2'}       2.5544    0.0062264    0.33478      1.9989     3.1026    2890.9    1.0004
    {'x3'}     -0.44516    0.0064468     0.3415     -1.0247    0.10198      2806         1
    {'x4'}    -0.062826    0.0028183    0.14112    -0.28345    0.17997    2507.5    1.0004


truePars =

     2
     3
     0
     0

Визуализация образцов

Исследовать такие проблемы, как сходимость и смешение, чтобы определить, представляют ли извлеченные выборки разумный набор случайных реализаций из целевого распределения. Чтобы проверить выходные данные, постройте графики трассировки выборок с использованием первой цепочки.

drawSamples способ отбрасывает загоревшие образцы с начала цепи Маркова для уменьшения влияния начальной точки отбора проб. Кроме того, графики трассировки выглядят как высокочастотный шум без какой-либо видимой корреляции между выборками. Это поведение указывает на то, что цепь хорошо перемешана.

figure;
subplot(2,2,1);
plot(chains{1}(:,1));
title(sprintf('Intercept, Chain 1'));
for p = 2:1+NumPredictors
    subplot(2,2,p);
    plot(chains{1}(:,p));
    title(sprintf('Beta(%d), Chain 1',p-1));
end
subplot(2,2,4);
plot(chains{1}(:,end));
title(sprintf('LogNoiseVariance, Chain 1'));

Объедините цепочки в одну матрицу и создайте графики рассеяния и гистограммы для визуализации 1-D и 2-D краевых задних распределений.

concatenatedSamples = vertcat(chains{:});
figure;
plotmatrix(concatenatedSamples);
title('All Chains Combined');

Функции для вычисления апостериорного распределения

logPosterior функция возвращает логарифм произведения нормального правдоподобия и нормы перед линейной моделью. Входной аргумент Parameter имеет формат [Intercept;Beta;LogNoiseVariance]. X и Y содержат значения предикторов и ответа соответственно.

normalPrior функция возвращает логарифм многомерной нормальной плотности вероятности со средним значением Mu и стандартные отклонения Sigma, указанные как скаляры или векторы столбцов той же длины, что и P. Вторым выходным аргументом является соответствующий градиент.

function [logpdf, gradlogpdf] = logPosterior(Parameters,X,Y, ...
    InterceptPriorMean,InterceptPriorSigma, ...
    BetaPriorMean,BetaPriorSigma, ...
    LogNoiseVarianceMean,LogNoiseVarianceSigma)


% Unpack the parameter vector
Intercept        = Parameters(1);
Beta             = Parameters(2:end-1);
LogNoiseVariance = Parameters(end);
% Compute the log likelihood and its gradient
Sigma                   = sqrt(exp(LogNoiseVariance));
Mu                      = X*Beta + Intercept;
Z                       = (Y - Mu)/Sigma;
loglik                  = sum(-log(Sigma) - .5*log(2*pi) - .5*Z.^2);
gradIntercept1          = sum(Z/Sigma);
gradBeta1               = X'*Z/Sigma;
gradLogNoiseVariance1	= sum(-.5 + .5*(Z.^2));
% Compute log priors and gradients
[LPIntercept, gradIntercept2]           = normalPrior(Intercept,InterceptPriorMean,InterceptPriorSigma);
[LPBeta, gradBeta2]                     = normalPrior(Beta,BetaPriorMean,BetaPriorSigma);
[LPLogNoiseVar, gradLogNoiseVariance2]  = normalPrior(LogNoiseVariance,LogNoiseVarianceMean,LogNoiseVarianceSigma);
logprior                                = LPIntercept + LPBeta + LPLogNoiseVar;
% Return the log posterior and its gradient
logpdf               = loglik + logprior;
gradIntercept        = gradIntercept1 + gradIntercept2;
gradBeta             = gradBeta1 + gradBeta2;
gradLogNoiseVariance = gradLogNoiseVariance1 + gradLogNoiseVariance2;
gradlogpdf           = [gradIntercept;gradBeta;gradLogNoiseVariance];
end

function [logpdf,gradlogpdf] = normalPrior(P,Mu,Sigma)
Z          = (P - Mu)./Sigma;
logpdf     = sum(-log(Sigma) - .5*log(2*pi) - .5*(Z.^2));
gradlogpdf = -Z./Sigma;
end

См. также

Функции

Классы

Документация по инструментам для статистического и машинного обучения

Поддержка