Отбор проб из двумерного распределения с известной плотностью с использованием MCMC

9

Я попытался смоделировать из двумерной плотности используя алгоритмы Метрополиса в R, и мне не повезло. Плотность может быть выражена как , где - распределение Сингх-Маддалап(Икс,Y)п(Y|Икс)п(Икс)п(Икс)

п(Икс)знак равноaQИксa-1бa(1+(Иксб)a)1+Q

с параметрами a , Q , б и п(Y|Икс) является логарифмически нормальным, с логарифмическим значением как дробью Икс , а log-sd - константой. Чтобы проверить, является ли мой образец тем, который я хочу, я посмотрел на предельную плотность Икс , которая должна быть п(Икс) . Я пробовал разные алгоритмы Metropolis из пакетов R MCMCpack, mcmc и dream. Я отказался от выжигания, использовал прореживание, использовал образцы размером до миллиона, но полученная предельная плотность никогда не была той, которую я поставил.

Вот последняя версия моего кода, который я использовал:

logvrls <- function(x,el,sdlog,a,scl,q.arg) {
    if(x[2]>0) {
         dlnorm(x[1],meanlog=el*log(x[2]),sdlog=sdlog,log=TRUE)+
         dsinmad(x[2],a=a,scale=scl,q.arg=q.arg,log=TRUE)
    }
    else -Inf    
}

a <- 1.35
q <- 3.3
scale <- 10/gamma(1 + 1/a)/gamma(q - 1/a)*  gamma(q) 

Initvrls <- function(pars,nseq,meanlog,sdlog,a,scale,q) {
    cbind(rlnorm(nseq,meanlog,sdlog),rsinmad(nseq,a,scale,q))
}

library(dream)
aa <- dream(logvrls,
        func.type="logposterior.density",
        pars=list(c(0,Inf),c(0,Inf)),
        FUN.pars=list(el=0.2,sdlog=0.2,a=a,scl=scale,q.arg=q),
        INIT=Initvrls,
        INIT.pars=list(meanlog=1,sdlog=0.1,a=a,scale=scale,q=q),
        control=list(nseq=3,thin.t=10)
        )

Я остановился на пакете мечты, так как он сэмплирует до схождения. Я проверил, имею ли я правильные результаты тремя способами. Используя статистику KS, сравнивая квантили и оценивая параметры распределения Сингх-Маддала с максимальной вероятностью из полученной выборки:

ks.test(as.numeric(aa$Seq[[2]][,2]),psinmad,a=a,scale=scale,q.arg=q)

lsinmad <- function(x,sample)
    sum(dsinmad(sample,a=x[1],scale=x[2],q.arg=x[3],log=TRUE))
 optim(c(2,20,2),lsinmad,method="BFGS",sample=aa$Seq[[1]][,2])

 qq <- eq(0.025,.975,by=0.025)   
 tst <- cbind(qq,
              sapply(aa$Seq,function(l)round(quantile(l[,2],qq),3)),
              round(qsinmad(qq,a,scale,q),3))
 colnames(tst) <- c("Quantile","S1","S2","S3","True")

 library(ggplot2)
 qplot(x=Quantile,y=value,
       data=melt(data.frame(tst),id=1), 
       colour=variable,group=variable,geom="line")

Когда я смотрю на результаты этих сравнений, статистика KS почти всегда отвергает нулевую гипотезу о том, что выборка из распределения Сингх-Маддала с предоставленными параметрами. Оцененные параметры максимального правдоподобия иногда приближаются к его истинным значениям, но обычно слишком далеко от зоны комфорта, чтобы допустить, что процедура отбора проб была успешной. То же самое для квантилей, эмпирических квантилей не слишком далеко, но слишком далеко.

У меня вопрос, что я делаю не так? Мои собственные гипотезы:

  1. MCMC не подходит для этого типа отбора проб
  2. MCMC не может сходиться по теоретическим причинам (функция распределения не удовлетворяет требуемым свойствам, какими бы они ни были)
  3. Я не правильно использую алгоритм Метрополис
  4. Мои тесты распространения не верны, так как у меня нет независимой выборки.
mpiktas
источник
В распределительной ссылке Сингх-Маддала , у pdf есть два параметра - {c, k}, но функция R dsinmadпринимает три параметра или я что-то упустил.
csgillespie
Извините, ссылка на википедию ссылается на неправильную формулу, на первый взгляд она выглядела нормально, когда я составлял вопрос. Я не нашел готовую ссылку, поэтому я просто поставил формулу в вопросе.
mpiktas

Ответы:

3

Я думаю, что порядок правильный, но метки, назначенные p (x) и p (y | x), были неправильными. Исходная проблема состояний p (y | x) лог-нормальна, а p (x) Сингх-Маддала. Так что это

  1. Создать X из Сингх-Маддала и

  2. генерировать Y из логарифмического нормального, имеющего среднее значение, которое является частью сгенерированного X.

Ян Галковски
источник
3

На самом деле, вы не должны делать MCMC, так как ваша проблема намного проще. Попробуйте этот алгоритм:

Шаг 1: сгенерируйте X из логарифма

Шаг 2: Сохраняя этот X фиксированным, сгенерируйте Y из Сингх Маддала.

Вуаля! Образец готов!

Мохит
источник
Я предполагаю, что вы имели в виду шаги, обратные. Но если это так просто, зачем нам выборка Гиббса?
mpiktas
1
Нет, я имел в виду шаги 1 и 2 в порядке, который я написал. В конце концов, распределение y задано условно на X, поэтому вы должны сгенерировать X перед Y. Что касается выборки Гиббса, то это более сложное решение, предназначенное для более сложных задач. Твоя, как ты это описал, довольно прямолинейна, ИМХО.
Мохит
1
п(Y|Икс)п(Икс|Y)п(Икс)