import numpy as np
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
%matplotlib inline
def p(x, y):
g1 = mlab.bivariate_normal(x, y, 1.0, 1.0, -1, -1, -0.8)
g2 = mlab.bivariate_normal(x, y, 1.5, 0.8, 1, 2, 0.6)
return 0.6*g1+28.4*g2/(0.6+28.4)
def q(r):
return r + np.random.normal(size=2)
N = 100000
s = 10 #thinning
r = np.zeros(2)
p0 = p(r[0], r[1])
print p0
samples = []
for i in xrange(N):
rn = q(r)
pn = p(rn[0], rn[1])
if pn >= p0:
p0 = pn
r = rn
else:
u = np.random.rand()
if u < pn/p0:
p0 = pn
r = rn
if i % s == 0:
samples.append(r)
samples = np.array(samples)
0.00632453617818
plt.scatter(samples[:, 0], samples[:, 1], alpha=0.5, s=1)
'''Plot target'''
dx = 0.01
x = np.arange(np.min(samples), np.max(samples), dx)
y = np.arange(np.min(samples), np.max(samples), dx)
X, Y = np.meshgrid(x, y)
Z = p(X, Y)
CS = plt.contour(X, Y, Z, 10)
plt.clabel(CS, inline=1, fontsize=10)
plt.show()
print np.mean(samples[:,0]), np.mean(samples[:,1])
0.210361197903 0.826203546417
np.histogram(samples[:,0])
(array([ 24, 251, 1463, 2693, 2517, 1821, 929, 256, 41, 5]), array([-4.90450004, -3.72145056, -2.53840108, -1.35535161, -0.17230213, 1.01074735, 2.19379683, 3.37684631, 4.55989579, 5.74294527, 6.92599475]))
plt.hist(samples[:,0], bins=80, histtype='step',normed=True, color='b', label='x')
plt.hist(samples[:,1], bins=80, histtype='step', normed=True, color='r', alpha=0.5, label='y')
plt.title("Marginals")
plt.xlabel("Value")
plt.ylabel("Probability")
plt.legend()
plt.show()
from sklearn import mixture
def fit_samples(samples):
gmix = mixture.GMM(n_components=2, covariance_type='full')
gmix.fit(samples)
print gmix.means_
colors = ['r' if i==0 else 'g' for i in gmix.predict(samples)]
ax = plt.gca()
ax.scatter(samples[:,0], samples[:,1], c=colors, alpha=0.8)
plt.show()
fit_samples(samples)
[[-0.96433308 -1.03744208] [ 0.93701863 1.97904132]]
?mlab.bivariate_normal