import sympy as S
S.var("x mu sigma", real=True)


f = 1/(sqrt(2*S.pi)*sigma)*S.exp(-(x-mu)**2/(2*sigma**2))
f


ll = log(f)
ll

ll = simplify(logcombine(ll))
ll

ll.diff(mu)


ll.diff(sigma)


ll.diff(mu, mu)


ll.diff(mu, sigma)

ll.diff(sigma, sigma)

simplify(_)

solve(f.diff(mu), mu)

solve(f.diff(sigma), sigma)

d2f_sigma2 = lambdify((mu, sigma, x), f.diff(sigma, sigma))

d2f_sigma2(1.2, 2, 3)

timeit d2f_sigma2(1.2, 2, 3)

from sympy.utilities import codegen

codegen.ccode(f)

 print codegen.codegen(('d2f_dsigma2', f), 'C', 'd2f_dsigma2')[0][1]

from sympy.utilities import autowrap

fw = autowrap.autowrap(f.diff(sigma, sigma))

timeit fw(1.2, 2, 3)

logpdf = autowrap.autowrap(ll)

logpdf(2.3, 6., 1)

from scipy.stats.distributions import norm
s_logpdf = norm.logpdf

s_logpdf(1, 2.3, 6.)

%timeit logpdf(2.3, 6., 1)

%timeit norm.logpdf(2.3, 6., 1)

x = np.arange(100)
vlogpdf = np.vectorize(logpdf)

timeit vlogpdf(2.3, 6, x)

timeit s_logpdf(2.3, 6., x)