Marginal likelihood - implementation of the Chib/Jeliazkov method

[florianhartig]

trafficstars

Question from a user:

ich habe eine Frage zu Ihrem Rcode für die Chib/Jeliazkov Methode in der Funktion marginalLikelihood. Dort steht, dass theta.star für dmvnorm und rmvnorm als Mittelwert verwendet wird. Ist dies richtig? Das bezieht sich doch auf q(theta,theta.star|y) bzw q(theta.star,theta|y) im Chib/Jeliazkov paper (Gleichung 8 und 9) und ich hatte die Notation so verstanden, dass dies nur im zweiten Fall so ist? Bin gerade auf Fehlersuche in meiner eigenen Praxisanwendung, und hatte dies erst anders implementiert.

[florianhartig]

Tankred / Max - if you have time for this, can you have a look at this? In doubt, leave it until I have time for this.

[MaximilianPi]

Maybe he's right, it's about Equatation 9 and they say that theta.g is sampled from the posterior distribution, which is done:

lik <- chain[,sampler$setup$numPars + 2]
g <- sample.int(nrow(x), numSamples, replace=TRUE)
alpha.g <- sapply(lik.g, function(l) min(1, exp(lik.star - l)))

# q(theta.g, theta.star | y) is calculated by:
q.g <- mvtnorm::dmvnorm(x[g,], mean = theta.star, sigma = V, log = FALSE)

theta.j should be drawn from q(theta.star, theta | y)

theta.j <- mvtnorm::rmvnorm(numSamples, mean = theta.star, sigma = V)
lik.j <- apply(theta.j, 1, sampler$setup$likelihood$density)

alpha.j <- sapply(lik.j, function(l) min(1, exp(l - lik.star)))

[florianhartig]

Does that make a difference for the calculated values?

[florianhartig]

OK, yes, it looks as if there is a difference of our implementation to the paper, but but I wrote a few further tests, and our calculated values seem totally fine, see https://github.com/florianhartig/BayesianTools/commit/e50350914736912c0a57c3c64abbe70e6886db01

so it seems it is not as urgent as I thought to fix this.

Max, can you implement the Chib/Jeliazkov version and make a comparison of the error for both options? I'll try to dig deeper into the motivation for the two options.