scipy 使von Mises KDE适应海运

sdnqo3pr 于 2023-03-08 发布在其他

关注(0)|答案(2)|浏览(149)

我尝试使用Seaborn在极投影上绘制一个二元（联合）KDE。Seaborn不支持这个，Scipy也不直接支持角（von Mises）KDE。
scipy gaussian_kde and circular data解决了一个相关但不同的情况。相似之处是-随机变量被定义在单位圆上线性间隔的Angular 上;绘制KDE。差异：我想用Seaborn的joint kernel density estimate support来绘制这样的等高线图-

但没有分类（“物种”）变异，并在极投影上。边缘图将是很好的拥有，但并不重要。
我的情况的直线版本将是

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 °的Angular 之间也不应该有递减的轮廓线。
当然，正如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)

我从一个monkeypatching方法开始：

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变换，这是一个非常扭曲和低质量的预览：

尽管径向轴线将从中心向外增加而不是减小。

scipy

来源：https://stackoverflow.com/questions/75637853/adapt-von-mises-kde-to-seaborn

2条答案

按热度按时间

c9x0cxw01#

下面是一个方法的想法：

通过sin和cos将极坐标转换为笛卡尔坐标
使用这些来生成一个普通的jointplot（或kdeplot，这可以包括hue）
隐藏笛卡尔图的轴
要添加极面：在顶部创建虚拟极坐标图，并适当地重新缩放
可以手动添加Angular 和功率的kdeplot
要启用Angular 的环绕，只需复制偏移-360和+360度的Angular ，然后限制显示范围

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，并让Angular 顺时针旋转，则需要在调用2D kdeplot时切换x=和y=。

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')

赞(0）回复(0）举报 2023-03-08

x759pob22#

我过去使用seaborn和matplotlib的极投影的组合来完成这个任务。下面是一个例子：

import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

n_data = 1000
data_phase = np.random.rand(n_data) * 1.2 * np.pi
data_amp = np.random.randn(n_data)

fig = plt.figure()
ax = fig.add_subplot(111, projection='polar')

ax.scatter(data_phase, data_amp, vmin=0, vmax=2 * np.pi, s=10, alpha=0.3)
ax.set_thetagrids(angles=np.linspace(0, 360, 5));
sns.kdeplot(x=data_phase, y=data_amp, n_levels=5, c='k', ax=ax)

希望你能在那里适应你的需要？

赞(0）回复(0）举报 2023-03-08

我来回答

scipy 使von Mises KDE适应海运

2条答案

相关问题

热门标签

最新问答