# 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)
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:

## 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))
```
```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: