Probenahme aus der bivariaten Verteilung mit bekannter Dichte unter Verwendung von MCMC

9

Ich habe versucht, mit Metropolis-Algorithmen in R aus einer bivariaten Dichte zu simulieren und hatte kein Glück. Die Dichte kann ausgedrückt werden als , wobei die Singh-Maddala-Verteilung istp ( y | x ) p ( x ) p ( x )p(x,y)p(y|x)p(x)p(x)

p(x)=einqxein- -1bein(1+(xb)ein)1+q

mit den Parametern , , und ist log-normal mit log-Mittelwert als Bruchteil von und log-sd eine Konstante. Um zu testen, ob meine Stichprobe die gewünschte ist, habe ich mir die Randdichte von , die . Ich habe verschiedene Metropolis-Algorithmen aus den R-Paketen MCMCpack, mcmc und dream ausprobiert. Ich habe das Einbrennen verworfen, die Ausdünnung verwendet und Proben mit einer Größe von bis zu Millionen verwendet, aber die resultierende Grenzdichte war nie die, die ich geliefert habe.q b p ( y | x ) x x p ( x )einqbp(y|x)xxp(x)

Hier ist die endgültige Ausgabe meines Codes, den ich verwendet habe:

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

Ich habe mich für ein Traumpaket entschieden, da es bis zur Konvergenz probiert. Ich habe auf drei Arten getestet, ob ich die richtigen Ergebnisse habe. Verwenden der KS-Statistik, Vergleichen von Quantilen und Schätzen der Parameter der Singh-Maddala-Verteilung mit maximaler Wahrscheinlichkeit aus der resultierenden Stichprobe:

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

Wenn ich mir die Ergebnisse dieser Vergleiche ansehe, lehnt die KS-Statistik fast immer die Nullhypothese ab, dass die Stichprobe aus der Singh-Maddala-Verteilung mit den angegebenen Parametern stammt. Die geschätzten Parameter für die maximale Wahrscheinlichkeit kommen manchmal ihren wahren Werten nahe, liegen jedoch normalerweise zu weit außerhalb der Komfortzone, um zu akzeptieren, dass das Probenahmeverfahren erfolgreich war. Ebenso wie die Quantile sind empirische Quantile nicht zu weit, sondern zu weit entfernt.

Meine Frage ist, was ich falsch mache? Meine eigenen Hypothesen:

  1. MCMC ist für diese Art der Probenahme nicht geeignet
  2. MCMC kann aus theoretischen Gründen nicht konvergieren (die Verteilungsfunktion erfüllt nicht die erforderlichen Eigenschaften, unabhängig davon, um welche es sich handelt).
  3. Ich verwende den Metropolis-Algorithmus nicht richtig
  4. Meine Verteilungstests sind nicht korrekt, da ich keine unabhängige Stichprobe habe.
mpiktas
quelle
In der Singh-Maddala-Verteilungsverbindung hat das PDF zwei Parameter - {c, k}, aber die R-Funktion dsinmadakzeptiert drei Parameter oder fehlt mir etwas.
Csgillespie
Entschuldigung, der Wikipedia-Link zitiert die falsche Formel. Auf den ersten Blick sah es in Ordnung aus, als ich die Frage verfasste. Ich habe keinen fertigen Link gefunden, also habe ich einfach die Formel in die Frage gestellt.
mpiktas

Antworten:

3

Ich denke, die Reihenfolge ist korrekt, aber die Beschriftungen, die p (x) und p (y | x) zugewiesen wurden, waren falsch. Die ursprünglichen Problemzustände p (y | x) sind logarithmisch normal und p (x) ist Singh-Maddala. So ist es

  1. Generiere ein X aus einem Singh-Maddala und

  2. Erzeugen Sie ein Y aus einer logarithmischen Normalen mit einem Mittelwert, der ein Bruchteil des erzeugten X ist.

Jan Galkowski
quelle
3

In der Tat sollten Sie MCMC nicht machen, da Ihr Problem so viel einfacher ist. Versuchen Sie diesen Algorithmus:

Schritt 1: Generieren Sie ein X aus Log Normal

Schritt 2: Halten Sie dieses X fest und generieren Sie ein Y aus der Singh Maddala.

Voilà! Probe fertig !!!

Mohit
quelle
Ich gehe davon aus, dass Sie die Schritte umgekehrt gemeint haben. Aber wenn dies so einfach ist, warum brauchen wir Gibbs-Sampling?
mpiktas
1
Nein, ich meinte die Schritte 1 und 2 in der Reihenfolge, in der ich sie geschrieben habe. Schließlich wird die Verteilung von y von X abhängig gemacht, sodass Sie vor Y ein X generieren müssen. Bei der Gibbs-Abtastung handelt es sich um eine kompliziertere Lösung für kompliziertere Probleme. Ihre, wie Sie es beschreiben, ist meiner Meinung nach ziemlich direkt.
Mohit
1
p(y|x)p(x|y)p(x)