数据的核密度估计及其可视化：Python实现
文章目录
- 实例数据介绍
- 精神疾病治疗时长数据表
- 间歇泉爆发持续时间数据表
- 文件放置介绍
- 频率分布直方图法
- 原理探讨与分析
- Python 代码实现
- 简单核密度估计法
- 原理探讨与分析
- Python代码实现
- 核密度估计法
- 原理探讨与分析
- Python代码实现
- 最小近邻方法的核密度估计
- 原理探讨与分析
- Python 代码实现
- 动态核密度方法
- 原理探讨与分析
- Python 代码实现
本文的参考文献主要如下：
- Silverman B.W. Density Estimation. Chapman and Hall, London, 1986
- 标准 ISO 13528
实例数据介绍
本博客采用的示例数据为本博客的参考文献一中的表 2.1 和 2.2，如下所示：
精神疾病治疗时长数据表
间歇泉爆发持续时间数据表
文件放置介绍
我们将上述的数据表放置在代码文件夹，所在路径的另一个命名为 “数据” 的文件夹中：
代码：
-------> xxx.py (代码)
数据：
-------> 间歇泉爆发持续时间表
-------> 自杀倾向治疗时间
频率分布直方图法
原理探讨与分析
设数据集为 X=(X1,X2,⋯,Xn)X=(X_1, X_2, \cdots, X_n)X=(X1,X2,⋯,Xn)：
频率分布直方图需要我们提供画图原点 x0x_0x0 以及一个带宽 hhh，将落在 [x0+(m−1)⋅h,x0+m⋅h][x_0 + (m-1)\cdot h, x_0 + m \cdot h][x0+(m−1)⋅h,x0+m⋅h] 上的点的数量进行计数，后除以带宽长度即可：
f^(x)=1nh(落在与x同一带宽上的Xi的数量)\hat{f}(x) = \frac{1}{nh} (\text{落在与} x \text{同一带宽上的} X_i \text{的数量}) f^(x)=nh1(落在与x同一带宽上的Xi的数量)
一般的在对数据进行分析时,常用频率分布直方图，对数据进行初步的可视化，我们将精神疾病治疗时长数据表，与间歇泉爆发持续时间数据表进行频率分布直方图画图，如下所示：
从上图可以看出，即便带宽相同，直方图的原点不同会导致直方图具有不同的视觉效果。对于间接泉的爆发持续时间直方图来说，虽然第1个图和第2个图都表现出明显的双峰分布，但是对于图2来说双峰的间隔点难以看出。
对于精神疾病治疗时长的直方图来说，第2个图在区 [150,250][150, 250][150,250] 时，可以看出一个次峰。对于非统计学者来说，很可能会造成数据呈现双峰分布的误解。毕竟从图1来看这很可能是由随机性导致的。
Python 代码实现
# -*- coding: utf-8 -*-
"""
Created on Sun Sep 19 10:11:15 2021@author: zhuozb本代码主要展示了 13528 参考文献 30 里的第二节中，有关 histogram 的用例讨论。本代码主要讨论了间歇泉爆发性持续时间，以及自杀精神治疗时间的频数分布画图。
"""import matplotlib.pyplot as plt
import numpy as np
import pandas as pddef extract_data(path):# 将Excel数据表读入Python中table = pd.read_excel(path, header=None)# 将数据表转换为 2D 数列data = table.values# 剔除非数字数据并将一维数列转换为二维数列data = data[~np.isnan(data)]return datadef plot_histogram(data, bins='auto', title=''):# fig, ax = plt.subplots()# 直接使用 plt.hist 画图freq, bins, patches = plt.hist(data, bins=bins,density=True,edgecolor='black')# 频数分布直方图的标签位置bin_centers = np.diff(bins)*0.5 + bins[:-1]for fr, x, patch in zip(freq, bin_centers, patches):height = round(fr, 4)plt.annotate("{}".format(height),xy = (x, height), # top left corner of the histogram barxytext = (0,0.2), # offsetting label position above its bartextcoords = "offset points", # Offset (in points) from the *xy* valueha = 'center', va = 'bottom',fontfamily='arial') # 还用a字体展示标签# 自动的调整X轴和Y轴的显示范围，虽然没有什么用ax.set_xlim(auto=True)ax.set_ylim(auto=True)# 设置X轴和Y轴的标签if title != '':# 如果没有给出标题，则不显示X轴的标签plt.xlabel(title, fontfamily='simhei')plt.ylabel('频数', fontfamily='simhei')plt.grid(axis='y')if __name__ == '__main__':eruption_lenght_path = r'../数据/间歇泉爆发持续时间表.xlsx'treatment_spell_path = r'../数据/自杀倾向治疗时间.xlsx'# 读取数据集# 读取数据表格 2.1treatment_spell_data = extract_data(treatment_spell_path)# 画出表格2.2eruption_lenght_data = extract_data(eruption_lenght_path)# 画出图2.1plot_histogram(eruption_lenght_data, bins=np.arange(1.25, 5.5, 0.5), title='爆发持续时间')plot_histogram(eruption_lenght_data, bins=np.arange(1.5, 5.5, 0.5), title='爆发持续时间')# 画出图2.2plot_histogram(treatment_spell_data, bins=np.arange(0, 1000, 50), title='自杀精神治疗时长')plot_histogram(treatment_spell_data, bins=np.arange(-25, 1000, 50), title='自杀精神治疗时长')
简单核密度估计法
原理探讨与分析
简单核密度的估计法的原理是将频率分布直方图的带宽设置为接近于0，对于一个概率密度函数来说都有：
f(x)=limh→012hP(x−h<X<x+h)f(x) = \lim_{h\to0} \frac{1}{2h} P(x-h < X < x+h) f(x)=h→0lim2h1P(x−h<X<x+h)
于是我们可以用带宽为无限小的频数分布图来估计概率密度：
f^(x)=12nh(落在带宽(x−h,x+h)上的Xi的数量)\hat{f}(x) = \frac{1}{2nh} (\text{落在带宽} (x-h, x+h) \text{上的} X_i \text{的数量}) f^(x)=2nh1(落在带宽(x−h,x+h)上的Xi的数量)
为了方便表述我们定义一个核密度函数如下所示：
w(x)={12if ∣x∣<10otherwise. w(x)= \begin{cases}\frac{1}{2} & \text { if }|x|<1 \\ 0 & \text { otherwise. }\end{cases} w(x)={210 if ∣x∣<1 otherwise.
于是，f^(x)\hat{f}(x)f^(x) 可以表示为：
f^(x)=1hn∑i=1nw(x−Xih)\hat{f}(x) = \frac{1}{hn} \sum_{i=1}^{n} w ( \frac{x-X_i}{h} ) f^(x)=hn1i=1∑nw(hx−Xi)
在实际应用中带宽当然可以不需要取到无限接近于0，以间歇泉爆发持续时间为例，我们可以将带宽设置为0.25画出的简单核密度图如下所示：
由于我们的经验密度函数f^(x)\hat{f}(x)f^(x)不是一个连续函数，在中断点处是不可导的。从上图也可以直观的看出，这种图像对于未经学习的学者来说，图像中的阶跃部分可能具有一定的迷惑性。
Python代码实现
# -*- coding: utf-8 -*-
"""
Created on Sun Sep 19 11:37:22 2021@author: zhuozb这是ISO13528中的参考文献，30里面的2.3小节的案例画图
"""import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from histogram_example_showing import extract_datadef digits_after_decimal(num):'''返回变量 num 的小数点位数'''num = str(num)# 小数点所在位置digit_loc = num.find('.')# 小数点后的位数digits_decimal = len(num) - digit_loc - 1return digits_decimaldef plot_naive_estimator(data, bandwidth=0.25, title=''):'''data ：待统计的数据bandwidth：画图带宽title：画出图像的标题，如果不给则默认不写出标题'''fig, ax = plt.subplots()# n 为样本容量n = len(data)# 返回数列中最小的数，并根据这个数来产生（极限）间隔范围min_num = min(data)# 找出最小数的小数点后的位数digits_decimal = digits_after_decimal(min_num)# 画图点间隔data_interval = 0.1**(digits_decimal) # 求出数列中的最大数，以此来产生画图范围max_num = max(data)# 产生X轴的画图点xtick_array = np.arange(0.9*min_num if (min_num>0) else 1.1*min_num,1.1*max_num, data_interval)# 产生Y轴的画图点ytick_array = []def calculate_y(x, x_data, bandwidth):# 根据公式2.1求解，X对应的Yy_sum = 0for x_i in data:# 根据 (x-X)/h 从数据数列中产生一个权重数列weight_array = np.abs((x - x_i)/bandwidth)# 根据公式2.1, 将权重数列按条件求和sum_weight = len(weight_array[(weight_array < 1)])/2y_sum = y_sum + sum_weight# y = y_sum/(n*bandwidth)return y ytick_array = [calculate_y(x, data, bandwidth) for x in xtick_array]# 设置X轴和Y轴的标签if title != '':# 如果没有给出标题，则不显示X轴的标签plt.xlabel(title, fontfamily='simhei')plt.ylabel('概率密度估计', fontfamily='simhei') plt.plot(xtick_array, ytick_array, linewidth=1)# 设置网格线plt.grid(axis='both')if __name__ == '__main__':eruption_lenght_path = r'../数据/间歇泉爆发持续时间表.xlsx'treatment_spell_path = r'../数据/自杀倾向治疗时间.xlsx'# 读取数据集# 读取数据表格 2.1treatment_spell_data = extract_data(treatment_spell_path)# 画出表格2.2eruption_lenght_data = extract_data(eruption_lenght_path) plot_naive_estimator(eruption_lenght_data, title=f'爆发持续时间,h = 0.25')
核密度估计法
原理探讨与分析
如上述的简单，核密度估计一般使用核密度估计时，首先需要定义一个核密度函数 K(x)K(x)K(x)，核密度函数一是一个概率密度，他要求即在坐标轴上的积分面积为单位面积，一般的可以取标准正态分布作为核密度函数：
K(x)=12πexp(−x22)K(x) = \frac{1}{\sqrt{2\pi}} \exp(-\frac{x^2}{2}) K(x)=2π1exp(−2x2)
并定义经验密度函数f^(x)\hat{f}(x)f^(x)如下：
f^(x)=1nh∑i=1nK(x−Xih)\hat{f}(x) = \frac{1}{nh} \sum_{i=1}^{n} K ( \frac{x-X_i}{h} ) f^(x)=nh1i=1∑nK(hx−Xi)
在核密度估计中，一般把参数 h 称为窗口宽度，窗口宽度的选择将影响到数据的可视化效果，我们以7个数据为例，虽然在实际应用中少量数据不适合用核密度的方法进行可视化：
示例数据为：example_data = [-1.95, -1.5, -0.7, -0.65, -0.62, 0.1, 0.9]，可以画出当窗口，宽度分别为0.2，0.4 和 0.8 时的效果如下所示：
可以看到的是，当窗口长度的选择较小时数据的一些细节就会被放大如窗口长度为 0.2 时，整个数据呈现一种类似于迪拉克分布的效果。当窗口长度选择较宽时，数据的一些细节就会被平滑过度，如窗口长度为 0.8 时数据的一些次峰分布就被覆盖。
下面我们用两个分布为 N(−1,0.65)N(-1, 0.65)N(−1,0.65) 和 N(1,0.65)N(1, 0.65)N(1,0.65) 累加的分布产生 200 个随机数据，很明显这两个正态分布累加起来的分布是一个双峰分布如下所示：
然而当我们采用随机数据，并根据不同的窗口长度进行核密度画图时，可视化效果分别如下：
从上述图像中，我们也验证了前面的结论及窗口长度选择过小时，一些细节可能会被放大，而窗口长度选择过大时，一些细节可能会被覆盖，而无法可见。
下面我们将用间歇泉爆发持续时间的数据，采用核密度的方法进行画图，可视化效果如下所示：
使用核密度方法画图的一个缺点，除了上述窗口长度的选择之外，对某些单向拖尾的数据，也可能会因为核密度的选择而造成某些细节表达上的错误。
如我们的精神疾病治疗时长的数据中，数据的取值必须是大于0的，而核密度函数却规定数据的取值范围为 [−∞,∞][-\infty, \infty][−∞,∞]。我们以示例数据为例，画出窗口长度为 20 和60 时的可视化效果：
我们可以从图一中看出，当窗口长度选择较小时会导致右边的拖尾部分出现一定的噪声，创口长度选择较大时，虽然能够平滑过渡这些噪声，但却导致主体部分的细节无法完全展示。
Python代码实现
# -*- coding: utf-8 -*-
"""
Created on Sun Sep 19 22:05:00 2021@author: zhuozb根据ISO13528的参考文件，三十中的2.4节和密度估计画图
"""import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import norm
from histogram_example_showing import extract_data
from naive_estimator_example_showing import digits_after_decimal
from scipy.stats import binom
import randomdef kernel_function(x):'''和密度函数'''return norm.pdf(x)def plot_kernel_estimate(data, fn, window_width=0.4, size=None,title='', interval=None):# 样本容量n = len(data)# 返回数列中最小的数，min_num = min(data)# 求出数列中的最大数，以此来产生画图范围max_num = max(data)if interval == None:# 若没有给出画图间隔， 则根据数据中的最小值的小数点后的位数自动生成# 找出最小数值的小数点后的位数digits_decimal = digits_after_decimal(min_num)# 产生画图范围interval = 0.1**(digits_decimal)if size == None:# 产生画图和画图范围xtick_array = np.arange(0.9*min_num if (min_num>0) else 1.1*min_num, 1.1*max_num, interval)else:xtick_array = np.arange(size[0], size[1], interval)ytick_array = []for x in xtick_array:x_trans = (x - data)/window_widthy_norm = fn(x_trans)y = np.sum(y_norm)/(n*window_width)ytick_array.append(y)plt.plot(xtick_array, ytick_array, linewidth=1)plt.ylim((0, 1.1*max(ytick_array)))# 设置X轴和Y轴的标签plt.xlabel(title+f'h={window_width}', fontfamily='simhei')plt.ylabel('概率密度', fontfamily='simhei') # 设置网格线plt.grid(axis='both')def example_showing(window_width=0.4, title=''):'''该函数主要用于画出实例图'''example_data = [-1.95, -1.5, -0.7, -0.65, -0.62, 0.1, 0.9]plot_kernel_estimate(example_data, kernel_function, title=title, window_width=window_width, size=(-4, 3))# 样本容量n = len(example_data)# 返回数列中最小的数，min_num = min(example_data)# 求出数列中的最大数，以此来产生画图范围max_num = max(example_data)size = (-4, 3)if size == None:# 产生画图和画图范围xtick_array = np.arange(0.9*min_num if (min_num>0) else 1.1*min_num, 1.1*max_num, 0.1)else:xtick_array = np.arange(size[0], size[1], 0.1)def plot_individual(x):# 画出原始数据，每一个点的正态分布情况# 求出X轴对应Y的已被画图之用ytick_array = (xtick_array - x)/window_widthytick_array = norm.pdf(ytick_array)/(n*window_width)# 画出原始数据的正态分布plt.plot(xtick_array, ytick_array)# 标出原始数据点，并画出一条平行于Y轴的直线plt.plot(x, 0, marker='x', color='black', markersize=10)plt.vlines(x, 0, norm.pdf(0)/(n*window_width), color='black')for x in example_data:# 画出每一个原始数据的正态分布图plot_individual(x)def bimodal_example_showing(data, kernel_function, window_width, interval=0.01, title=''):plot_kernel_estimate(data, kernel_function, window_width,size=(-3, 3), title=title, interval=0.01)if __name__ == '__main__':# 画出参考文献中的图，2.42.5for window_width in (0.2, 0.4, 0.8):fig = plt.figure()example_showing(window_width, title='简单案例')plt.show()# 产生一个双峰分布的数据及数据量为200个# 数据的峰值分别为-1和1数据使用正态分布，产生正态分布的方差为0.65np.random.seed(seed=3)def generate_bimodal(low1=-2.5, high1=0.5, mode1=-1,low2=-0.5, high2=2.5, mode2=1):data = [norm.rvs(-1, 0.65) if random.choice((0,1))else norm.rvs(1, 0.65) for i in range(200)]return datadata = generate_bimodal()# 画出参考文献中的图2.6for window_width in (0.2, 0.4, 0.8):fig = plt.figure()bimodal_example_showing(data, kernel_function, window_width, title='双峰分布图')# 画出双峰分布的原始概率密度分布图def plot_bimodal_kernal():data = [norm.rvs(-1, 0.65) if random.choice((0,1))else norm.rvs(1, 0.65) for x in range(20000)]fig = plt.figure()plot_kernel_estimate(data, kernel_function, window_width=0.4, size=(-3,3),interval=0.1)plt.ylabel('概率密度', fontfamily='simhei')plt.xlabel('双峰分布原图', fontfamily='simhei')plt.show()plot_bimodal_kernal()# 读取两个表格的示例数据，并画出它们的和密度图# 读取间歇泉爆发持续时间数据表eruption_lenght_path = r'../数据/间歇泉爆发持续时间表.xlsx'eruption_lenght_data = extract_data(eruption_lenght_path)# 画出图2.8fig = plt.figure()plot_kernel_estimate(eruption_lenght_data, kernel_function, window_width=0.25,size=(0,6), title='间歇泉爆发持续时间 ')plt.show()# 读取自杀治疗时间时长数据treatment_spell_path = r'../数据/自杀倾向治疗时间.xlsx'treatment_spell_data = extract_data(treatment_spell_path)# 画出图2.9a和2.9bfor h in (20, 60):fig = plt.figure()plot_kernel_estimate(treatment_spell_data, kernel_function, window_width=h,size=(-200, 1000), title='间歇泉爆发持续时间 ')plt.show()
最小近邻方法的核密度估计
原理探讨与分析
从上述方法的讨论与分析中，可以看到固定的窗口宽度或者说带宽宽度，会导致可视化效果不是有噪声，就是会把主体部分的细节覆盖掉。因此考虑能否根据数据与数据之间的距离，来动态的决定窗口长度呢？最小近邻方法便解决了这个问题：
定义当前数据 xxx 与原始数据集 Xi,i∈(1,2,⋯,n)X_i, i\in(1,2,\cdots, n)Xi,i∈(1,2,⋯,n)的距离为 di=∣x−Xi∣d_i = |x- X_i|di=∣x−Xi∣，将数据按升序进行排序，并找到第 k 个最小距离 dkd_kdk。其中K为近邻个数是一个大于0的整数，从而可得经验密度函数如下：
f^(x)=k2ndk\hat{f}(x) = \frac{k}{2n d_k} f^(x)=2ndkk
这个公式其实很好理解，假设数据点的概率密度为f(x)f(x)f(x)，那么当带宽 rrr 较小时，我们期望有 2rnf(x)2 r n f(x)2rnf(x) 个数据落在带宽 [x−r,x+r][x-r, x+r][x−r,x+r]。过来说，若有 k 个数据落在区间 [x−dk,x+dk][x-d_k, x+d_k][x−dk,x+dk] 中，即 k=2dknf(x)k = 2d_k n f(x)k=2dknf(x)则可以估计概率密度函数为：f^(x)=k2ndk\hat{f}(x) = \frac{k}{2n d_k}f^(x)=2ndkk。
对于 x<Xminx < X_{min}x<Xmin 的点，由于需要保证我们的窗口宽度不能太大，故需要将 dk=X(k)−xd_k = X_{(k)} - xdk=X(k)−x；同样的，对于 x>Xmaxx > X_{max}x>Xmax，需要将 dk=x−X(n−k+1)d_k = x - X_{(n-k+1)}dk=x−X(n−k+1)，其中 X(i),i∈(1,2,⋯,n)X_{(i)}, i\in(1,2,\cdots,n)X(i),i∈(1,2,⋯,n) 为原始数据 Xi,i∈(1,2,⋯,n)X_i, i\in(1,2,\cdots,n)Xi,i∈(1,2,⋯,n) 的次序统计量。
使用最小近邻方法进行密度估计的一个缺点是，最小近邻方法的核密度函数，不是一个概率密度函数，它对于X轴上的积分面积不是单位面积。因此在使用此种方法进行可视化时，会导致数据的拖尾非常严重。
此外使用该方法进行密度估计时，在核密度 dkd_kdk 的中断点上，核密度是不可导的，这就会导致可视化的效果具有较多的毛刺，我们已见星权爆发持续时间的示例数据为例，采用最小近邻方法进行密度估计可视化效果如下所示：
为了缓解这个问题，我们可以采用正态分布的核密度函数与最小近邻方法相结合，其经验密度函数的估计公式如下：
f^(x)=1ndk∑i=1nK(x−Xidk)\hat{f}(x)=\frac{1}{n d_{k}} \sum_{i=1}^{n} K\left(\frac{x-X_{i}}{d_{k}}\right) f^(x)=ndk1i=1∑nK(dkx−Xi)
用示例数据的可视化效果如下，可以看到虽然所在的问题有所缓和，但依旧没有解决到本质上的问题，也即：
- 核密度函数不是概率密度函数
- 在某些点不可导
Python 代码实现
# -*- coding: utf-8 -*-
"""
Created on Mon Sep 20 09:12:34 2021@author: zhuozb本代码主要展示了依据标准，ISO13528的参考文献，30的2.5小节中的最小邻近值法画出合密度图
"""import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from histogram_example_showing import extract_data
from naive_estimator_example_showing import digits_after_decimal
from kernel_estimate_example_showing import kernel_functiondef nearest_neighbour_method(data, k=20, fn='uniform', title='', size=None,interval=None):# 样本容量n = len(data)# 返回数列中最小的数，min_num = min(data)# 求出数列中的最大数，以此来产生画图范围max_num = max(data)if interval == None:# 若没有给出画图间隔， 则根据数据中的最小值的小数点后的位数自动生成# 找出最小数值的小数点后的位数digits_decimal = digits_after_decimal(min_num)# 产生画图范围interval = 0.1**(digits_decimal)if size == None:# 产生画图和画图范围xtick_array = np.arange(0.9*min_num if (min_num>0) else 1.1*min_num, 1.1*max_num, interval)else:xtick_array = np.arange(size[0], size[1], interval)ytick_array = []# 对原始数据进行排序data_sort = np.sort(data)data_k = data_sort[0]data_n = data_sort[-1]for x in xtick_array:# 若当前的X点小于原始数据的一K个最小值或者大于原始数据的最大值，则 d 用另外一种方式去求解if x < data_k:diff_k = data_sort[k-1] - xelif x > data_n:diff_k = x - data_sort[n-k]else:# 求出当前的X与数据点的距离diff = np.abs(x - data)# 找出第K个最小距离diff_sort = np.sort(diff)diff_k = diff_sort[k-1]if fn == 'uniform':# 工具是2.3，求出Fy = k/(2*n*diff_k)else:y_fn = fn((x-data)/diff_k)y = np.sum(y_fn)/(n*diff_k)ytick_array.append(y)plt.plot(xtick_array, ytick_array, linewidth=1)plt.ylim((0, 1.1*max(ytick_array)))# 设置X轴和Y轴的标签plt.xlabel(title+f'k={k}', fontfamily='simhei')plt.ylabel('概率密度', fontfamily='simhei') # 设置网格线plt.grid(axis='both')if __name__ == '__main__':# 读取间歇泉爆发持续时间数据表eruption_lenght_path = r'../数据/间歇泉爆发持续时间表.xlsx'eruption_lenght_data = extract_data(eruption_lenght_path)# 画出图2.10nearest_neighbour_method(eruption_lenght_data, k=20,size=(0,6), title='间歇泉爆发持续时间 ')# 结合核密度方法和最小近邻方法nearest_neighbour_method(eruption_lenght_data, k=20, fn=kernel_function, size=(0,6), title='间歇泉爆发持续时间 ')
动态核密度方法
原理探讨与分析
受到最小近邻方法的启发，可以结合原始数据中数据与数据之间的间距，来动态的调整窗口的宽度，设 djkd_jkdjk 为数据集 XiX_iXi 与剩余数据 XjX_jXj 的第 [k/2][k/2][k/2] 个最小距离，从而有经验密度函数如下所示：
f^(t)=1n∑j=1n1hdj,kK(t−Xjhdj,k)\hat{f}(t)=\frac{1}{n} \sum_{j=1}^{n} \frac{1}{h d_{j, k}} K\left(\frac{t-X_{j}}{h d_{j, k}}\right) f^(t)=n1j=1∑nhdj,k1K(hdj,kt−Xj)
采用此种方法需要选择最小近邻个数 k 以及窗口长度（真心觉得解决了一个问题，又产生了另外一个问题），以精神治疗时长的数据为例我们采用动态和密度方法设置 k=8,h=5k=8, h=5k=8,h=5，可视化效果如下所示：
从上图可以看出采用动态核密度方法能够将右边的拖尾部分的噪声平滑过渡掉，并且能够体现出主体部分的具体细节。并且相较于最小净灵方法的核密度，估计来说也解决了核密度为概率密度以及可导的问题，因此可谓是一种好方法。但一个缺点就在于需要选择两个参数，具有较大的主观性。
Python 代码实现
# -*- coding: utf-8 -*-
"""
Created on Mon Sep 20 10:02:00 2021@author: zhuozb本代码实现了依据标准ISO13528里面的参考文献30里的公示，2.6展示的可变的核密度估计
"""import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from histogram_example_showing import extract_data
from naive_estimator_example_showing import digits_after_decimal
from kernel_estimate_example_showing import kernel_functiondef variable_kernel_estimate(data, fn=kernel_function, window_width=5, k=8, title='',size=None, interval=None):# 样本容量n = len(data)# 返回数列中最小的数，min_num = min(data)# 求出数列中的最大数，以此来产生画图范围max_num = max(data)if interval == None:# 若没有给出画图间隔， 则根据数据中的最小值的小数点后的位数自动生成# 找出最小数值的小数点后的位数digits_decimal = digits_after_decimal(min_num)# 产生画图范围interval = 0.1**(digits_decimal)if size == None:# 产生画图和画图范围xtick_array = np.arange(0.9*min_num if (min_num>0) else 1.1*min_num, 1.1*max_num, interval)else:xtick_array = np.arange(size[0], size[1], interval)ytick_array = []for x in xtick_array:# 最美一个坐标轴X上的数据根据公式2.6求出 Y轴的数据y_sum = 0for x_data in data:# 对原始数据的每一个点，找出其与其他 n- 1个点的第 k 个各最小近邻距离diff = np.abs(data-x_data)diff_sort = np.sort(diff)diff_jk = diff_sort[round(k/2)]# 根据公式2.6求出密度估计值x_trans = (x-x_data)/(window_width*diff_jk)y_fn = fn(x_trans)/(window_width*diff_jk)y_sum = y_sum + y_fny = y_sum/nytick_array.append(y)plt.plot(xtick_array, ytick_array, linewidth=1)plt.ylim((0, 1.1*max(ytick_array)))# 设置X轴和Y轴的标签plt.xlabel(title+f'window_width={window_width}, k={k}', fontfamily='simhei')plt.ylabel('概率密度', fontfamily='simhei') # 设置网格线plt.grid(axis='both')if __name__ == '__main__':# 读取自杀治疗时间时长数据treatment_spell_path = r'../数据/自杀倾向治疗时间.xlsx'treatment_spell_data = extract_data(treatment_spell_path)# 画出图2.11variable_kernel_estimate(treatment_spell_data, window_width=5, k=8, fn=kernel_function, size=(-200,1000), title='自杀倾向治疗时间 ',interval=1)
数据的核密度估计及其可视化：Python实现相关推荐
- 毕业设计 大数据B站数据分析与可视化 - python 数据分析 大数据
文章目录 0 前言 1 课题背景 2 实现效果 3 数据获取 4 数据可视化 5 最后 0 前言
- 【毕业设计】大数据北京二手房数据分析与可视化 - python 数据挖掘
文章目录 0 前言 1 探索性分析与文本数据预处理 2 数据可视化分析 2.1 Region特征分析 2.2 Year 特征分析 2.3 Floor 特征分析 3 最后 0 前言
- 大数据毕业设计 住房数据分析与可视化 - python
文章目录 0 前言 分析展示 一.北上广租房房源分布可视化 二.北上广内区域租金分布可视化 三.房源距地铁口租金的关系可视化 四.房屋大小与租金关系可视化 结论 租个人房源好还是公寓好 北上广深租房时 ...
- 数据分析毕业设计 大数据商城人流数据分析与可视化 - python
文章目录 0 前言 课题背景 分析方法与过程 初步分析: 总体流程: 1.数据探索分析 2.数据预处理 3.构建模型 总结 0 前言
- 【毕业设计】 大数据二手房数据爬取与分析可视化 -python 数据分析 可视化
1 前言
- 数据分析毕业设计 大数据京东消费行为分析与可视化 - python 机器学习
1 前言
- 【毕业设计】 大数据上海租房数据爬取与分析可视化 -python 数据分析 可视化
1 前言
- 大数据毕业设计 招聘网站数据分析可视化 - python flask 网络爬虫
文章目录 0 前言 1 课题背景 2 实现效果 3 Flask框架 4 Echarts 5 爬虫 6 最后 0 前言
- 核密度估计及其Python实践
一.参数估计简介 很多情况下,我们只有有限的样本集,而类条件概率密度函数p(x|ωi)和先验概率P(ωi)是未知的,需要根据已有样本进行参数估计,然后将估计值当作真实值来使用. 由给定样本集求解随机变 ...
最新文章
- ES6新特性之Promise
- 一个好用的小工具 thefuck
- python telnetlib详解 执行循环命令_Python3+telnetlib实现telnet客户端
- C#处理控制台关闭事件
- apache SSL配置
- “智慧城市”如火如荼 与“数字城市”又有何差别？
- WOL局域网与外网远程唤醒概要
- iOS 银行卡号识别
- 十分好用的拓扑图插件JTopo
- Cdn英文的读音音标_宋sir的美式音标教程 Unit 1 /i/ tea
- 很多网友问那个磁力搜索站好用，就由本君说说吧！
- 解释下ArrayList集合为啥允许值为null
- Ubuntu 使用firefox插件下载百度云文件
- 算法专题训练（1）股票问题
- vi,vim文本编辑器
- CSUSTOJ-藤原千花不想知道数学成绩（数组及无数组解法）
- 关闭localized intellij idea切换语言提醒
- 计算机毕设结束语致谢,毕业设计结束语和致谢
- linux打开文件乱码
- 扫雷 洛谷p2327