Keywords: discrete sampling | mcmc | metropolis | poisson distribution | proposal matrix | | Download Notebook
In simulated annealing, we carried out combinatorical oprimization by sampling from a state space where each state was a vector of baseball simulation features.
Since Metropolis MCMC is the same algorithm, it should be clear that its possible to simulate discrete possibilities in MCMC as long as you choose proposals which satisfy detailed balance.
As an example, consider simulating a poisson distribution. Since its discrete, the proposal wont be a continuous $q(x,y)$ (the proposal probability to go from y to x), but rather a matrix indexed by a variable that corresponds to (indexes) the various states that can be obtained.
def metropolis(p, qdraw, nsamp, xinit): samples=np.empty(nsamp) x_prev = xinit accepted=0 for i in range(nsamp): x_star = qdraw(x_prev) p_star = p(x_star) p_prev = p(x_prev) pdfratio = p_star/p_prev if np.random.uniform() < min(1, pdfratio): samples[i] = x_star x_prev = x_star accepted+=1 else:#we always get a sample samples[i]= x_prev return samples, accepted
Example: Sampling a Poisson
The poisson pmf is:
from scipy.stats import poisson xxx= np.arange(1,20,1) plt.plot(xxx, poisson.pmf(xxx, mu=5), 'o');
To sample from this distribution, we must create a proposal matrix which allows us to go from any integer output to any other in a finite number of steps. This matrix must be symmetric, since we wish to use Metropolis.
A simple such matrix, which is although a bit slow, would be one which has immediate off-diagonal elements (from Stats 580 at Iowa state..)
def prop_pdf(ito, ifrom): if ito == ifrom - 1: return 0.5 elif ito == ifrom + 1: return 0.5 elif ito == ifrom and ito == 0:#needed to make first row sum to 1 return 0.5 else: return 0
def prop_draw(ifrom): u = np.random.uniform() if ifrom !=0: if u < 1/2: ito = ifrom -1 else: ito = ifrom + 1 else: if u < 1/2: ito=0 else: ito=1 return ito
rv = poisson(5) samps, acc = metropolis(rv.pmf, prop_draw, 50000, 1)
xxx = np.arange(0,samps.max()) plt.hist(samps, bins=xxx, normed=True, align='left'); plt.plot(xxx, rv.pmf(xxx),'o');