2
$\begingroup$

I am drawing a sample Y of size n from a p-dimensional Normal ($\mu, \Sigma$). Typically, p is 5. I have $\bar{Y}$ and $V = YY'$, the sum of squares. Now I want to draw samples from this $\bar{Y}$, given V. The density of $\bar{Y}$ given V is proportional to $det(v - n \bar{y} \bar{y}')^{(n - p - 2)/2}$ with support$ \bar{y} : v - n \bar{y} \bar{y}'$ is positive definite.

This is not a Bayesian framework, but possibly (at least I wasn't able to figure out) MCMC is the only way to draw samples from this not-so-well-known density. I used Metropolis-Hasting's algorithm with a normal proposal density, mean being this $\bar{y}$ and proposal variance being the observed variance-covariance matrix, divided by n. This is the unconditional distribution of $\bar{y}$. This gives me good-looking traceplots but an acceptance rate of about 2%! I tried lowering the proposal variance (say, multiplying by 0.2), but that makes the acceptance rate even low! Isn't acceptance rate supposed to increase if you lower the proposal variance?

Is there any simpler way out other than using adaptive MCMC? I've also been thinking of using other proposal distributions but couldn't come up with anything reasonable.

$\endgroup$
  • $\begingroup$ Are you sure of this conditional distribution? I would have thought $\bar Y$ and $V$ to be independent. $\endgroup$ – Xi'an Jul 30 '18 at 9:19
  • $\begingroup$ @Xi'an As functions of the same sample of random variables, I don't think they can be independent, but I'm also curious how the OP derived this density. $\endgroup$ – deasmhumnha Jul 30 '18 at 15:54
0
$\begingroup$

I can't diagnosis why the rejection rate increased without knowing more specifics about your implementation, but I can help with a better sampling method.

Using the matrix determinant lemma, we can rewrite the density as $$\mathop{\mathrm{det}} (V)^{\frac{n-p}{2}-1}(1-n\bar{y}^TV^{-1}\bar{y})^{\frac{n-p}{2}-1}$$ and we see that our desired support refers to those points such that $0\leq \bar{y}^TV^{-1}\bar{y} \leq \frac{1}{n}$ since the determinant of a PD matrix is positive and a quadratric form with a PD matrix is also positive. We see then that all vectors with equal values for $\bar{y}^TV^{-1}\bar{y}$ are equally probable, so we can first sample a value $c$ according to $(1-nc)^{\frac{n-p}{2}-1}$ and then choose a vector such that $\bar{y}^TV^{-1}\bar{y}=c$. Note that $c$ can be simulated as $x/n$ where $x\sim \mathrm{Beta}\left(1, \frac{n-p}{2}\right)$. Next we simulate a random vector $y\sim \mathcal{N}(0,V)$ and set $$\bar{y}=\left(\frac{c}{y^TV^{-1}y}\right)^{\frac{1}{2}}y$$

For the sake of efficiency, you will probably want to use Cholesky $V=U^TU$ to calculate the determinant, the matrix-vector products with $V^{-1}$ and to generate the random normal vectors as $y=U^Tw$ where $w\sim\mathcal{N}(0,I)$. Using Cholesky also leads to the simplification that $y^TV^{-1}y=w^Tw$ and thus $$\bar{y}=\frac{\sqrt{c}}{\|w\|_2}U^Tw$$

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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