Demo iagengauss

Illustrate the generation of a d-dimensional Gaussian image

Description

The sequence below shows a technique to a d-dimensional Gaussian image, understanding the difficulties in computing an equation with vector and matrix notation.

Script

1 dimensional case

The Gaussian function is a symmetric bell shaped function that is characterized by two parameters: mean and variance. The one-dimensional Gaussian function at point is given by the following equation, with mean and variance . The function is maximum at point and it falls by the factor (aprox. 0.6) at point away from the mean.

Equation

As this function is scalar, it is possible to compute this function on M samples represented as a 1 x M vector :

Code

 1 from numpy import *
 2 from numpy.linalg import *
 3 import ia636 as ia
 4 
 5 # First case: unidimensional
 6 # x: single value (single sample) or a row of values (many samples)
 7 # mu and sigma are scalar
 8 def fun1(x, mu, sigma):
 9   d = x - mu
10   k = (1. * d * d)/ (sigma * sigma)
11   return (1./(sqrt(2 * pi) * sigma)) * exp(-1./2 * k)
12 
13 print 'Computing the Gaussian function at a single point'
14 ex1 = "fun1( 10, 10, 5)"
15 print ex1,"=>", eval(ex1)
16 ex2 = "fun1( 15, 10, 5)"
17 print ex2,"=>", eval(ex2)
18 
19 print '\nComputing the Gaussian function at many points, using the same code'
20 ex3 = "fun1( array([10,15,20]), 10, 5)"
21 print ex3,"=>", eval(ex3)
22 
23 x = arange(-5,26)
24 y = fun1(x, 10, 5)
25 adshow(ia.iaplot([y],[x]))
26 print fun1(array([-4.,-3.0,-2.0,-1.0,0,1.0,2.0,3.0,4.0]),0,50)
Computing the Gaussian function at a single point
fun1( 10, 10, 5) => 0.0797884560803
fun1( 15, 10, 5) => 0.0483941449038

Computing the Gaussian function at many points, using the same code
fun1( array([10,15,20]), 10, 5) => [ 0.07978846  0.04839414  0.01079819]
[ 0.00795335  0.0079645   0.00797247  0.00797725  0.00797885  0.00797725
  0.00797247  0.0079645   0.00795335]

d-dimensional Case

If a sample point is a vector of dimension d: , the d-dimensional Gaussian function is characterized by the mean vector: and the symmetric square covariance matrix:

Equation

Code

 1 # Second case: d-dimensional, single sample
 2 # x: single column vector (single sample with d characteristics)
 3 # mu: column vector, 1 x d
 4 # sigma: covariance matrix, square and symmetric, d x d
 5 def funn(X, MU, COV):
 6   d = len(X)
 7   Xc = X - MU
 8   k = 1. * dot( dot(transpose(Xc), inv(COV)), Xc)
 9   return (1./((2 * pi)**(d/2.) * sqrt(det(COV)))) * exp(-1./2 * k)
10 
11 print '\ncomputing the Gaussian function at a single 3-D sample'
12 X1 = array([[10],
13             [5],
14             [3]])
15 MU = X1
16 COV = array([[10*10, 0,   0],
17              [0,     5*5, 0],
18              [0,     0,   3*3]])
19 print 'X1=',X1
20 print 'MU=',MU
21 print 'COV=',COV
22 ex4 = "funn( X1, MU, COV)"
23 print ex4,"=>", eval(ex4)
24 
25 print '\nComputing the Gaussian function at two 3-D samples'
26 print '\nNote that it does not work'
27 X2 = 1. * X1/2
28 X = hstack([X1,X2])
29 print 'X=',X
30 ex5 = "funn( X, MU, COV)"
31 print ex5,"=>", eval(ex5)
computing the Gaussian function at a single 3-D sample
X1= [[10]
 [ 5]
 [ 3]]
MU= [[10]
 [ 5]
 [ 3]]
COV= [[100   0   0]
 [  0  25   0]
 [  0   0   9]]
funn( X1, MU, COV) => [[ 0.00042329]]

Computing the Gaussian function at two 3-D samples

Note that it does not work
X= [[ 10.    5. ]
 [  5.    2.5]
 [  3.    1.5]]
funn( X, MU, COV) => [[ 0.00042329  0.00042329]
 [ 0.00042329  0.00029092]]

Computing d-dimensional Gaussian function on n sample points directly

The exponent part of the d-dimensional equation is an inner product with the covariance matrix in the center. When the data is arranged as a block matrix of n d-dimensional points, we need to apply the inner product to each d-dimensional point. This is equivalent to use only the diagonal results of the matrix product. More information can be seen at Matrix Mulplication wikipedia page.

 1 # Third case: m n-dimensional computing
 2 # X: n column vectors (n samples with d characteristics)
 3 # MU: column vector, 1 x M
 4 # COV: covariance matrix, square and symmetric, d x d
 5 def funm(X, MU, COV):
 6   d = len(MU)
 7   Xc = X - MU
 8   k = 1. * diagonal(dot(transpose(Xc), dot(inv(COV), Xc)))
 9   return (1./((2 * pi)**(d/2.) * sqrt(det(COV)))) * exp(-1./2 * k)
10 
11 print '\ncomputing the Gaussian function on two 3-D samples'
12 X = array([[10, 5],
13            [ 5, 2.5],
14            [ 3, 1.5]])
15 MU = transpose(array([[10, 5, 3]]))
16 COV = array([[10*10, 0,   0],
17              [0,     5*5, 0],
18              [0,     0,   3*3]])
19 print 'X=',X
20 print 'MU=',MU
21 print 'COV=',COV
22 ex6 = "funm( X, MU, COV)"
23 print ex6,"=>", eval(ex6)
computing the Gaussian function on two 3-D samples
X= [[ 10.    5. ]
 [  5.    2.5]
 [  3.    1.5]]
MU= [[10]
 [ 5]
 [ 3]]
COV= [[100   0   0]
 [  0  25   0]
 [  0   0   9]]
funm( X, MU, COV) => [ 0.00042329  0.00029092]
 1 from ia636 import ianormalize
 2 # Forth case: optimized m n-dimensional computing
 3 # X: n column vectors (n samples with d characteristics)
 4 # MU: column vector, 1 x M
 5 # COV: covariance matrix, square and symmetric, d x d
 6 def funm1(X, MU, COV):
 7   d = len(MU)
 8   Xc = X - MU
 9   k = 1. * Xc * dot(inv(COV), Xc)
10   k = sum(k,axis=0) #the sum is only applied to the rows
11   return (1./((2 * pi)**(d/2.) * sqrt(det(COV)))) * exp(-1./2 * k)
12 
13 print '\ncomputing the Gaussian function on two 3-D samples'
14 X = array([[10, 5],
15            [ 5, 2.5],
16            [ 3, 1.5]])
17 MU = transpose(array([[10, 5, 3]]))
18 COV = array([[10*10, 0,   0],
19              [0,     5*5, 0],
20              [0,     0,   3*3]])
21 print 'X=',X
22 print 'MU=',MU
23 print 'COV=',COV
24 ex6 = "funm1( X, MU, COV)"
25 print ex6,"=>", eval(ex6)
26 
27 i,j = indices((50,100))
28 x = vstack((ravel(i),ravel(j)))
29 MU = transpose(array([[25, 40]]))
30 COV = array([[15*15, 0, ],
31              [0,     10*10]])
32 y = funm1(x, MU, COV).reshape((50,100))
33 adshow(ianormalize(y).astype(uint8))
computing the Gaussian function on two 3-D samples
X= [[ 10.    5. ]
 [  5.    2.5]
 [  3.    1.5]]
MU= [[10]
 [ 5]
 [ 3]]
COV= [[100   0   0]
 [  0  25   0]
 [  0   0   9]]
funm1( X, MU, COV) => [ 0.00042329  0.00029092]

Suggested Exercises

  1. Compute the 2 dimensional Gaussian image using the equation of a decomposible Gaussian below. You can use the meshgrid style programming where the image (coordinates and value) can be stored in an image like format.
  1. Show that the previous 2 dimensional Gaussian image can be composed as an outer product of two one dimensional signals
  2. Find the equation of a 2 dimensional Gaussian image when the covariance matrix is of the form: