# Code from Chapter 15 of Machine Learning: An Algorithmic Perspective (2nd Edition) # by Stephen Marsland (http://stephenmonika.net) # You are free to use, change, or redistribute the code in any way you wish for # non-commercial purposes, but please maintain the name of the original author. # This code comes with no warranty of any kind. # Stephen Marsland, 2008, 2014 # A simple Gibbs sampler import pylab as pl import numpy as np def pxgiveny(y,mx,my,s1,s2): return np.random.normal(mx + (y-my)/s2,s1) #return random.binomial(16,y,1) def pygivenx(x,mx,my,s1,s2): return np.random.normal(my + (x-mx)/s1,s2) #return random.beta(x+2,16-x+4,1) def gibbs(N=500): k=10 x0 = np.zeros(N,dtype=float) m1 = 10 m2 = 20 s1 = 2 s2 = 3 for i in range(N): y = np.random.rand(1) for j in range(k): x = pxgiveny(y,m1,m2,s1,s2) y = pygivenx(x,m1,m2,s1,s2) x0[i] = x return x0 #def f(x): # n = 16 # alph = 2 # bet = 4 # return 20.0*(np.factorial(n)/(np.factorial(x)*np.factorial(n-x)))*np.factorial(x+1)*np.factorial(19-x)/np.factorial(21) # #def factorial(n): # x = 1 # for i in range(n): # x *= (i+1) # return x def f(x): return np.exp(-(x-10)**2/10) N=500 s=gibbs(N) x1 = np.arange(0,17,1) pl.hist(s,bins=x1,fc='k') x1 = np.arange(0,17,0.1) px1 = np.zeros(len(x1)) for i in range(len(x1)): px1[i] = f(x1[i]) pl.plot(x1, px1*N*10/np.sum(px1), color='k',linewidth=3) pl.show()