python - 使 von Mises KDE 适应 Seaborn

Question

我正在尝试使用 Seaborn 在极坐标投影上绘制双变量(联合)KDE。 Seaborn 不支持这一点，Scipy 也不直接支持 Angular (von Mises) KDE。

scipy gaussian_kde and circular data解决了一个相关但不同的案例。相似之处是 - 随机变量是在单位圆上线性间隔的角度上定义的； KDE 已绘制。区别:我想使用Seaborn的joint kernel density estimate support生成此类等高线图 -

但没有分类(“物种”)变化，并且是极坐标投影。边缘地 block 固然很好，但并不重要。

我的情况的直线版本是

import matplotlib
import pandas as pd
from matplotlib import pyplot as plt
import numpy as np
import seaborn as sns
from numpy.random._generator import default_rng

angle = np.repeat(
    np.deg2rad(
        np.arange(0, 360, 10)
    ),
    100,
)
rand = default_rng(seed=0)
data = pd.Series(
    rand.normal(loc=50, scale=10, size=angle.size),
    index=pd.Index(angle, name='angle'),
    name='power',
)

matplotlib.use(backend='TkAgg')
joint = sns.JointGrid(
    data.reset_index(),
    x='angle', y='power'
)
joint.plot_joint(sns.kdeplot, bw_adjust=0.7, linewidths=1)

plt.show()

但是这是在错误的投影中显示的，而且在 0 到 360 度之间不应该有递减的轮廓线。

当然，如Creating a circular density plot using matplotlib and seaborn解释说，在极坐标投影中使用现有高斯 KDE 的简单方法是无效的，即使我想也不能，因为 axisgrid.py 硬编码了没有参数的子图设置:

        f = plt.figure(figsize=(height, height))
        gs = plt.GridSpec(ratio + 1, ratio + 1)

        ax_joint = f.add_subplot(gs[1:, :-1])
        ax_marg_x = f.add_subplot(gs[0, :-1], sharex=ax_joint)
        ax_marg_y = f.add_subplot(gs[1:, -1], sharey=ax_joint)

我开始使用猴子补丁方法:

import scipy.stats._kde
import numpy as np


def von_mises_estimate(
    points: np.ndarray,
    values: np.ndarray,
    xi: np.ndarray,
    cho_cov: np.ndarray,
    dtype: np.dtype,
    real: int = 0
) -> np.ndarray:
    """
    Mimics the signature of gaussian_kernel_estimate
    https://github.com/scipy/scipy/blob/main/scipy/stats/_stats.pyx#L740
    """

    # https://stackoverflow.com/a/44783738
    # Will make this a parameter
    kappa = 20

    # I am unclear on how 'values' would be used here


class VonMisesKDE(scipy.stats._kde.gaussian_kde):
    def __call__(self, points: np.ndarray) -> np.ndarray:
        points = np.atleast_2d(points)

        result = von_mises_estimate(
            self.dataset.T,
            self.weights[:, None],
            points.T,
            self.inv_cov,
            points.dtype,
        )
        return result[:, 0]


import seaborn._statistics
seaborn._statistics.gaussian_kde = VonMisesKDE

这确实成功地被调用来代替默认的高斯函数，但是(1)它不完整，并且(2)我不清楚是否可以说服联合绘图方法使用新的投影。

通过 Gimp 转换得到的非常扭曲且低质量的预览:

虽然径向轴从中心向外会增加而不是减少。

Answer 1

这是一种方法的想法:

• 通过sin和cos将极坐标转换为笛卡尔坐标
• 使用它们生成正常的jointplot(或kdeplot，这可以包括hue)
• 隐藏笛卡尔图的轴
• 要添加极坐标方面:在顶部创建一个虚拟极坐标图并适本地重新缩放
• 可以手动添加角度和功率的 kdeplot
• 要启用角度环绕，只需复制移动 -360 度和 +360 度的角度，然后限制显示范围

from matplotlib import pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np

# test data from https://www.kaggle.com/datasets/muthuj7/weather-dataset
df = pd.read_csv('weatherhistory.csv')[['Wind Speed (km/h)', 'Wind Bearing (degrees)']].rename(
    columns={'Wind Bearing (degrees)': 'angle', 'Wind Speed (km/h)': 'power'})
df['angle'] = np.radians(df['angle'])

df['x'] = df['power'] * np.cos(df['angle'])
df['y'] = df['power'] * np.sin(df['angle'])

fig = plt.figure(figsize=(10, 10))
grid_ratio = 5
gs = plt.GridSpec(grid_ratio + 1, grid_ratio + 1)

ax_joint = fig.add_subplot(gs[1:, :-1])
ax_marg_x = fig.add_subplot(gs[0, :-1])
ax_marg_y = fig.add_subplot(gs[1:, -1])

sns.kdeplot(data=df, x='x', y='y', bw_adjust=0.7, linewidths=1, ax=ax_joint)

ax_joint.set_aspect('equal', adjustable='box')  # equal aspect ratio is needed for a polar plot
ax_joint.axis('off')
xmin, xmax = ax_joint.get_xlim()
xrange = max(-xmin, xmax)
ax_joint.set_xlim(-xrange, xrange)  # force 0 at center
ymin, ymax = ax_joint.get_ylim()
yrange = max(-ymin, ymax)
ax_joint.set_ylim(-yrange, yrange)  # force 0 at center

ax_polar = fig.add_subplot(projection='polar')
ax_polar.set_facecolor('none')  # make transparent
ax_polar.set_position(pos=ax_joint.get_position())
ax_polar.set_rlim(0, max(xrange, yrange))

# add kdeplot of power as marginal y
sns.kdeplot(y=df['power'], ax=ax_marg_y)
ax_marg_y.set_ylim(0, df['power'].max() * 1.1)
ax_marg_y.set_xlabel('')
ax_marg_y.set_ylabel('')
ax_marg_y.text(1, 0.5, 'power', transform=ax_marg_y.transAxes, ha='center', va='center')
sns.despine(ax=ax_marg_y, bottom=True)

# add kdeplot of angles as marginal x, repeat the angles shifted -360 and 360 degrees to enable wrap-around
angles = np.degrees(df['angle'])
angles_trippled = np.concatenate([angles - 360, angles, angles + 360])
sns.kdeplot(x=angles_trippled, ax=ax_marg_x)
ax_marg_x.set_xlim(0, 360)
ax_marg_x.set_xticks(np.arange(0, 361, 45))
ax_marg_x.set_xlabel('')
ax_marg_x.set_ylabel('')
ax_marg_x.text(0.5, 1, 'angle', transform=ax_marg_x.transAxes, ha='center', va='center')
sns.despine(ax=ax_marg_x, left=True)

plt.show()

PS:这是填充版本的样子(使用 cmap='turbo'):

如果您希望顶部为 0，并让角度顺时针旋转，则需要在对 2D 的调用中切换 x= 和 y= >kdeplot.

sns.kdeplot(data=df, x='y', y='x', bw_adjust=0.7, fill=True, cmap='turbo', ax=ax_joint)

# ....
ax_polar = fig.add_subplot(projection='polar')
ax_polar.set_theta_zero_location('N')
ax_polar.set_theta_direction('clockwise')