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)

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])

np.histogram(samples[:,0])

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)
    

?mlab.bivariate_normal