Rotate 2D Gaussian given parameters a, b and c

[schniefen]

TL;DR Summary

Rotate mixtured Gaussian with variances unknown.

I have a Gaussian function of the form:

def f(x,y):
  a=some number
  b=...
  c=...
  return 3*np.exp(-a*(-0.5 + x)**2+b*(x-0.5)*(y-0.5)-c*(-0.5 + y)**2)

This is a Gaussian function symmetric around y=x, and I'd like to rotate it 45 degrees (counter)clockwise. Wikipedia gives an overdetermined system of equations for the variances of x and y respectively, but it looks cumbersome. Is there a simple way to do this?

Parameters a,b and c

 

[jambaugh]

The "best" way to view this is in linear algebra terms. The quadratic exponent in the Gaussian can be express via matrix representation of a quadratic form using (and defining) the covariance matrix.

$$f(x,y) = N \exp\left( -\frac{1}{2} (\mathbf{x}-\boldsymbol{\mu})^T K (\mathbf{x}-\boldsymbol{\mu})\right)$$
where, with your numbers...
$$\mathbf{x}=\left(\begin{array}{c} x \\ y\end{array}\right),\quad \boldsymbol{\mu}=\left(\begin{array}{c} \mu_x \\ \mu_y\end{array}\right)=\left(\begin{array}{c} 0.5\\ 0.5\end{array}\right), \text{ and }\quad K=\left(\begin{array}{cc} \sigma^2_x & \sigma_{xy}\\ \sigma_{xy} & \sigma^2_y\end{array}\right)=\left(\begin{array}{cc} 2a & b\\ b & 2c\end{array}\right)$$
The factors of 2 in the covariance matrix are due to your lack of a factored 1/2 in your code. It doesn't occur on the b terms because they will be added twice.

Then to rotate the distribution you rotate the covariance matrix (which is a rank 2 tensor) and the vector mean using a rotation matrix.

$$R_{45^\circ} = \frac{\sqrt{2}}{2}\left( \begin{array}{rr} 1 & -1\\ 1 & 1\end{array}\right),\quad \boldsymbol{\mu}' = R\boldsymbol{\mu}, \quad K' = RKR^T$$
With your numbers...
$$\boldsymbol{\mu}'=\frac{\sqrt{2}}{2}\left( \begin{array}{rr} 1 & -1\\ 1 & 1\end{array}\right)\cdot\left(\begin{array}{c} 0.5\\ 0.5\end{array}\right)=\frac{\sqrt{2}}{2}\left(\begin{array}{c} 0.5-0.5\\ 0.5+0.5\end{array}\right)=\left(\begin{array}{c} 0 \\ \sqrt{2}/2 \end{array}\right)$$
and
$$K' = \frac{2}{4}\left( \begin{array}{rr} 1 & -1\\ 1 & 1\end{array}\right)\left(\begin{array}{cc} 2a & b\\ b & 2c\end{array}\right)\left( \begin{array}{rr} 1 & 1\\ -1 & 1\end{array}\right)= \cdots =\left(\begin{array} ((a-b+c) & (a-c)\\(a-c) & (a+b+c)\end{array}\right)$$

It's "simpler" this way but that usually comes with more sophisticated approach. So brush up on your matrix multiplication.

Now this was an active rotation of the distribution ccw 45degrees. For the passive version, where you are rotating the axes use the inverse rotation matrix. For rotations this is just the transpose.
Also you should check my work (both to understand better and to catch any errors.)

In general a 2x2 rotation matrix takes the form:
$$R_\theta = \left(\begin{array}{cc} \cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array}\right)$$
Note that transposing has the same effect as reversing direction $\theta \mapsto -\theta$.

 

[schniefen]

If one would like to keep the Gaussian centered at 0.5,0.5, it is really just plugging in the numbers 0.5,0.5 into $\mu'$ ? Thanks for the helpful repetion on linear algebra.

 

[jambaugh]

Yes, basically just don't transform the means.

On a further note, you can dive deeper into this in principle component analysis. There you find a general linear transformation (not necessarily a rotation) which will make the covariance matrix diagional. This means you're resolving the distribution into independent components. In linear algebra terms you're finding the eigen-values and eigen-vectors. In probability terms you're determining the independent set of random variables.

 

[schniefen]

On a related note, consider the 2D Gaussian function below, for which a level curve has been sketched. Clearly it has a width and a length, and are either of these related to the variances? In 1D, the width is another name for the variance. In 2D, however, which is which? In this particular case I'm looking for the semi-minor axes of the ellipse.

 

[schniefen]

I don't know why, but somehow I don't get the Gaussian properly rotated. Consider the following code:

import numpy as np
import matplotlib.pyplot as plt
a=1.25
b=0   
c=10000
d=(a-b+c)/2
e=a-c
f=(a+b+c)/2
fig, ax = plt.subplots()
x,y=np.meshgrid(np.linspace(0,1,50),np.linspace(0,1,50))
z=3*np.exp(-a*(-0.5 + x)**2-b*(-0.5 + x)*(-0.5 + y)-c*(-0.5 + y)**2)
w=3*np.exp(-d*(-0.5 + x)**2-e*(-0.5 + x)*(-0.5 + y)-f*(-0.5 + y)**2) #rotated by 45 degrees
cs=ax.contour(x,y,z,levels=[0.8],colors='k',linestyles='dashed');
cs=ax.contour(x,y,w,levels=[0.8],colors='k',linestyles='dashed');
print(d,e,f)

Then I get 5000.625, -9998.75 , 5000.625, and the following plot:

Also, consider the integrals evaluated on WolframAlpha:

double integrate 3*exp(-1.25*(-0.5 + x)**2-0(x-0.5)(y-0.5)-10000(-0.5 + y)**2) from 0 to 1 - Wolfram|Alpha
double integrate 3*exp(-5000.625*(-0.5 + x)**2-9998.75(x-0.5)(y-0.5)-5000.625(-0.5 + y)**2) from 0 to 1 - Wolfram|Alpha

 

[schniefen]

How does one rotate the square $0 \leq x,y \leq 1$ by 45 degrees, over which the Gaussian is defined? And what are the new bounds for $x,y$?

 

[jambaugh]

To your last question, That's easily enough done geometrically. The origin point remains fixed, the opposite diagional (1,1) will rotate to $(0,\sqrt{2})$

As to your difficulties with the rotated gaussian, I'll check my math. BRB.

 

[schniefen]

Concerning my previous posts, I think the level curves are correct, it is just when one rotates the Gaussian, its central axis changes to be the diagonal of the square, and so more of the function fits into the square and the level curve changes correspondingly, right? However, the integral is still not correct, since the domain of integration for the rotated function should be the rotated square. What are the integration limits for that rotated square?

 

[jambaugh]

Concerning the rotated quadratic. I confirmed the matrix multiplication. It is easier to just ignore the usual 1/2 in the gaussian and rotate the quadratic form as is:
$$a' = a-b+c, \quad b'/2 = a-c, c'=a+b+c$$
So you're parameter transform looks good. Hmmm...
your c=10000 is huge that's a standard deviation of like 141.4 and you're "playing" in the unit square. Possibly both ellipses are off the grid but the pre-rotated one looks different due to truncation. Try c=1.

To integrate the rotated square you would normally just change coordinates, (effectively rotate back and integrate the original.) Or you can integrate over a suitable large standard square. You want the boundaries to be on the order of 10 standard deviations out from the mean. (I use 10 as a rule of thumb; it's easy to simply shift the decimal, but 6 should be sufficient.)