1
$\begingroup$

I understand the Metropolis algorithm. Where I get confused is the MH algorithm where asymmetric proposal distributions may be used.

Formula

I understand that P(x) and P(x') represent the likelihood/probability density of x and x' according to the target distribution. Likewise, I understand that g(x|x')/g(x'|x) is a term used to correct an asymmetrical proposal distribution. I am not confused by its purpose; I don't understand its execution.

As a toy problem, I've developed an exponential distribution sampler. There are two variations, one that uses a symmetric proposal distribution, the uniform dist. And one that does not: Namely, Beta(a=3,b=2) - 0.5. I've chosen this distribution because (A) it is asymmetric and mostly positive (but occasionally negative, due to the -0.5 term.)

I have no idea how to find g(x|x')/g(x'|x).

Code:

def target(x,lam):
    return int(x>0) * lam * np.exp(-x * lam)

def exponential_MH(hops,lam=3):
    states = []
    burn_in = int(hops*0.2)
    current = lam
    
    for i in range(hops):
        states.append(current)

#         movement = current + random.uniform(-1,1) # does not require asymmetric correction
        movement = current + np.random.beta(a=3,b=2)-0.5 # requires asymmetric correction

        acceptance = target(x=movement,lam=lam)/target(x=current,lam=lam)
        event = random.uniform(0,1)
        if acceptance > event:
            current = movement
            
    return states[burn_in:]        
        

lam = 1
exp_samples = exponential_MH(hops=10_000,lam=lam)
lines = np.linspace(0,5,10_000)
exp_curve = [lam*np.exp(-l*lam) for l in lines]

plt.hist(exp_samples,normed=1,bins=20) 
plt.plot(lines,exp_curve)

Asymmetric proposal distribution plot:

Asymmetric proposal distribution

Symmetric proposal distribution plot (with commented lines flipped):

Symmetric proposal distribution plot

To answer this question, please edit code to reflect a valid g(x|x')/g(x'|x) for the proposal distribution Beta(a=3,b=2) - 0.5 which perturbations are drawn from.

$\endgroup$

1 Answer 1

3
$\begingroup$

Take a look at the updated code and plot below. Notice that g(x|x')/g(x'|x) is essentially a measure of how likely these perturbations are to be seen under proposal distribution, which has been defined as Beta(a=3,b=2) -0.5.

First, find the difference between the current and proposed events. Second, adjust for the -0.5; we'll call these unbiased perturbations (where -0.5 is a bias.) third find the likelihood of each perturbation (curr->prop & prop->curr). Lastly, return the ratio as correction.

We'll use this correction and multiply it with other terms in the acceptance variable definition. That's pretty much it!

def target(x,lam):
    return int(x>0) * lam * np.exp(-x * lam)

def correct(prop,curr,a=3,b=2):
    x0 = curr - prop + 0.5
    x1 = prop - curr + 0.5
    b0 = beta.pdf(x=x0, a=a, b=b)
    b1 = beta.pdf(x=x1, a=a, b=b)
    return b0/b1 

def exponential_MH(hops,lam=3):
    states = []
    burn_in = int(hops*0.2)
    current = 1
    
    for i in range(hops):
        states.append(current)
        movement = current + np.random.beta(a=3,b=2)-0.5 # requires assymetric correction        
        correction = correct(curr=current,prop=movement)
        acceptance = target(x=movement,lam=lam)/target(x=current,lam=lam)*correction
        event = random.uniform(0,1)
        if acceptance > event:
            current = movement
            
    return states[burn_in:]        

Updated sampler

$\endgroup$
1
  • $\begingroup$ Out of curiosity, what would this solution look like for a multivariable/joint distribution? I know that Gibbs can do well in these situations; but it's my understanding that MH can handle more than a single variable. $\endgroup$
    – mjake
    Aug 18, 2020 at 3:18

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.