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

9

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)=aqxa1ba(1+(xb)a)1+q

z parametrami a, q, b, i p(y|x) jest log-normal z log-średnią jako ułamek xi log-sd stała. Aby sprawdzić, czy moja próbka jest tą, której chcę, przyjrzałem się marginalnej gęstościx, 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.
mpiktas
źródło
W linku dystrybucyjnym Singh-Maddala plik pdf ma dwa parametry - {c, k}, ale funkcja R dsinmadprzyjmuje trzy parametry lub czegoś mi brakuje.
csgillespie,
Niestety, link do Wikipedii cytuje niewłaściwą formułę, na pierwszy rzut oka wyglądał dobrze, kiedy tworzyłem pytanie. Nie znalazłem gotowego linku, więc po prostu postawiłem formułę w pytaniu.
mpiktas,

Odpowiedzi:

3

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.

Jan Galkowski
źródło
3

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

Mohit
źródło
Zakładam, że miałeś na myśli odwrócenie kroków. Ale jeśli jest to tak proste, dlaczego potrzebujemy próbkowania Gibbsa?
mpiktas,
1
Nie, miałem na myśli kroki 1 i 2 w kolejności, w której napisałem. W końcu rozkład y jest określony warunkowo na X, więc musisz wygenerować X przed Y. Jeśli chodzi o próbkowanie Gibbs, jest to bardziej skomplikowane rozwiązanie przeznaczone dla bardziej skomplikowanych problemów. Twój, jak to opisujesz, jest dość prosty, IMHO.
Mohit,
1
Gdybyś wiedział, użyłbyś próbkowania Gibbs p(y|x) i p(x|y), ale nie, jeśli znasz margines p(x)
prawdopodobieństwo