Logarithmische Gleichverteilung

27.02.2012, 19:34

vogelstrauss

Logarithmische Gleichverteilung

Hallo,

ich habe ein Problem beim samplen von einer logarithmischen Gleichverteilung in R mit Hilfe des Metropolis-Hasting Algorithmus.
Bei der nichtlogarithmischen Gleichverteilung im Intervall [a,b] mit der Dichtefunktion

$\begin{eqnarray*} pdf=\left\{ \begin{array}{rl} 1/(b-a), x \in [a,b]\\ 0,sonst \end{array}\right. \end{eqnarray*}$

klappt es (siehe Code unten), ich sample x 10000 mal und bekomme einen Mittelwert von ca 0.5 bei a=0 und b=1.
Bei der logarithmischen Glecihverteilung mit Dichtefunktion

$\begin{eqnarray*} pdf=\left\{ \begin{array}{rl} 1/(x*(ln(b)-ln(a))), x \in [a,b]\\ 0,sonst \end{array}\right. \end{eqnarray*}$

sample ich x 10000 mal, aber der Mittelwert von exp(aller x-Werte) ist ungefaehr 1.
Wie kann das sein? Warum kann ich nicht meine Verteilung als logarithmierte Werte samplen und dann mit der e-Funktion wieder auf eine nicht-logarithmische Skala bringen?
Irgendwas wichtiges bei logarithmischen Verteilungen scheint mir noch nicht klar zu sein...

Hier mein Code in R:

code:

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:

require('mcmc')

cost_uni <- function(p){
  a = 0.0001
  b = 1
  cost <- -Inf ##die metrop() funktion will -Inf wenn W-keit=0
  if (a<=p && p<=b)
    cost <- 1./(b-a)
  return(cost)
}

cost_loguni <- function(p){
  a = 0.0001
  b = 1
  cost <- -Inf
  if (a<=p && p<=b)
    cost <- 1./(p*((log(b)) - log(a)))
  return(cost)
}

mean(metrop(cost_uni, 1, 10000, scale=1, debug=TRUE)$batch) ## ~0.5
mean(exp(metrop(cost_loguni, 1, 10000, scale=1, debug=TRUE)$batch)) ## ~1, aber warum?

Ich hoffe, mir kann jemand helfen, ich bin schon sehr lange am Raetseln!!

Vielen Dank und viele Gruesse,

vogelstrauss

27.02.2012, 19:59

HAL 9000

Auf diesen Beitrag antworten »

Zitat:

Original von vogelstrauss

code:
1:	cost <- -Inf ##die metrop() funktion will -Inf wenn W-keit=0

Kann es sein, dass deine Funktionen gar nicht die Dichten, sondern die Logarithmen der Dichten zurückliefern müssen? Aus welchen Gründen sonst soll hier $\begin{eqnarray*} -\infty \end{eqnarray*}$ statt 0 stehen? verwirrt

Bei der Gleichverteilung fällt der Unterschied nicht auf, da MCMC die Dichten sowieso nur bis auf einen Vorfaktor genau benötigt - bei der anderen Verteilung schon!

28.02.2012, 12:33

vogelstrauss

Auf diesen Beitrag antworten »

Ich glaube nicht, dass die Funktion die Logarithmen der Dichten liefern muss...
Hier ein Auszug aus der Dokumentation der metrop-Funktion:

code:
1: 2: 3: 4:	an R function that evaluates the log unnormalized probability density of the desired equilibrium distribution of the Markov chain.

Deswegen scheint auch das Beispiel fuer die non-logarithmische Skala mir den erwarteten Mittelwert zu geben (~0.5), aber ich verstehe immer noch nicht, wie ich es anstellen kann, im Logspace zu samplen damit mean(exp(alle_samples)) wieder 0.5 ist...

Danke fuer die Antwort!

28.02.2012, 13:42

HAL 9000

Auf diesen Beitrag antworten »

Ohne weiteren Kommentar

Zitat:

an R function that evaluates the log unnormalized probability density of the desired equilibrium distribution of the Markov chain.

-------------------------------------

Nachdem diese Sache geklärt ist, dann doch noch ein Wort zu deiner "logarithmischen Gleichverteilung", dieser Begriff ist ja nicht gerade allgemein gebräuchlich:

Du sprichst davon, dass $\begin{eqnarray*} X \end{eqnarray*}$ auf $\begin{eqnarray*} [a,b] \end{eqnarray*}$ logarithmisch gleichverteilt ist, wenn $\begin{eqnarray*} \ln(X) \end{eqnarray*}$ auf $\begin{eqnarray*} [\ln(a),\ln(b)] \end{eqnarray*}$ stetig gleichverteilt ist - oder?

In dem Fall ist aber

$\begin{eqnarray*} E(X) = \int\limits_a^b x\cdot f_X(x) ~ \dd x = \frac{b-a}{\ln(b)-\ln(a)} \end{eqnarray*}$ ,

was im Fall $\begin{eqnarray*} a = 0.0001,\;b = 1 \end{eqnarray*}$ einen Wert $\begin{eqnarray*} E(X) \approx 0.10856 \end{eqnarray*}$ erwarten lässt, d.h. gewiss nicht die von dir erwarteten 0.5. unglücklich

EDIT: An sich hätte ich ja noch irgendeine Reaktion erwartet - nicht die des Vogels Strauß, der seinen Kopf in den Sand steckt...

07.03.2012, 15:37

vogelstrauss

Auf diesen Beitrag antworten »

Hallo HAL,

Entschuldigung fuer das Kopf-in-den-Sand-Stecken; nach der Antwort hab ich mich nicht getraut, nochmal nachzufragen, obwohl ich es immer noch nicht verstanden habe (wegen dem Titel "Ohne weiteren Kommentar", ich dachte das Problem ist wohl zu trivial...). Zudem bekomme ich leider keine email-Benachrichtigung, wenn ein Kommentar nochmal editiert wird, deshalb habe ich nicht gewusst, dass noch was kam...

"Logarithmische Gleichverteilung" war wohl nicht das, was ich meinte.
Ich beschreibe mal mein ganzes Problem:

Ich habe einen Vektor aus 5 parametern, die ich mit Metropolis gleichzeitig samplen will. Diese 5 Parameter will ich mit meinen Messdaten vergleichen und mit least-squares die Kosten fuer diese Parameter berechnen. Ausserdem weiss ich aus der Literatur, dass Parameter #5 normalerweise Werte zwischen 0 und 1 hat und sehr selten Werte ein wenig ueber 1.
Also moechte ich, dass fuer Parameter #5 Werte zwischen 0 und 1 mit hoeherer W-keit gesampled werden. dazu verbinde ich eine Gleichverteilung mit einer Lognormalverteilung (mit mean=1).
Da es keinen Sinn macht, dass ein Parameter <0 ist, moechte ich gerne im Logspace samplen, die least-squares-Berechnung soll aber mit nicht-logarithmierten Werten arbeiten.
Also habe ich eine 'normale' und eine log-Kostenfunktion:

code:

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:

cost = function(params){
  x <- params[5]
  ret <- -Inf
  if(log(x)<=1 && log(x)>=0)
    ret <- (((dunif(log(x), min=0, max=1, log=TRUE))))
  else
    ret <- ((dlnorm(x, meanlog=1, sdlog=0.3, log=TRUE)))
  ##leastsq_cost <- my_leastsq_cost_function(params) ##Hier werden die Parameter mit den Daten verglichen                                                                           
  ##return(ret - leastsq_cost)                                                                                                                                                      
  return(ret)
}

cost_log <- function(params){
  return (cost(exp(params)))
}

Die Berechnung der Least-squares kosten ist hier weggelassen, aber das ist hier glaube ich auch nicht das Problem.
Wenn ich jetzt sample und die Verteilung des 5. Parameters plotte:

code:
1: 2: 3: 4: 5:	require(mcmc) pstart <- c(1,1,1,1,0.5) out <- metrop(cost_log, log(pstart), 10000) hist(out$batch[,5])

bekomme ich auch genau die Verteilung, die ich wollte.
Aber hier liegt mein grosses Problem:
Ich will ja (und hab ja) im Logspace gesampled!!
Und deshalb koennen in meinem batch die Werte von den anderen 4 Parametern
negativ werden, was ich ja verhindern wollte!
Was ich genau will, ist dass exp(out$batch[,5]) mir die gewuenschte Verteilung (uniform in [0,1] und lognormal fuer >1) gibt!

Ich verstehe einfach nicht, wie das funktionieren kann und ich verstehe auch nicht, was in der Hilfe mit "log unnormalized probability density" gemeint ist...
Zum Beispiel hier

code:
1:	dunif(log(x), min=0, max=1, log=TRUE))

gebe ich doch den Logarithmus der Dichtefunktion zurueck, das ist doch dann nicht "log unnormalized"??

Sorry fuer den langen Beitrag, aber jetzt hab ich das komplette Problem beschrieben.
Falls Sie noch interessiert sind, mir zu helfen, waere ich sehr dankbar. Ich denke, ich verstehe die Essenz der logarithmischen Verteilungen noch nicht.

viele Gruesse

vogelstrauss

07.03.2012, 16:42

HAL 9000

Auf diesen Beitrag antworten »

Ich bin nicht völlig unbewandert mit MCMC, insbesondere Metropolis-Hastings. Und nachdem ich zwischenzeitlich die Dokumentation von deinem R-Paket "mcmc" überflogen habe - leider konnte ich es auf meiner Windows-Version von R nicht zum Laufen bringen, weil ich keine vorkomplierte Windows-DLL-gefunden habe und keine Lust hatte, mich selber durch die Erzeugung einer solchen zu kämpfen - bin ich mir vollkommen sicher, dass mit "log unnormalized density" einfach der Logarithmus der Dichte deiner zu simulierenden Verteilung gemeint ist. Allerdings musst du diese Dichte $\begin{eqnarray*} f(x) \end{eqnarray*}$ nur bis auf einen Vorfaktor genau kennen, d.h. du musst nur ein $\begin{eqnarray*} g \end{eqnarray*}$ kennen, für das es einen Proportionalitätsfaktor $\begin{eqnarray*} c>0 \end{eqnarray*}$ gibt, so dass $\begin{eqnarray*} f(x)=cg(x) \end{eqnarray*}$ für alle $\begin{eqnarray*} x \end{eqnarray*}$ gilt. Man symbolisiert das gern mit $\begin{eqnarray*} f(x) \propto g(x) \end{eqnarray*}$ .

In dem Sinne kann man für $\begin{eqnarray*} a<b \end{eqnarray*}$ eine Gleichverteilung auf dem Intervall $\begin{eqnarray*} [a,b] \end{eqnarray*}$ schlicht und einfach durch die nicht-normalisierte Dichte

$\begin{eqnarray*} g(x) = \begin{cases} 1 & \;\mbox{für}\; a\leq x\leq b \\ 0 & \;\mbox{sonst} \end{cases} \end{eqnarray*}$

beschreiben, in dem Sinne wäre dann der Proportionalitätsfaktor $\begin{eqnarray*} c=\frac{1}{b-a} \end{eqnarray*}$ .

Deine R-Prozedur muss nun den Logarithmus einer solchen nicht-normalisierten Dichte liefert, hier bei der Gleichverteilung wäre das demnach

$\begin{eqnarray*} \ln(g(x)) = \begin{cases} 0 & \;\mbox{für}\; a\leq x\leq b \\ -\infty & \;\mbox{sonst} \end{cases} \end{eqnarray*}$ .

Während $\begin{eqnarray*} g(x) \end{eqnarray*}$ also nur bis auf einen Faktor eindeutig bestimmt ist, bewirkt der Logarithmus, dass $\begin{eqnarray*} \ln(g(x)) \end{eqnarray*}$ nur bis auf eine additive Konstante eindeutig bestimmt ist. D.h., der Metropolis-Algorithmus arbeitet für jede beliebige Wahl von $\begin{eqnarray*} K \end{eqnarray*}$ in

$\begin{eqnarray*} \ln(g(x)) = \begin{cases} K & \;\mbox{für}\; a\leq x\leq b \\ -\infty & \;\mbox{sonst} \end{cases} \end{eqnarray*}$

gleich - alle diese "log unnormalized densities" beschreiben ein- und dieselbe Verteilung, eben jene stetige Gleichverteilung auf $\begin{eqnarray*} [a,b] \end{eqnarray*}$ . Augenzwinkern

07.03.2012, 18:02

vogelstrauss

Auf diesen Beitrag antworten »

Vielen Dank fuer die Hilfe,
Das mit den log-unnormalized densities habe ich jetzt glaub ich verstanden, auch das mit der Konstanten $\begin{eqnarray*} K \end{eqnarray*}$ ist klar, da Die Dichte in $\begin{eqnarray*} [a,b] \end{eqnarray*}$ ueberall gleich ist.
Da ich fuer Werte >1 die (logarithmierte) Dichte der Normalverteilung zurueckgeben will, sieht meine logarithmierte Dichte nun so aus:

$\begin{eqnarray*} ln(g(x))=\begin{cases} 0, & \text{fuer } a\leq x\leq b\\ ln(N(1, 0.3)), & \text{fuer } x > b \\ -\infty & \text{fuer } x<0 \end{cases} \end{eqnarray*}$

Wenn ich das mit mcmc sample funktioniert das auch sehr gut, das Histogramm der Verteilung sieht so aus, wie ich es erwarte.

Mein Problem ist, dass beim Samplen recht viele negative Werte vorgeschlagen werden, bei der die Dichte immer $\begin{eqnarray*} - \infty \end{eqnarray*}$ ist, was natuerlich dann mein acceptance-ratio verringert. Negative Werte machen in meinen Berechnungen keinen Sinn, da ich Messdaten simuliere.
Deshalb wuerde ich gerne im Logspace samplen, so dass alle vorgeschlagenen Parameterwerte positiv sind.
Ich wuerde also gerne anstatt

code:

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:

cost <- function(params){
  x <- params[5]
  if(min(params)<0)
    return (-Inf)
  if (x>=0 && x<=1)
    return (0)
  else
    return(dnorm(x, mean=1, sd=0.3, log=TRUE))
} 
out = metrop(cost, c(1,1,1,1,0.5), 10000)

die Logarithmen meiner Startwerte in die metrop()-Funktion geben und die gesamplete Verteilung danach wieder exponieren:

code:
1: 2: 3: 4: 5:	cost_log <- function(params){ return(cost(exp(params))) } out = metrop(cost_log, log(c(1,1,1,1,0.5)), 10000) hist(exp(out$batch[,5]))

und dann die Verteilung erhalten, die ich auch beim normalen samplen, ohne Logspace bekommen wuerde. Macht das Sinn?
In dem code Beispiel hier funktioniert das natuerlich nicht und ich verstehe nicht, warum.
Wie kann ich im logspace so samplen, dass der Exponent meiner posterior-Verteilung (exp(batch)) die oben genannte Dichte hat?

Ich hoffe, mein Problem ist verstaendlich,

Gruesse

vogelstrauss

07.03.2012, 18:18

HAL 9000

Auf diesen Beitrag antworten »

Nein, so richtig verstehe ich nicht, was du da genau betrachtest. Das mag u.a. auch an der befremdlichen Symbolik liegen - man muss sich alles erstmal aus dem R-Quelltext mühsam rücküberlegen, weil die mathematischen Formeln nur sehr schwer zu gebrauchen sind:

Was soll z.B. $\begin{eqnarray*} N(1,0.3) \end{eqnarray*}$ sein? Soll das die Dichte der $\begin{eqnarray*} N(1,0.3) \end{eqnarray*}$ -Verteilung an der Stelle $\begin{eqnarray*} x \end{eqnarray*}$ sein? Da hätte ich zumindest $\begin{eqnarray*} x \end{eqnarray*}$ in der Parameterliste jender Funktion erwartet, also etwa $\begin{eqnarray*} N(x,1,0.3) \end{eqnarray*}$ , wenn schon nicht in der allgemein gebräuchlichen Variante $\begin{eqnarray*} \frac{1}{0.3}\cdot \varphi\left( \frac{x-1}{0.3} \right) \end{eqnarray*}$ mit Dichte $\begin{eqnarray*} \varphi \end{eqnarray*}$ der Standardnormalverteilung.

Sätze wie

Zitat:

Original von vogelstrauss
Wie kann ich im logspace so samplen, dass der Exponent meiner posterior-Verteilung (exp(batch)) die oben genannte Dichte hat?

finde ich erst recht unverständlich. Du darfst nicht davon ausgehen, dass jeder die Sonderbezeichnungen ("batch" ???) deines R-Pakets kennt. unglücklich

Zitat:

Original von vogelstrauss
Mein Problem ist, dass beim Samplen recht viele negative Werte vorgeschlagen werden, bei der die Dichte immer $\begin{eqnarray*} - \infty \end{eqnarray*}$ ist, was natuerlich dann mein acceptance-ratio verringert.

Das liegt ja nicht an der zu sampelnden Dichte, sondern an der Wahl der Vorschlagsverteilung ("proposal distribution"). Im einfachen Metropolis-Algorithmus ist das ja schlicht Random-Walk, aber auch da kann man wenigstens die Schrittweite wählen. Ich vermag nicht zu erkennen, wo das bei dir im R-Quelltext geschieht, dazu kenne ich das Paket "mcmc" zu wenig. unglücklich

08.03.2012, 19:53

vogelstrauss

Auf diesen Beitrag antworten »

Ich versuche gerade, das alles in korrekter Notation aufzuschreiben, habe leider gestern und heute nicht so viel Zeit gehabt, ich bleib aber dran!!

15.03.2012, 16:55

vogelstrauss

Auf diesen Beitrag antworten »

Also, ich probier nochmal, mein Problem zu erlaeutern:

Ich moechte meine 5 Parameter gleichzeitig mit Metropolis-Hasting samplen, wobei der fuenfte Parameter die Dichte
$\begin{eqnarray*} g(x)=\begin{cases} 1, & \text{wenn }a \leq x \leq b \\ \frac{1}{\sigma \sqrt{2 \pi}} exp(- \frac{(x- \mu)^2}{2 \sigma^2}) \text{mit }\mu=1,\sigma=0.3, & \text{wenn }x> a \\ 0 , & \text{wenn } x < 0 \\ \end{cases} \end{eqnarray*}$
haben soll.

Weil die Energie-Funktion in der Implementation im mcmc-Paket den Logarithmus der nichtnormalisierten Dichte verlangt, uebergebe ich der Metropolis-Funktion die folgende Dichte:

$\begin{eqnarray*} ln(g(x))=\begin{cases} 0, & \text{wenn }a \leq x \leq b \\ ln(\frac{1}{\sigma \sqrt{2 \pi}} exp(- \frac{(x- \mu)^2}{2 \sigma^2})) \text{mit }\mu=1,\sigma=0.3, & \text{wenn }x> a \\ - \infty , & \text{wenn } x < 0 \\ \end{cases} \end{eqnarray*}$
,also den Logarithmus von $\begin{eqnarray*} g(x) \end{eqnarray*}$ .

Wenn ich diese Verteilung num mit dem R-Paket sample, sieht die Posterior-Verteilung des Parameters so aus wie ich es erwarte: Ein Histogramm der Posterior-Verteilung zeigt, dass der Parameter gleichverteilt zwischen 0 und 1 ist ( $\begin{eqnarray*} a=0, b=1 \end{eqnarray*}$ ), fuer Werte > 1 wird die Dichte geringer.

Zitat:

Du darfst nicht davon ausgehen, dass jeder die Sonderbezeichnungen ("batch" ???) deines R-Pakets kennt

Was mir die Metropolis-Funktion als "batch" zurueckgibt, ist genau die o.g. Posterior-Verteilung der Parameter.

So weit, so gut.

Jetzt weiss ich zusaezlich, dass keiner meiner Parameter <=0 sein kann. Fuer den fuenften Parameter mit der oben beschriebenen Posterior-Verteilung ist das auch kein Problem, es steht ja in $\begin{eqnarray*} g(x) \end{eqnarray*}$ , dass die Wahrscheinlichkeit bei Werten <0 gleich 0 ist.

Zitat:

Man kann in dem R-Paket die Schrittweite waehlen, die Sache ist nur, dass wenn ich die Schrittweite zu klein mache, werden vielleicht nicht so viele negative Werte vorgeschlagen, aber der Algorithmus brauch laenger, zu konvergieren.
Daher wuerde ich gerne im logspace samplen. Ich wuerde als Startwerte die Logarithmen meiner Startwerte der Metropolis-Funktion uebergeben und die Werte einer Posterior-Verteilung $\begin{eqnarray*} g'(x) \end{eqnarray*}$ erhalten. Diese Posterior-Verteilung sollte so aussehen, dass $\begin{eqnarray*} ln(g'(x))=g(x) \end{eqnarray*}$ .

Mein Problem ist nur, wie muss ich dann $\begin{eqnarray*} g'(x) \end{eqnarray*}$ definieren??
Wenn ich die Wahrscheinlichkeit nicht von x, sondern von exp(x) berechne, also wenn der Logarithmus der nichtnormalisierten Dichtefunktion so aussieht:

$\begin{eqnarray*} ln(g'(x))=\begin{cases} 0, & \text{wenn }a \leq exp(x) \leq b \\ ln(\frac{1}{\sigma \sqrt{2 \pi}} exp(- \frac{(exp(x)- \mu)^2}{2 \sigma^2})) \text{mit }\mu=1,\sigma=0.3, & \text{wenn }exp(x)> a \\ \end{cases} \end{eqnarray*}$

bekomme ich nicht die gewuenschte Dichte $\begin{eqnarray*} g'(x) \end{eqnarray*}$ so dass $\begin{eqnarray*} ln(g'(x)) = g(x) \end{eqnarray*}$ . Ich habe das in R ausprobiert und bekomme eine andere Verteilung.

Ich hoffe mein Problem ist ein bisschen genauer beschrieben...

15.03.2012, 17:11

HAL 9000

Auf diesen Beitrag antworten »

Zitat:

Original von vogelstrauss
Man kann in dem R-Paket die Schrittweite waehlen, die Sache ist nur, dass wenn ich die Schrittweite zu klein mache, werden vielleicht nicht so viele negative Werte vorgeschlagen, aber der Algorithmus brauch laenger, zu konvergieren.

Das ist der Preis:

kleine Schrittweite = langsame Konvergenz, hohe Akzeptanzraten

große Schrittweite = schnelle Konvergenz, falls man hohe Akzeptanzraten hätte - hat man da aber leider nicht

Irgendwo gibt es ein Optimum, zumindest etwas, was dem nahekommt.

Auf deine Ausführungen will ich nicht im einzelnen eingehen, zumal ich immer noch nicht den Begriff "im logspace sampeln" verstehe. Was ich nicht begreife, ist folgendes:

Du willst doch eine feste Verteilung simulieren, da ist das $\begin{eqnarray*} g \end{eqnarray*}$ doch aber (bis auf einen Vorfaktor) festgelegt, da kannst du doch nicht dran herumdrehen. Das einzige, woran du in Hinblick auf Konvergenzrate da drehen darfst, ist doch bei Metropolis-Hastings die Vorschlagsverteilung (bei Metropolis reduziert sich diese Stellschraube leider nur auf die Wahl der Schrittweite). verwirrt

Insofern ist mir dein Herumdoktorn am $\begin{eqnarray*} g \end{eqnarray*}$ ein Buch mit sieben Siegeln. unglücklich

Neue Frage »

Antworten »

Logarithmische Gleichverteilung

Verwandte Themen