mhsample

Выборка Гастингса столицы

Синтаксис

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 реализует Случайный Обход выборка Гастингса Столицы. Если proppdf или logproppdf удовлетворяют q (x, y) = q (x), то есть, распределение предложения независимо от текущих значений, mhsample реализует Независимую выборку Гастингса Столицы.

smpl = mhsample(...,'symmetric',sym) чертит случайные выборки nsamples от целевого стационарного распределения pdf с помощью алгоритма Гастингса Столицы. sym является логическим значением, которое указывает, симметрично ли распределение предложения. Значение по умолчанию является ложным, который соответствует асимметричному распределению предложения. Если sym верен, например, распределение предложения симметрично, proppdf и logproppdf являются дополнительными.

smpl = mhsample(...,'burnin',K) генерирует Цепь Маркова со значениями между отправной точкой и k th точка, не использованная в сгенерированной последовательности. Значения вне k th точка сохранены. k является неотрицательным целым числом со значением по умолчанию 0.

smpl = mhsample(...,'thin',m) генерирует Цепь Маркова с m-1 из значений m, не использованных в сгенерированной последовательности. m является положительным целым числом со значением по умолчанию 1.

smpl = mhsample(...,'nchain',n) генерирует Цепи Маркова n с помощью алгоритма Гастингса Столицы. n является положительным целым числом со значением по умолчанию 1. smpl является матрицей, содержащей выборки. Последняя размерность содержит индексы для отдельных цепочек.

[smpl,accept] = mhsample(...) также возвращает accept, пропускную способность предложенного распределения. accept является скаляром, если одна цепочка сгенерирована и является вектором, если кратные цепи сгенерированы.

Примеры

свернуть все

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

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)

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

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];

Представленный в R2006a