14  核密度函数图

14.1 核密度估计

核密度估计 (Kernel Density Estimation, KDE) 是一种用于估计未知概率密度函数的非参数方法，适用于连续型数据且不依赖于事先指定的分布形式。其基本思想是：在密度函数的每一个估计点上，根据样本点到该点的距离，使用核函数分配权重并加权平均，从而构建平滑的密度曲线。

设样本为 \(x_1, x_2, \dots, x_n\)，其密度函数在任意点 \(x\) 上的估计形式为：

\[ \hat{f}_h(x) = \frac{1}{n h} \sum_{i=1}^{n} K\left( \frac{x - x_i}{h} \right) \]

其中：

  • \(K(\cdot)\) 是核函数（kernel function），通常是一个对称的概率密度函数；
  • \(h > 0\) 是带宽参数（bandwidth），控制核函数的缩放程度和平滑水平；
  • \(\hat{f}_h(x)\) 是点 \(x\) 处的密度估计值。

14.1.1 核函数

在实际应用中，核函数的选择对估计结果的影响相对较小，而带宽的设置对估计曲线的光滑程度影响较大。

核函数的作用可以理解为：在估计点 \(x\) 处，根据样本点 \(x_i\) 与 \(x\) 之间的距离，赋予不同的权重。距离 \(x\) 越近的样本点，其权重越大；距离越远，权重越小。通过对所有样本点的加权平均，得到该点的密度估计。将所有位置的估计值拼接起来，即可得到整体的密度函数曲线。

为了更清楚地理解核函数的加权机制，我们可以对距离进行标准化处理，设：

\[ u_i = \frac{X_i - c}{h} \]

则以下两式等价：

\[ |u_i| \leq 1 \Longleftrightarrow |X_i - c| \leq h \]

记 \(D_i = |X_i - c|\)，表示第 \(i\) 个观察值与估计点 \(c\) 的距离。核函数的任务就是为每个 \(D_i\) 分配权重。

如下图所示，三种典型核函数的权重分配机制具有显著差异：

  • Uniform 核：在 \(|u| \leq 1\) 范围内赋予所有观察值相同的权重，超出范围的样本点权重为 0 (相当于弃之不用)。对应的密度估计不具有平滑性，常用于教学演示。
  • Triangle 核：采用线性下降的加权方式，距离估计点越近权重越大，边界处权重为 0，估计结果具有一定的连续性。
  • Epanechnikov 核：采用抛物线型权重函数，在 \(u=0\) 处取得最大值，具有最小均方误差（MSE）性质，估计曲线光滑、效率较高。
  • Gaussian 核：采用正态分布函数，所有样本点均有非零权重，平滑程度高，适用于大多数实际应用场景。

14.1.2 核函数的性质

常见核函数及其表达式：

  • Uniform 核函数 \(K(u) = \frac{1}{2} \cdot \mathbf{1}\{\left|u\right| \leq 1\}\) （也称为 Rectangular 核函数）

  • Triangle 核函数 \(K(u) = (1 - \left|u\right|) \cdot \mathbf{1}\{\left|u\right| \leq 1\}\)

  • Epanechnikov 核函数 \(K(u) = \frac{3}{4}(1 - u^2) \cdot \mathbf{1}\{\left|u\right| \leq 1\}\)

  • Quartic 核函数 \(K(u) = \frac{15}{16}(1 - u^2)^2 \cdot \mathbf{1}\{\left|u\right| \leq 1\}\)

  • Triweight 核函数 \(K(u) = \frac{35}{32}(1 - u^2)^3 \cdot \mathbf{1}\{\left|u\right| \leq 1\}\)

  • Gaussian 核函数 \(K(u) = \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{u^2}{2}\right)\)

  • Cosinus 核函数 \(K(u) = \frac{\pi}{4} \cos\left(\frac{\pi}{2} u\right) \cdot \mathbf{1}\{\left|u\right| \leq 1\}\)

核函数通常需要满足以下数学性质：

  1. 非负性：\(K(u) \geq 0\)
  2. 单位积分：\(\int_{-\infty}^{\infty} K(u) \, du = 1\)
  3. 对称性：\(K(u) = K(-u)\)
  4. 有限的二阶矩：\(\int u^2 K(u) \, du < \infty\)

实际使用中，还有一些细节需要注意。例如，部分文献或软件将 \(\mathbf{1}\{|u| \leq 1\}\) 写为 \(\mathbf{1}\{|u| < 1\}\)。对于连续变量，两者几乎没有区别；但若数据是离散型的（如整数型变量），则可能影响边界值是否被纳入计算。

核密度估计的构造可以理解为：以每一个样本点为中心放置一个缩放后的核函数，然后在每一个估计位置 \(x\) 上，取所有样本点的核值加权平均。因此，它是一种基于样本加权“局部贡献”的整体平滑过程。

总结而言：

  • 核函数定义了如何根据样本点与估计点之间的距离分配权重；
  • 带宽参数决定了每个样本点的影响范围；
  • 合理选择核函数和带宽参数是核密度估计中最关键的步骤；
  • 核密度估计为我们提供了一种平滑、灵活且无需模型假设的分布估计方法，广泛应用于经济学、金融学、机器学习等领域的探索性数据分析任务中。

14.2 单变量核密度函数图

14.3 多变量核密度函数图

  • 不同时期的收入分布 - 时序
  • 不同种族的收入分布 - 截面
  • 联合分布
#### 实例 1：山脊图

https://seaborn.pydata.org/examples/kde_ridgeplot.html
import requests
import pyreadstat

# 设置 User-Agent 模拟浏览器访问
url = "https://www.stata-press.com/data/r17/nlsw88.dta"
headers = {
    "User-Agent": "Mozilla/5.0"
}

# 下载文件并保存到本地
r = requests.get(url, headers=headers)
with open("data/nlsw88.dta", "wb") as f:
    f.write(r.content)

# 用 pyreadstat 读取本地文件
df, meta = pyreadstat.read_dta("data/nlsw88.dta")

# 显示前几行
print(df.head())
   idcode  age  race  married  never_married grade  collgrad  south  smsa  \
0       1   37     2        0              0    12         0      0     1   
1       2   37     2        0              0    12         0      0     1   
2       3   42     2        0              1    12         0      0     1   
3       4   43     1        1              0    17         1      0     1   
4       6   42     1        1              0    12         0      0     1   

   c_city industry occupation union       wage hours    ttl_exp     tenure  
0       0        5          6     1  11.739125    48  10.333334   5.333333  
1       1        4          5     1   6.400963    40  13.621795   5.250000  
2       1        4          3   NaN   5.016723    40  17.730770   1.250000  
3       0       11         13     1   9.033813    42  13.211537   1.750000  
4       0        4          6     0   8.083731    48  17.820513  17.750000  
import pandas as pd
import seaborn as sns

import matplotlib.pyplot as plt

# 读取 Stata 数据
url = "https://www.stata-press.com/data/r17/nlsw88.dta"
url = r"https://www.stata-press.com/data/r17/nlsw88.dta"
#data = pd.read_stata(url)

import pyreadstat

# 读取 Stata 文件
#data, meta = pyreadstat.read_dta(r"data\nlsw88.dta")

data, meta = pyreadstat.read_dta(url)



# 绘制山脊图
plt.figure(figsize=(10, 6))
sns.kdeplot(data=data, x="wage", hue="occupation", multiple="fill", common_norm=False)
plt.title("Ridge Plot of wage by Industry")
plt.xlabel("lnwage")
plt.ylabel("Density")
plt.show()

帮我读入 Stata 官网数据 webuse network, clear ，然后使用 seaborn 包绘制 lnwage 的山脊图，分组变量为 industry

目的：绘制多只股票 xxx 年日收益率的山脊图
数据来源：你选择合适的 package，适合中国大陆用户，自动在线下载
股票：中国移动，贵州茅台，万科A，比亚迪, 宁德时代, 南方航空, 格力电器
年份：2015 年 或 2024 年 或 2025 年 (可以指定一个 year 变量，有我随后根据需要填写)
import akshare as ak
import pandas as pd
import seaborn as sns

import matplotlib.pyplot as plt

# Define the stocks and year
stocks = ["中国移动", "贵州茅台", "万科A", "比亚迪", "宁德时代", "南方航空", "格力电器"]
year = 2024  # Replace with the desired year (e.g., 2015, 2024, or 2025)

# Initialize an empty DataFrame to store all stock data
all_data = pd.DataFrame()

# Fetch data for each stock
for stock in stocks:
    stock_data = ak.stock_zh_a_hist(symbol=stock, period="daily", start_date=f"{year}0101", end_date=f"{year}1231", adjust="qfq")
    stock_data["stock"] = stock  # Add a column for stock name
    # Inspect the column names to find the correct column for the closing price
    print(stock_data.columns)
    stock_data["daily_return"] = stock_data["Close"].pct_change()  # Calculate daily return
    all_data = pd.concat([all_data, stock_data], ignore_index=True)

# Filter the data for the specified year
all_data["date"] = pd.to_datetime(all_data["日期"])
all_data["year"] = all_data["date"].dt.year
filtered_data = all_data[all_data["year"] == year]

# Create a ridge plot
plt.figure(figsize=(10, 6))
sns.kdeplot(data=filtered_data, x="daily_return", hue="stock", multiple="fill", common_norm=False)
plt.title(f"Ridge Plot of Daily Returns for Selected Stocks in {year}")
plt.xlabel("Daily Return")
plt.ylabel("Density")
plt.show()

14.4 进阶用法

14.4.1 散点图 + 核密度函数图

penguins 数据集是一个经典的数据集，包含了三种企鹅物种（Adelie、Chinstrap 和 Gentoo）的生物测量数据。该数据集常用于统计分析和机器学习建模的教学和研究。本例中，我们抽取了三个变量：

  • species：企鹅的物种（Adelie、Chinstrap 或 Gentoo）。
  • flipper_length_mm：鳍的长度（单位：毫米）。
  • bill_length_mm：喙的长度（单位：毫米）。
  • body_mass_g：体重（单位：克）。

下图中，我们使用 seaborn 中的 jionplot 函数绘制了散点图和核密度函数图。散点图展示了每个物种的鳍长和喙长的分布情况，而核密度函数图则呈现了横轴变量 (flipper_length_mm) 和纵轴变量 (bill_length_mm) 的分布情况。

散点图反映了企鹅的喙长 (\(y\)) 与鳍长 (\(x\)) 之间的关系。可以看出，虽然三种企鹅的 \(y \sim x\) 之间都是正相关的，但 Chinstrap 类和 Gentoo 类企鹅的喙长和鳍长之间的正相关关系更强一些，而 Adelie 类企鹅的喙长和鳍长之间的正相关关系相对较弱。

密度函数图则揭示了单变量的分布特征。可以看出，平均而言，Adelie 类企鹅的喙长和鳍长都较小，而 Gentoo 类企鹅的喙长和鳍长都较大。Chinstrap 类企鹅的喙长 (\(y\)) 与 Gentoo 类企鹅的长度范围差不多，但平均而言，前者的喙更长一些；有趣的是，Chinstrap 类企鹅的鳍长 (\(x\)) 与 Adelie 类企鹅的长度范围差不多，但平均而言，前者的鳍更长一些。

虽然尚未看到三种企鹅的真实照片，但我们已经能够大致判断它们的体型特征了。

  • Adelie：体型较小，喙长和鳍长都较小。
  • Chinstrap：体型中等，喙长与 Gentoo 差不多；而鳍长则与 Adelie 相近。
  • Gentoo：体型较大，喙长和鳍长都较大。

简言之，Adelie 是「短嘴短翅」，Chinstrap 是「长嘴短翅」，Gentoo 是「长嘴长翅」。若进一步结合散点图来看，则可以推断出：Gentoo 类企鹅群体身材较为均匀，而 Adelie 类企鹅个体之间的身材差异较大。

# Source: https://seaborn.pydata.org/tutorial/introduction.html
import seaborn as sns
sns.jointplot(data=penguins, 
              x="flipper_length_mm", 
              y="bill_length_mm", 
              hue="species", 
              height=6)

它们到底长啥样？

你可以在 这里 找到它们的详细介绍。有位热心的网友特意绘制了它们的合影 (从左到右依次为 Adelie、Chinstrap 和 Gentoo)，如下图所示：

14.4.2 多个变量的情形：pairplot

在上面的例子中，我们只考虑了两个变量之间的关系。实际上，数据集中可能有多个变量，此时可以使用 seaborn 中的 pairplot 函数来绘制多个变量之间的关系图。

下面的例子中，我们同时呈现了三个变量 (flipper_length_mm、bill_length_mm 和 body_mass_g) 之间的两两配对散点图，以及单个变量的核密度函数图。与此同时，我们还使用了不同的颜色来区分不同的物种。

# Source：https://seaborn.pydata.org/tutorial/introduction.html
vlist = ["flipper_length_mm", "bill_length_mm", "body_mass_g"]
sns.pairplot(data=penguins, vars=vlist, 
             hue="species", 
             corner=True)

14.4.3 联合分布图

https://seaborn.pydata.org/generated/seaborn.kdeplot.html

14.4.4 3D 核密度函数图

Source: matplotlib - Fill under 3D line graphs

14.4.5 密度函数图+条形码

seaborn.rugplot

import seaborn as sns; sns.set_theme()
tips = sns.load_dataset("tips")
sns.kdeplot(data=tips, x="total_bill")
sns.rugplot(data=tips, x="total_bill")