Образец Метрополиса-Гастингса
smpl = mhsample(start,nsamples,'pdf',pdf,'proppdf',proppdf, 'proprnd',proprnd)
smpl = mhsample(...,'symmetric',sym)
smpl = mhsample(...,'burnin',K)
smpl = mhsample(...,'thin',m)
smpl = mhsample(...,'nchain',n)
[smpl,accept] = mhsample(...)
smpl = mhsample(start,nsamples,'pdf',pdf,'proppdf',proppdf, 'proprnd',proprnd) тянет nsamples случайные выборки из целевого стационарного распределения pdf используя алгоритм Метрополиса-Гастингса.
start - вектор строки, содержащий начальное значение цепи Маркова, nsamples - целое число, указывающее количество генерируемых выборок, и pdf, proppdf, и proprnd являются дескрипторами функций, созданными с помощью @. proppdf определяет плотность распределения предложения, и proprnd определяет генератор случайных чисел для распределения предложений. pdf и proprnd взять один аргумент в качестве входных данных с тем же типом и размером, что и start. proppdf принимает два аргумента в качестве входных данных с тем же типом и размером, что и start.
smpl - вектор столбца или матрица, содержащая образцы. Если предпочтительна функция логарифмической плотности, 'pdf' и 'proppdf' может быть заменен на 'logpdf' и 'logproppdf'. Функции плотности, используемые в алгоритме Метрополиса-Гастингса, не обязательно нормализуются.
Предлагаемое распределение q (x, y) дает плотность вероятности для выбора x в качестве следующей точки, когда y является текущей точкой. Иногда записывается как q (x 'y).
Если proppdf или logproppdf удовлетворяет q (x, y) = q (y, x), то есть предложенное распределение симметрично,mhsample реализует выборку Random Walk Metropolis-Hastings. Если proppdf или logproppdf удовлетворяет q (x, y) = q (x), то есть распределение предложения не зависит от текущих значений,mhsample реализует независимую выборку Metropolis-Hastings.
smpl = mhsample(...,'symmetric',sym) тянет nsamples случайные выборки из целевого стационарного распределения pdf используя алгоритм Метрополиса-Гастингса. sym - логическое значение, указывающее, является ли распределение предложений симметричным. Значением по умолчанию является false, что соответствует асимметричному распределению предложений. Если sym верно, например, распределение предложений симметрично, proppdf и logproppdf являются необязательными.
smpl = mhsample(...,'burnin',K) создает цепочку Маркова со значениями между начальной точкой и kВ сформированной последовательности пропущена третья точка. Значения за пределами kСохраняют десятую точку. k - неотрицательное целое число со значением по умолчанию 0.
smpl = mhsample(...,'thin',m) генерирует цепочку Маркова с m-1 вне m значения опущены в сгенерированной последовательности. m является положительным целым числом со значением по умолчанию 1.
smpl = mhsample(...,'nchain',n) производит n цепи Маркова с использованием алгоритма Метрополиса-Гастингса. n является положительным целым числом со значением по умолчанию 1. smpl - матрица, содержащая образцы. Последнее измерение содержит индексы для отдельных цепей.
[smpl,accept] = mhsample(...) также возвращает accept, коэффициент принятия предлагаемого распределения. accept является скаляром, если генерируется одна цепь, и вектором, если генерируется несколько цепей.
Используйте независимую выборку Metropolis-Hastings для оценки момента второго порядка гамма-распределения.
rng default; % For reproducibility alpha = 2.43; beta = 1; pdf = @(x)gampdf(x,alpha,beta); % Target distribution proppdf = @(x,y)gampdf(x,floor(alpha),floor(alpha)/alpha); proprnd = @(x)sum(... exprnd(floor(alpha)/alpha,floor(alpha),1)); nsamples = 5000; smpl = mhsample(1,nsamples,'pdf',pdf,'proprnd',proprnd,... 'proppdf',proppdf);
Постройте график результатов.
xxhat = cumsum(smpl.^2)./(1:nsamples)'; figure; plot(1:nsamples,xxhat)
Используйте выборку Random Walk Metropolis-Hastings для генерации данных выборки из стандартного нормального распределения.
rng default % For reproducibility delta = .5; pdf = @(x) normpdf(x); proppdf = @(x,y) unifpdf(y-x,-delta,delta); proprnd = @(x) x + rand*2*delta - delta; nsamples = 15000; x = mhsample(1,nsamples,'pdf',pdf,'proprnd',proprnd,'symmetric',1);
Постройте график данных образца.
figure; h = histfit(x,50); h(1).FaceColor = [.8 .8 1];
