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)