Pobieranie próbek z rozkładu dwuwymiarowego o znanej gęstości za pomocą MCMC

Próbowałem symulować z dwuwymiarowej gęstości $p(x,y)$używając algorytmów Metropolis w R i nie miał szczęścia. Gęstość można wyrazić jako $p(y|x)p(x)$, gdzie $p(x)$ jest dystrybucją Singh-Maddala

$p(x)=\dfrac{aq x^{a-1}}{b^a (1 + (\frac{x}{b})^a)^{1+q}}$

z parametrami $a$, $q$, $b$, i $p(y|x)$ jest log-normal z log-średnią jako ułamek $x$i log-sd stała. Aby sprawdzić, czy moja próbka jest tą, której chcę, przyjrzałem się marginalnej gęstości$x$, które powinno być $p(x)$. Próbowałem różnych algorytmów Metropolis z pakietów R MCMCpack, mcmc i dream. Odrzuciłem wypalenie, użyłem rozcieńczania, użyłem próbek o wielkości do miliona, ale wynikowa gęstość krańcowa nigdy nie była tą, którą podałem.

Oto ostatnia edycja mojego kodu, którego użyłem:

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)
        )

Zdecydowałem się na pakiet marzeń, ponieważ pobiera próbki do zbieżności. Testowałem, czy mam prawidłowe wyniki na trzy sposoby. Wykorzystanie statystyki KS, porównanie kwantyli i oszacowanie parametrów rozkładu Singh-Maddala z maksymalnym prawdopodobieństwem z otrzymanej próbki:

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")

Kiedy patrzę na wyniki tych porównań, statystyka KS prawie zawsze odrzuca hipotezę zerową, że próbka pochodzi z rozkładu Singha-Maddali z podanymi parametrami. Szacowane parametry maksymalnego prawdopodobieństwa czasami zbliżają się do swoich prawdziwych wartości, ale zwykle znajdują się zbyt daleko poza strefą komfortu, aby zaakceptować, że procedura pobierania próbek zakończyła się powodzeniem. To samo dotyczy kwantyli, kwantyle empiryczne nie są zbyt daleko, ale zbyt daleko.

Moje pytanie brzmi: co robię źle? Moje własne hipotezy:

1. MCMC nie jest odpowiednie dla tego rodzaju pobierania próbek
2. MCMC nie może zbiegać się z przyczyn teoretycznych (funkcja rozkładu nie spełnia wymaganych właściwości, niezależnie od tego, jakie są)
3. Nie używam poprawnie algorytmu Metropolis
4. Moje testy dystrybucji są nieprawidłowe, ponieważ nie mam niezależnej próbki.

Myślę, że kolejność jest prawidłowa, ale etykiety przypisane do p (x) i p (y | x) były nieprawidłowe. Pierwotny problem mówi, że p (y | x) jest log-normalny, a p (x) to Singh-Maddala. Więc jest to

1. Wygeneruj X z Singh-Maddala i

2. generuje Y z logarytmu normalnego o średniej, która jest ułamkiem wygenerowanego X.

W rzeczywistości nie powinieneś robić MCMC, ponieważ twój problem jest o wiele prostszy. Wypróbuj ten algorytm:

Krok 1: Wygeneruj X z Log Normal

Krok 2: Utrzymując ten X na stałym poziomie, wygeneruj Y z Singh Maddala.

Voilà! Próbka gotowa !!!