Covariance matrix for circular variables?

Question

In my current project, I have a collection of three-dimensional samples such as [-0.5,-0.1,0.2]*pi, [0.8,-0.1,-0.4]*pi. These variables are circular/periodic, with their values ranging from -pi to +pi. It is my goal to calculate a 3-by-3 covariance matrix for these circular variables.

Python has an in-built function to calculate circular standard deviations, which I can use to calculate the standard deviations along each dimension, then use them to create a diagonal covariance matrix (i.e., without any correlation). Ideally, however, I would like to consider correlations between the parameters as well. Is there a way to calculate correlations between circular variables, or to directly compute the covariance matrix between them?

import numpy as np
import scipy.stats

# A collection of N circular samples
samples = np.asarray(
    [[0.384917, 1.28862,    -2.034],
    [0.384917,  1.28862,    -2.034],
    [0.759245,  1.16033,    -2.57942],
    [0.45797,   1.31103,    2.9846],
    [0.898047,  1.20955,    -3.02987],
    [1.25694,   1.74957,    2.46946],
    [1.02173,   1.26477,    1.83757],
    [1.22435,   1.62939,    1.99264]])

# Calculate the circular standard deviations
stds    = scipy.stats.circstd(samples, high = np.pi, low = -np.pi, axis = 0)

# Create a diagonal covariance matrix
cov     = np.identity(3)
np.fill_diagonal(cov,stds**2)

Answer 1

I had the same issue where my end goal was to perform a sort of "circular" PCA. In the end my approach was to use the first part from the formula to compute covariance between complex variables and a modeling of angles in the complex plane (Ref 1 & (Wikipedia))

Given your samples this is what I obtain:

import numpy as np
from scipy.stats import circvar, circstd, circmean
samples = np.asarray(
    [[0.384917, 1.28862,    -2.034],
    [0.384917,  1.28862,    -2.034],
    [0.759245,  1.16033,    -2.57942],
    [0.45797,   1.31103,    2.9846],
    [0.898047,  1.20955,    -3.02987],
    [1.25694,   1.74957,    2.46946],
    [1.02173,   1.26477,    1.83757],
    [1.22435,   1.62939,    1.99264]])

# Model our angles in the complex space
z = np.exp(1j*samples)
# Circular mean by feature
z_mean = np.mean(z, axis=0) 
# Computes the deviation to the mean by feature
z_deviation = z - z_mean
# Covariance matrix is computed using first term from formula from (1) and properties of matrix multiplication
cov_matrix = (z_deviation.conj().T @ z_deviation)/samples.shape[0] 
print(cov_matrix.real)

[[0.10913549 0.03437553 0.12631992] [0.03437553 0.03787767 0.00218626] [0.12631992 0.00218626 0.58049681]]

I cannot explain why but its diagonal is not equal from what you obtain with scipy circvar (which is not a built-in function) though their evolutions are strikingly similar (in the limits of small angles?)

print(cov_matrix.real.diagonal()/circvar(samples, high=np.pi, low=-np.pi, axis=0))

[1.94385619 1.98087835 1.64769066]

I lack the mathematical background and culture to explain this phenomenon but in my case (biology) this is just good enough. This approach is very consistent with what I obtain by performing a PCA on linearized sines and cosines of the angles and restores the interpretation of signs in the principal components.

Ref 1) Park, Kun Il (2018). Fundamentals of Probability and Stochastic Processes with Applications to Communications. Springer. ISBN 9783319680743.