В этом примере показано, как выполнить Байесов вывод на модели линейной регрессии использование сэмплера Гамильтонова Монте-Карло (HMC).
В Байесовом выводе параметра цель состоит в том, чтобы анализировать статистические модели с объединением предварительных знаний параметров модели. Апостериорное распределение свободных параметров комбинирует функцию правдоподобия с предшествующим распределением, с помощью теоремы Бейеса:
Обычно, лучший способ обобщить апостериорное распределение состоит в том, чтобы получить выборки из того распределения с помощью Методов Монте-Карло. Используя эти выборки, можно оценить крайние апостериорные распределения и выведенную статистику, такие как следующее среднее значение, медиана и стандартное отклонение. HMC является основанной на градиенте Цепью Маркова сэмплер Монте-Карло, который может быть более эффективным, чем стандартные сэмплеры, специально для средних размерных и высоко-размерных проблем.
Анализируйте модель линейной регрессии с прерыванием, линейные коэффициенты (вектор-столбец), и шумовое отклонение распределения данных как свободные параметры. Примите, что каждая точка данных имеет независимое Распределение Гаусса:
Смоделируйте среднее значение Распределения Гаусса как функция предикторов и параметров модели как
В Байесовом анализе также необходимо присвоить предшествующие распределения всем свободным параметрам. Присвойте независимое Гауссово уголовное прошлое на прерывании и линейных коэффициентах:
,
.
Чтобы использовать HMC, все переменные выборки должны быть неограничены, означая, что следующая плотность и ее градиент должны быть четко определены для всех действительных значений параметров. Если у вас есть параметр, который ограничивается к интервалу, то необходимо преобразовать этот параметр в неограниченный. Чтобы сохранить вероятность, необходимо умножить предшествующее распределение на соответствующий якобиевский фактор. Кроме того, примите этот фактор во внимание при вычислении градиента следующего.
Шумовое отклонение является масштабным коэффициентом (в квадрате), который может только быть положительным. Это затем может быть легче и более естественным рассмотреть свой логарифм как свободный параметр, который неограничен. Присвойте нормальное до логарифма шумового отклонения:
.
Запишите логарифм следующей плотности свободных параметров как
Проигнорируйте постоянный термин и вызовите сумму последних двух сроков. Чтобы использовать HMC, создайте указатель на функцию, который оценивает и его градиент для любого значения. Функции раньше вычисляли, расположены в конце скрипта.
Задайте истинные значения параметров для прерывания, линейные коэффициенты 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);
Задайте начальную точку, чтобы начать произвести от, и затем вызывать 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, и затем использовать его в качестве начальной точки для настройки сэмплера и рисования выборок. Оцените и отобразите точку 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
Неродной размер
, и 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.483198e+02 | 2.954e-01 | 10 | 9.167e-01 | 0 | | 600 | -1.490021e+02 | 2.483e-02 | 3 | 9.100e-01 | 0 | | 900 | -1.509799e+02 | 2.533e-01 | 4 | 9.133e-01 | 0 | | 1200 | -1.496693e+02 | 1.316e-01 | 14 | 9.050e-01 | 0 | | 1500 | -1.486535e+02 | 2.954e-01 | 10 | 9.107e-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.3852 0.0055959 0.27807 1.9185 2.8472 2469.3 1.0006 {'x2'} 2.5495 0.0063275 0.33239 1.9995 3.0968 2759.6 1.0001 {'x3'} -0.44887 0.00663 0.34165 -1.0257 0.10043 2655.4 1.0001 {'x4'} -0.060602 0.002823 0.14217 -0.28334 0.18381 2536.3 1.0003 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'));
Объедините цепи в одну матрицу и создайте графики рассеивания и гистограммы, чтобы визуализировать 1D и 2D крайние апостериорные распределения.
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