6
$\begingroup$

Motivating question

I have a high-dimensional state space $\Omega \subseteq \mathbb R^n$ with an admissible subset $S\subseteq \Omega$, which is connected. I would like to draw a uniform random sample from $S$.

In my application, it is easy to verify whether a state $\vec x$ is in $S$, but difficult to find a state in $S$ ad hoc. However, it is known that $\vec 0 \in S$.

Solution idea

I think the problem should be relatively easy to solve with a Metropolis-Hastings algorithm.

  1. Start at $x_0 = \vec 0$. Set $i:=0$.
  2. Set $i:=i+1$.
  3. Randomly generate a close-by state $\vec x_i:= \vec X( \vec x_{i-1})$ based on $\vec x_{i-1}$.
  4. Accept the step if $\vec x_i \in S$, else set $\vec x_i = \vec x_{i-1}$.
  5. Add $\vec x_i$ to the sample.
  6. Repeat the procedure after step 2.

We may need to through away a large number of steps from the burn-in period.

Question

I am wondering which properties the random variable $\vec X(\vec x)$ needs to have to lead to a uniform sample over the state space. Does the step distribution need to satisfy detailed balance? Why?

Example

Draw a uniform sample from the unit disk in 2D and use the bivariate standard normal distribution centred at $\vec x$ for the steps, i.e., $\vec X(\vec x) \sim \vec N(\vec x, \underline{1})$. I think this should work, but I also think that biasing the steps towards the centre would lead to a different sample.


The question is so basic that I would expect to find lecture notes that help me, but so far I was unsuccessful in finding something where the Metropolis algorithm was interpreted in a Bayesian framework with prior and posterior distribution. A corresponding reference might be perfectly fine as an answer.

$\endgroup$
8
  • $\begingroup$ (1) Could you please explain how a subset of a product whose components include the disconnected natural numbers could be "connected" unless each of those components is constant (and therefore can be ignored)? I suspect you have some specialized concept of "connected" in mind. (2) What is a "unit cycle"? Unit circle, perhaps? If so, do you really mean the circle or do you mean the disk? $\endgroup$
    – whuber
    Apr 14, 2023 at 15:45
  • 1
    $\begingroup$ Excellent comments @whuber. The statement of connectedness was only thought to be useful to argue why searching for admissible states in a neighbourhood of an admissible state makes sense. For that, a weaker concept of connectedness may be applicable. Because I am missing the right term at the moment, I have adjusted the question, as I think it applies in the same way if I use a product set of real numbers. And yes, I mean the unit disk (question adjusted). $\endgroup$
    – Samufi
    Apr 14, 2023 at 16:27
  • $\begingroup$ What you are describing, I believe, is exactly the Metropolis algorithm with a uniform prior. The likelihood of your data if not in S is 0, and 1 if in S, therefore the probability of acceptance is 1 anywhere you draw a point in S, and 0 otherwise. I am not sure why you expect this problem to be described in lecture notes. It's ok to use a Normal for the proposal, all you need is a symmetric distribution $\endgroup$
    – andrea m.
    Apr 14, 2023 at 20:01
  • $\begingroup$ @andream. Yes, I constructed the algorithm that way by intention. What I am missing is under which conditions the prior is uniform. You say that the proposal distribution needs to be symmetric. Which concept of symmetry (especially when working on more complicated spaces)? If it is not symmetric, how does it relate to the prior (how could I create a proposal for a specific prior I desire)? (The second question is beyond my original question but would be great for my understanding.) $\endgroup$
    – Samufi
    Apr 15, 2023 at 9:48
  • $\begingroup$ Symmetric means that the probability of going to point B from point A is the same as going from B to A. $\endgroup$
    – andrea m.
    Apr 16, 2023 at 14:09

2 Answers 2

8
+50
$\begingroup$

Since you want to sample uniformly, $p(\xi) = c I(\xi \in S)$ so the ratio $\frac{p(\xi^*)}{p(\xi)} = I(\xi^* \in S)$ for all proposals $\xi^*$. In MH you accept with $\frac{p^*}{p}\frac{q}{q^*}$ but, with a symmetrical proposal, you only care about $p$ since the factor $\tfrac{q(\xi\leftarrow \xi^*)}{q(\xi^*\leftarrow \xi)} = 1$ cancels out (This is Metropolis without Hastings inspiring work). Therefore all proposals in $S$ will be accepted. The stationary distribution of $\xi$ will follow $p$.

from numpy.random import randn as q
xi = {}
G  = 10**4
x  = y = 0
for g in range(G):
    xs, ys = x+q(), y+q()                         # xi_star ~ q
    x, y   = ((x,y),(xs,ys))[xs**2 + ys**2 <= 1]  # xi = xi_star if xi_star\in S
    xi[g]  = x, y
import numpy as np
from matplotlib.pyplot import plot,subplot
xi = np.vstack(xi.values())
subplot(211)
plot(*xi.T, '.')
subplot(212)
plot(xi)
Notation:
$p$ target density
$q$ proposal
$\endgroup$
2
  • $\begingroup$ I think I get the point but am unsure about some parts. Could you add some details on the symbols (what are $p$ and $q$) and where the $q$ parts cancels out (like, in which formula, and why is the quotient important)? How is that related to all proposals in S being accepted? If you could add the missing details, this might also be of great benefito to future readers of your answer. $\endgroup$
    – Samufi
    Apr 21, 2023 at 20:05
  • $\begingroup$ Just one note on a good answer: even though the acceptance probability is 1 for proposals in the set $S$, in high-dimensional spaces you can still get stuck in corners where most of the potential steps take you outside $S$. $\endgroup$ May 25, 2023 at 2:49
10
$\begingroup$

What you're looking for is almost exactly what "Metropolis-Hastings for Ratio-of-Uniforms" tries to do. There are unfortunately no known papers on this topic, but what is available are some slides by Luke Tierney: http://homepage.stat.uiowa.edu/~luke/talks/bormio05.pdf

Ratio-of-Uniforms is a way to generate random samples from a given distribution $F$ with density $f$ (need not be defined on any bounded set). However, in order to implement Ratio-of-Uniforms, one needs to obtain uniform draws from a bounded set $\mathscr{A}$ where this is: $$ \mathscr{A} = \left\{ (u,v): 0 \leq u \leq \sqrt{f(v/u)} \right\}\,. $$

If $(U,V) \sim \text{Unif}(\mathscr{A})$, the $V/U \sim F$. This is typically done by enclosing the set $\mathscr{A}$ within a rectangle and then implementing iid Accept-Reject sampling. However, this can be fairly inefficient if $\mathscr{A}$ is weird and high-dimensional.

Tierney presents a few different ideas on how Markov chains can be constructed to sample uniformly on this set $\mathscr{A}$ instead of rejection sampling.

There is also this Masters Thesis project that explore this for a specific problem: https://etda.libraries.psu.edu/files/final_submissions/1226

$\endgroup$

Your Answer

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

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