如何让Seaborn kdeplot生成类似R的双峰数据密度轮廓图？

我经常需要绘制二元数据的密度轮廓图，有时数据呈**双峰（bimodal）**分布。示例代码及问题如下：

1. 生成双峰分布数据 #

import numpy as np
x = np.concatenate([np.random.normal(size=1000, scale=.5), 
                    np.random.normal(size=100, loc=10, scale=.1)])
y = np.random.uniform(size=1100)
plt.scatter(x, y)

（散点图展示数据的双峰分布特征）

2. Seaborn kdeplot的过度平滑问题 #

使用Seaborn绘制二元密度轮廓时，默认参数得到的结果过度平滑，无法体现数据的双峰细节：

import seaborn as sns
import pandas as pd
plt.scatter(x, y)
sns.kdeplot(pd.DataFrame({'x': x, 'y':y}), x='x', y='y', c='orange', levels=[.1,.5])

（轮廓图仅显示过度模糊的密度层级）

同样，一维密度估计也存在类似问题：

plt.hist(x, density=True, bins=50)
sns.kdeplot(x)

（一维密度图过度平滑，双峰特征被掩盖）

3. R中的等效实现（符合预期） #

R中使用MASS::kde2d和density()可以得到清晰体现双峰细节的结果，代码如下：

x = c(rnorm(1000, sd=.5), rnorm(100, 10, .1))
y = runif(1100)
plot(x, y)

plot(x, y, cex=.5, xlim=c(-4, 12), ylim=c(-.2, 1.2))
contour(MASS::kde2d(x, y, n=100, lims=c(c(-4, 12), c(-.2, 1.2))), col='red', 
        add=T, levels = c(.1, .5))

hist(x, freq=F, breaks=50)
lines(density(x), col='red')

（结果能精准展示双峰分布的密度轮廓和细节）

需求：如何让Seaborn的kdeplot生成类似R的轮廓图？或Python中有其他更合适的成熟绘图工具？

一、调整Seaborn kdeplot的带宽参数 #

Seaborn默认使用Scott规则计算带宽，而R的density()和kde2d默认使用bw.nrd0规则，两者带宽计算逻辑不同导致平滑程度差异。可以手动指定带宽来贴近R的效果：

1. 一维密度估计 #

先实现R中bw.nrd0的等效带宽计算，再传入Seaborn：

import scipy.stats as stats
import matplotlib.pyplot as plt

# 实现R的bw.nrd0带宽计算逻辑
def bw_nrd0(x):
    x = np.asarray(x)
    if len(x) < 2:
        return 0
    hi = np.std(x, ddof=1)
    q75, q25 = np.percentile(x, [75, 25])
    lo = (q75 - q25) / 1.34
    return min(hi, lo) * len(x)**(-1/5)

# 计算带宽并绘图
bw = bw_nrd0(x)
plt.hist(x, density=True, bins=50)
# 用bw_adjust适配Seaborn的默认带宽基准
sns.kdeplot(x, bw_adjust=bw / stats.gaussian_kde(x).scotts_factor())
plt.show()

2. 二维密度轮廓图 #

分别指定x、y方向的R风格带宽：

import seaborn as sns
import pandas as pd

bw_x = bw_nrd0(x)
bw_y = bw_nrd0(y)

plt.scatter(x, y)
sns.kdeplot(
    pd.DataFrame({'x': x, 'y': y}),
    x='x', y='y',
    c='orange',
    levels=[.1, .5],
    bw_method=[bw_x, bw_y]  # 分别设置x、y方向的带宽
)
plt.xlim(-4, 12)
plt.ylim(-0.2, 1.2)
plt.show()

二、使用其他Python工具实现R风格的密度估计 #

1. 用scipy.stats.gaussian_kde自定义实现 #

Scipy的高斯核密度估计支持灵活调整带宽，完全复刻R的效果：

from scipy.stats import gaussian_kde
import matplotlib.pyplot as plt

# 构建二维核密度估计模型
kde = gaussian_kde(np.vstack([x, y]), bw_method='silverman')  # 用Silverman规则贴近R的默认逻辑

# 生成绘图网格
xgrid = np.linspace(-4, 12, 100)
ygrid = np.linspace(-0.2, 1.2, 100)
X, Y = np.meshgrid(xgrid, ygrid)
Z = kde(np.vstack([X.ravel(), Y.ravel()])).reshape(X.shape)

# 绘制散点+密度轮廓
plt.scatter(x, y, c='blue', s=10)
plt.contour(X, Y, Z, levels=[0.1, 0.5], colors='orange')
plt.xlim(-4, 12)
plt.ylim(-0.2, 1.2)
plt.show()

2. 用rpy2直接调用R的MASS包 #

如果想完全复用R的实现逻辑，可以通过rpy2调用R的kde2d：

import rpy2.robjects as robjects
from rpy2.robjects.packages import importr
import matplotlib.pyplot as plt
import numpy as np

MASS = importr('MASS')

# 将Python数据转换为R格式
r_x = robjects.FloatVector(x)
r_y = robjects.FloatVector(y)

# 调用R的kde2d函数
kde_result = MASS.kde2d(r_x, r_y, n=100, lims=robjects.FloatVector([-4,12,-0.2,1.2]))
xgrid = np.array(kde_result[0])
ygrid = np.array(kde_result[1])
Z = np.array(kde_result[2])

# 绘制结果
plt.scatter(x, y, s=10)
plt.contour(xgrid, ygrid, Z, levels=[0.1, 0.5], colors='red')
plt.xlim(-4,12)
plt.ylim(-0.2,1.2)
plt.show()

内容的提问来源于stack exchange，提问作者Devin F