from pylab import *

sigma = 4
mu = 3
N=1000
x = np.random.randn(N)*sigma + mu

H = hist(x, bins=N/30)

from scipy.stats import distributions

xx = np.linspace(-20,20)
yy = np.linspace(1, 10)
XX, YY = np.meshgrid(xx,yy)

@np.vectorize
def loglik(mu, sigma):
    return sum(distributions.norm.logpdf(x, mu, sigma))

contour(xx, yy, loglik(XX,YY), 200)

import scipy.optimize 

scipy.optimize.minimize(lambda x: -loglik(*x), [2,2])

np.mean(x)

np.std(x)