Help with Methodology (vectors' angles from Kramers-Moyal coefficients) about kramersmoyal HOT 5 OPEN

Building on @LRydin's answer, note that a power [m, n] enters the formula, as in the paper's Data Analysis, as

D_{x^m y^n} = \int (x - x')^m (y - y')^n P(...)

So to get D_x and D_y, you'd need the power [1,0] and [0,1], respectively. Note that these have index 1 and 2 respectively if powers=[[0,0], [1,0], [0,1]]. Hence,

import numpy as np
import kramersmoyal

y = np.random.randn(10000, 2) # Some fake data
D, edges = kramersmoyal.km(y, powers=[[0,0], [1,0], [0,1]], bins=[20,10])

D_x = D[1, ...]
D_y = D[2, ...]

To obtain a proper velocity vector we should stack these

D = np.stack((D_x, D_y))

Note that this "2d-vector" is is defined for all points of the underlying grid {x, y} and hence has shape (2, bins[0], bins[1]). To get the angle, you'd need a correct function

def angle_between(v1, v2):
    norms = np.linalg.norm(v1) * np.linalg.norm(v2)

    # Avoid division by zero and numerical instabilities with clamp
    if norms == 0:
        return 0
    return np.arccos(np.clip(np.dot(v1, v2) / norms, -1.0, 1.0))

So now the paper states that you are computing the angle of this velocity vector with its neighbours. This is hard-coded for 2d but trivial to generalise for 3 dimensions.

def calculate_angles(v):
    _, n_x, n_y = v.shape
    angles = np.zeros((n_x, n_y))

    # Iterate over each point in the lattice
    for x in range(n_x):
        for y in range(n_y):
            v_current = v[:, x, y]

            # Pre-determine the valid neighbor offsets (luckily the underlying grid is box-like)
            neighbor_offsets = [(dx, dy) for dx in [-1, 0, 1] for dy in [-1, 0, 1] 
                                if 0 <= x + dx < n_x and 0 <= y + dy < n_y and not (dx == 0 and dy == 0)]

            # Calculate angles with valid neighbors and SUM THEM
            for dx, dy in neighbor_offsets:
                v_neighbor = v[:, x + dx, y + dy]
                angles[x, y] += angle_between(v_current, v_neighbor)

    return angles

Hence,

theta_xy = calculate_angles(D)

To visualise

import matplotlib.pyplot as plt
X, Y = np.meshgrid(edges[0], edges[1])
plt.pcolormesh(X, Y, theta_xy.T)
plt.show()

To now solve your problem, you'd need just to adapt the powers to 3d, that is powers = [[0,0,0], [1,0,0], [0,1,0], [0,0,1] to obtain (D_x, D_y, D_z) and then generalise calculate_angles to all the different angles in a higher dimensional space

from kramersmoyal.

Hey @athantas95,
Nice to know that you are interested in this. I skimmed through the paper and, from what I can understand what you have in mind should be possible - and not too much of a nightmare. Let's start with a 2D example because 3D is simply harder to visualise.

You should have access to a 2D time series of your dynamical (stochastic) system. Let's call it `y'. You should just have to use

import kramersmoyal as kmc
# Y is a (2,N) time series
powers = [[0,0], [1,0], [0,1]]
kmc, edges = kmc.km(y, powers=2)

Now you are interested in the kmc[1,...] and kmc[2,...] which correspond to the drift in the first and the first in the second dimension (here in our 2D problem). You will need something like quiver from matplotlib to visualise the results. I always get terribly confused with the way quiver works, but I think something like this should do:

import matplotlib.pyplot as plt
import numpy as np

X, Y = np.meshgrid(edges[0], edges[1])
fig, ax = plt.subplots()
ax.quiver(X, Y, kmc[1,...], kmc[2,...], units='width')

Can you give that a try?

from kramersmoyal.

It just occurred to me that angle_between may not be 100% suitable since it's just the angle magnitude, defined in [0, pi], and since we are adding angles, it would be better to use a formula valid for [0, 2pi). I found out this formula, which works in 2d

def angle_between(v1, v2):
    angle = np.arctan2(v2[1], v2[0]) - np.arctan2(v1[1], v1[0])
    return np.mod(angle, 2 * np.pi)

Similarly, since we are adding so many angles, it's important the final result is mod 2\pi

def calculate_angles(v):
  # ...
  return np.mod(angles, 2*np.pi)

Please take the script with a rock of salt, there might be other "problems" lingering, we thought we would just give you a head start. Good luck!

from kramersmoyal.

I only glanced at the paper and the formula for the angle should actually be arccos((u . v)/(|u| |v|)).

Do you have any minimal example you can share with us? Note that the bins are just the discretization, per dimension, of the underlying real-space -- that is, the space where your vectors are defined.

from kramersmoyal.

Hi again! Thanks a lot for the very quick responses and help. I tried the first snippet by @LRydin and I guess now I understand a bit better the output of the km function. I will now try the extended code by @fmeirinhos and let you know. Thanks a ton for your help :D

from kramersmoyal.