import sympy as S
S.var("x mu sigma", real=True)
Here's the likelihood for a 1-d normal distribution, in terms of two parameters: the mean $\mu$ and standard deviation $\sigma$
f = 1/(sqrt(2*S.pi)*sigma)*S.exp(-(x-mu)**2/(2*sigma**2))
f
And here's the log likelihood, which we can tidy up:
ll = log(f)
ll
ll = simplify(logcombine(ll))
ll
Let's pull out the derivatives of the log-likelihood: $\frac{\partial\mathcal(l)}{\partial\mu}$, $\frac{\partial\mathcal(l)}{\partial\sigma}$
ll.diff(mu)
ll.diff(sigma)
Second derivatives:
ll.diff(mu, mu)
ll.diff(mu, sigma)
ll.diff(sigma, sigma)
simplify(_)
Let's find the point where the gradient is zero for $\mu$ and $\sigma$
solve(f.diff(mu), mu)
solve(f.diff(sigma), sigma)
What if we want to use these calculations in numerical optimization code?
d2f_sigma2 = lambdify((mu, sigma, x), f.diff(sigma, sigma))
d2f_sigma2(1.2, 2, 3)
timeit d2f_sigma2(1.2, 2, 3)
1000000 loops, best of 3: 1.61 µs per loop
from sympy.utilities import codegen
codegen.ccode(f)
'(1.0L/2.0L)*sqrt(2)*exp(-1.0L/2.0L*pow(-mu + x, 2)/pow(sigma, 2))/(sqrt(M_PI)*sigma)'
print codegen.codegen(('d2f_dsigma2', f), 'C', 'd2f_dsigma2')[0][1]
/****************************************************************************** * Code generated with sympy 0.7.4.1 * * * * See http://www.sympy.org/ for more information. * * * * This file is part of 'project' * ******************************************************************************/ #include "d2f_dsigma2.h" #include <math.h> double d2f_dsigma2(double mu, double sigma, double x) { return (1.0L/2.0L)*sqrt(2)*exp(-1.0L/2.0L*pow(-mu + x, 2)/pow(sigma, 2))/(sqrt(M_PI)*sigma); }
Too lazy to compile a C extension. Do it for me!
from sympy.utilities import autowrap
fw = autowrap.autowrap(f.diff(sigma, sigma))
timeit fw(1.2, 2, 3)
1000000 loops, best of 3: 214 ns per loop
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)
1000000 loops, best of 3: 197 ns per loop
%timeit norm.logpdf(2.3, 6., 1)
10000 loops, best of 3: 86 µs per loop
x = np.arange(100)
vlogpdf = np.vectorize(logpdf)
timeit vlogpdf(2.3, 6, x)
10000 loops, best of 3: 51.8 µs per loop
timeit s_logpdf(2.3, 6., x)
10000 loops, best of 3: 108 µs per loop