Download this notebook

MCMC Metropolis sampling

Table of Contents

Introduction

The Metropolis algorihtm is a simple and powerfull sampling method in the framework of Markov Chain Monte Carlo procedures. It works fairly similar to rejection sampling with the execption of using the previous result to adjust the next step and therefore build up a Markov Chain while sampling. In this notebook you will implement the Metropolis algorihtm and draw samples for two given probability distribution.

Requiremenst

Knowledge

To complete this exercise notebook, you should possess knowledge about the following topics.

  • MCMC Sampling (Metropolis algorithm)

The following material can help you to acquire this knowledge:

  • Bishop [BIS06], Chapter 11.2

Python Modules

import numpy as np
from matplotlib import pyplot as plt
import scipy.stats as sc
np.random.seed(42)
%matplotlib inline

Exercise

Consider the following two distribution functions,$ p_1 = 0.4\mathcal N(2,0.3) + 0.6\mathcal N(4,0.4) $ and$ p_2 = \frac{1}{2 \pi} (cos(x)+1) $ with$ x \in (0,2 \pi) $.

def desired_prop_distrubution1(x):
    
    gaus1 = sc.norm(loc=2,scale=0.3)
    gaus2 = sc.norm(loc=4,scale=0.4)
    double_gaussian = 0.4*gaus1.pdf(x) + 0.6*gaus2.pdf(x)
    
    return double_gaussian

def desired_prop_distrubution2(x):
    
    
    if np.size(x)>2:
        y = np.zeros_like(x)
        for i in np.arange(np.size(x)):
            if (x[i]<=2*np.pi) and (x[i]>=0):
                   y[i] = (np.cos(x[i])+1)/(2*np.pi)
            else:
                   y[i] = 0
    else:
        if (x<=2*np.pi) and (x>=0):
            y = (np.cos(x)+1)/(2*np.pi)
        else:
            y = 0
        
        
    return y

xx=np.linspace(-1,7,1000)
plt.figure(figsize=(17,5))
plt.subplot(121)
plt.plot(xx,desired_prop_distrubution1(xx),'r-',label='$p_1$')
plt.legend()
plt.subplot(122)
plt.plot(xx,desired_prop_distrubution2(xx),'g-',label='$p_2$')
plt.legend()
plt.show()

Task

Implement a MCMC sampling using the Metropolis algorihtm to sample get samples from the two distributions given above. Also visualize the distribution functions and their samples in plots.

Hint

Use the below defined Gaussian as proposal density function. 

def proposal_density(mu,sigma=1): 
    return np.random.normal(mu,sigma)
def MCMC_Metropolis(proposal,start,distribution,size=1000):
    """
    args:
        proposal: porposal density function
        start:    start value for the porposal density function
        distribution: desired distribution function we want a sample from
        
    return:
        sample: returns a sample for our  desired distribution
    """
    
    raise NotImplementedError()

Literature

Licenses

Notebook License (CC-BY-SA 4.0)

The following license applies to the complete notebook, including code cells. It does however not apply to any referenced external media (e.g., images).

exercise-MCMC-Metropolis-sampling
Oliver Fischer
is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://gitlab.com/deep.TEACHING.

Code License (MIT)

The following license only applies to code cells of the notebook.

Copyright 2019 Oliver Fischer

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.