Выборка Гастингса столицы
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];
