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 распределения использование алгоритма Metropolis-Hastings.
start - вектор-строка, содержащая начальное значение марковской цепи, nsamples - целое число, задающее количество сгенерируемых выборок и pdf, proppdf, и proprnd являются указателями на функцию, созданными с помощью @. proppdf определяет плотность распределения предложений и proprnd определяет генератор случайных чисел для распределения предложения. pdf и proprnd взять один аргумент в качестве входа с тем же типом и размером, что и start. proppdf принимает два аргумента в качестве входов с тем же типом и размером, что и start.
smpl - вектор-столбец или матрица, содержащая выборки. Если функция плотности журнала является предпочтительной, 'pdf' и 'proppdf' можно заменить на 'logpdf' и 'logproppdf'. Функции плотности, используемые в алгоритме Metropolis-Hastings, не обязательно нормализованы.
Предложение распределения 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 реализует Независимую выборку Митрополиса-Гастингса.
smpl = mhsample(...,'symmetric',sym) рисует nsamples случайные выборки из целевого стационарного pdf распределения использование алгоритма Metropolis-Hastings. sym является логическим значением, которое указывает, является ли распределение предложений симметричным. Значение по умолчанию false, что соответствует асимметричному распределению предложений. Если sym верно, например, распределение предложений симметрично, proppdf и logproppdf являются необязательными.
smpl = mhsample(...,'burnin',K) генерирует марковскую цепь со значениями между начальной точкой и kth точка опущена в сгенерированной последовательности. Значения за пределами kth сохраняются точки. k является неотрицательным целым числом со значением по умолчанию 0.
smpl = mhsample(...,'thin',m) генерирует марковскую цепь с m-1 вне m значения, опущенные в сгенерированной последовательности. m является положительным целым числом со значением по умолчанию 1.
smpl = mhsample(...,'nchain',n) генерирует n Марковские цепи с помощью алгоритма Metropolis-Hastings. 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)
Случайная прогулка Metropolis-Hastings Выборка
Используйте выборку 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];
