Plotting single points of a covariance ellipse using a Lissajous principle

covariance_ellipse_as_lissajous_curve.py

import numpy as np
from numpy import cos, sin
import matplotlib.pyplot as plt

"""
Looking at the covariance ellipse as a Lissajous figure.

In this gist: https://gist.github.com/CarstenSchelp/b992645537660bda692f218b562d0712
I plot a 2-D covariance ellipse by 
- normalizing the covariance matrix to contain only two ones, on the diagonal, and a pearson coefficient
- creating the ellipse as a matplotlib patch from that normalized covariance matrix
- transforming the ellipse:
    - rotate by 45 degrees
    - scale it back (reverting the normalization)
    - translate it to its original center point
    
 However, there is still an alternative method that lends itself well to generate single
 points instead of a matplotlib-patch.
 Remember Lissajous curves? https://en.wikipedia.org/wiki/Lissajous_curve
 The curve where the frequencies for x and y are exactly the same and the only difference is
 a shift in phase, is an ellipse.
 It "hangs in the graph" at a slant (of 45 degrees when both amplitudes are the same, too).
 
 You can use the code from the above mentioned gist.
 From line 45 on you can insert something like the following code:
 
 - Use the pearson-coefficient to calculate the phase-angle
 - choose some angles to apply the phase to (the array "phi" in this example)
 - calculate x and y
 - rescale x and y (not given in the code example here)
 - plot the points [x, y]
 
 You can detect the angles of the vertices of the ellipse by using the formula I give in the code
 ("phi_extreme")
 
 My apologies, these days I do not have the time to explain how I came to this.
 It will have to wait for another moment.
 But by all means, feel free to write a comment to this gist if you are curious or if
 you want to discuss things.
"""

# The pearson coefficient (only one in 2D).
# Put any value from -1 to 1,here:
pearson = .7
phase = np.arcsin(pearson)

# creating 360 different values of angles
# between 0 and 2*PI. You can create less,
# or more and different ones.
phi = np.linspace(0, 2*np.pi, num=360)
x = cos(phi)
y = sin(phi + phase)
# plot the "one stdev box"
plt.hlines((-1, 1), xmin=-1, xmax=1)
plt.vlines((-1, 1), ymin=-1, ymax=1)
# plot the lines that run through the
# tangent points (box-ellipse)
plt.hlines((pearson,-pearson), xmin=-1, xmax=1, color='gray')
plt.vlines((pearson,-pearson), ymin=-1, ymax=1, color='gray')
plt.axis('equal')
plt.plot(x, y)

# Show: local maxima of r are at 45deg and 135deg.
phi_extreme = np.arctan((1 - sin(phase))/cos(phase))
plt.scatter(cos(phi_extreme), sin(phi_extreme+phase), c='orange')
plt.scatter(-sin(phi_extreme), cos(phi_extreme+phase), c='purple')
plt.show()