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$
2
$\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 '20 at 3:18

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.