1# Copyright Jim Bosch 2010-2012. 2# Distributed under the Boost Software License, Version 1.0. 3# (See accompanying file LICENSE_1_0.txt or copy at 4# http://www.boost.org/LICENSE_1_0.txt) 5 6import numpy 7import gaussian 8 9mu = numpy.zeros(2, dtype=float) 10sigma = numpy.identity(2, dtype=float) 11sigma[0, 1] = 0.15 12sigma[1, 0] = 0.15 13 14g = gaussian.bivariate_gaussian(mu, sigma) 15 16r = numpy.linspace(-40, 40, 1001) 17x, y = numpy.meshgrid(r, r) 18 19z = g(x, y) 20 21s = z.sum() * (r[1] - r[0])**2 22print "sum (should be ~ 1):", s 23 24xc = (z * x).sum() / z.sum() 25print "x centroid (should be ~ %f): %f" % (mu[0], xc) 26 27yc = (z * y).sum() / z.sum() 28print "y centroid (should be ~ %f): %f" % (mu[1], yc) 29 30xx = (z * (x - xc)**2).sum() / z.sum() 31print "xx moment (should be ~ %f): %f" % (sigma[0,0], xx) 32 33yy = (z * (y - yc)**2).sum() / z.sum() 34print "yy moment (should be ~ %f): %f" % (sigma[1,1], yy) 35 36xy = 0.5 * (z * (x - xc) * (y - yc)).sum() / z.sum() 37print "xy moment (should be ~ %f): %f" % (sigma[0,1], xy) 38