5
$\begingroup$

I have a multivariate Gaussian defined as follows: $$ p(x) = \omega(x)\gamma(x) $$ where $\omega$ and $\gamma$ both are multivariate Gaussians and from which I can sample very efficiently given due to special structure in their covariances: one is diagonal, of the other I know the diagonalization, which is very sparse.

I wonder if there is a way to sample from $p$, given samples from $\omega$ and $\gamma$. The standard way of doing a Cholesky on the joint covariance is not efficient enough.

Improve this question
$\endgroup$
1
  • 1
    $\begingroup$ Presumably there's an implicit normalization constant in the product and $p$, $\omega$, and $\gamma$ refer to PDFs. In that case $p$ also is a multivariate Gaussian. Whether its variance has any special structure depends on the means and eigensystems of $\omega$ and $\gamma$. As far as I can tell, then, absent any specific information about those relationships, the only possible (valid) answers to this question will simply explain how to sample from an arbitrary multivariate Gaussian. $\endgroup$
    – whuber ♦
    Sep 18, 2014 at 17:14

2 Answers 2

Reset to default
4
$\begingroup$

Rejection sampling might work: sample from one distribution and accept based on the other density

  1. Sample $x$ from $\omega(x)$.
  2. Return $x$ with probability $\gamma(x) / \gamma(\mu_\gamma)$, otherwise goto 1.

The roles of $\gamma$ and $\omega$ may be switched due to symmetry, so one could try both to see which is more efficient.

However, if neither component distribution is 'close' to the joint distribution, this will not be efficient as most samples get rejected.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Since $p$ is also multivariate Gaussian, it seems unlikely that this rejection sampling method would be superior to any generic method of sampling from multivariate Gaussians. Under what conditions on $\gamma$ and $\omega$ do you suppose rejection sampling would have better performance? $\endgroup$
    – whuber ♦
    Sep 18, 2014 at 17:16
1
$\begingroup$

First, note that the product $$ p(x) = \omega(x)\gamma(x), $$ is not a density function, unless you normalise it.

What you can do is to take $\omega(x)$ as an instrumental "likelihood" and $\gamma(x)$ as an instrumental "prior" (not in a strict sense, it is only a trick to construct a tractable sampler). This leads you to a conjugate model (since the normal prior is the conjugate prior of the mean of a normal sampling model):

http://en.wikipedia.org/wiki/Conjugate_prior#Continuous_distributions

This product should be normal, with a certain mean and covariance structure, which is tractable due to conjugacy.

Improve this answer
$\endgroup$
5
  • $\begingroup$ The (normalized) product is indeed normal, but the last sentence of the question mentions that computing the mean and covariance of the 'posterior' and sampling from that in the standard way is not efficient enough. $\endgroup$
    – Juho Kokkala
    Sep 18, 2014 at 10:09
  • $\begingroup$ @JuhoKokkala I wonder why is it that difficult since the variance of the conjugate posterior is kind of simple. $\endgroup$
    – Coco
    Sep 18, 2014 at 10:20
  • $\begingroup$ That's what I meant with the last paragraph. $\endgroup$
    – bayerj
    Sep 18, 2014 at 11:07
  • $\begingroup$ @bayerj What's the dimension of $x$? $\endgroup$
    – Coco
    Sep 18, 2014 at 11:15
  • 1
    $\begingroup$ The dimension can be anything, but practical are a few hundred. However, There will be thousands of these problems in an inner loop of an optimisation, which should run as fast as possible. $\endgroup$
    – bayerj
    Sep 18, 2014 at 11:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.