精华内容
下载资源
问答
  • 文章目录实例数据介绍精神疾病治疗时长数据表间歇泉爆发持续时间数据表文件放置介绍频率分布直方图法原理探讨与分析Python 代码实现简单核密度估计法原理探讨与分析Python代码实现核密度估计法原理探讨与分析Python...


    本文的参考文献主要如下:

    1. Silverman B.W. Density Estimation. Chapman and Hall, London, 1986
    2. 标准 ISO 13528

    实例数据介绍

    本博客采用的示例数据为本博客的参考文献一中的表 2.1 和 2.2,如下所示:

    精神疾病治疗时长数据表

    精神疾病治疗时长数据表

    间歇泉爆发持续时间数据表

    间歇泉爆发持续时间数据表

    文件放置介绍

    我们将上述的数据表放置在代码文件夹,所在路径的另一个命名为 “数据” 的文件夹中:

    代码:
    -------> xxx.py (代码)
    数据:
    -------> 间歇泉爆发持续时间表
    -------> 自杀倾向治疗时间

    频率分布直方图法

    原理探讨与分析

    设数据集为 X = ( X 1 , X 2 , ⋯   , X n ) X=(X_1, X_2, \cdots, X_n) X=(X1,X2,,Xn)

    频率分布直方图需要我们提供画图原点 x 0 x_0 x0 以及一个带宽 h h h,将落在 [ x 0 + ( m − 1 ) ⋅ h , x 0 + m ⋅ h ] [x_0 + (m-1)\cdot h, x_0 + m \cdot h] [x0+(m1)h,x0+mh] 上的点的数量进行计数,后除以带宽长度即可:
    f ^ ( x ) = 1 n h ( 落在与 x 同一带宽上的 X i 的数量 ) \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 pd
    
    def extract_data(path):
        # 将Excel数据表读入Python中
        table = pd.read_excel(path, header=None)
        # 将数据表转换为 2D 数列
        data = table.values
        # 剔除非数字数据并将一维数列转换为二维数列
        data = data[~np.isnan(data)]
        return data
    
    
    
    
    def 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 bar
                       xytext = (0,0.2),             # offsetting label position above its bar
                       textcoords = "offset points", # Offset (in points) from the *xy* value
                       ha = '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.1
        treatment_spell_data = extract_data(treatment_spell_path)
        # 画出表格2.2
        eruption_lenght_data = extract_data(eruption_lenght_path)
        
        # 画出图2.1
        plot_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.2
        plot_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 ) = lim ⁡ h → 0 1 2 h P ( x − h < X < x + h ) f(x) = \lim_{h\to0} \frac{1}{2h} P(x-h < X < x+h) f(x)=h0lim2h1P(xh<X<x+h)
    于是我们可以用带宽为无限小的频数分布图来估计概率密度:
    f ^ ( x ) = 1 2 n h ( 落在带宽 ( x − h , x + h ) 上的 X i 的数量 ) \hat{f}(x) = \frac{1}{2nh} (\text{落在带宽} (x-h, x+h) \text{上的} X_i \text{的数量}) f^(x)=2nh1(落在带宽(xh,x+h)上的Xi的数量)
    为了方便表述我们定义一个核密度函数如下所示:
    w ( x ) = { 1 2  if  ∣ x ∣ < 1 0  otherwise.  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 ) = 1 h n ∑ i = 1 n w ( x − X i h ) \hat{f}(x) = \frac{1}{hn} \sum_{i=1}^{n} w ( \frac{x-X_i}{h} ) f^(x)=hn1i=1nw(hxXi)
    在实际应用中带宽当然可以不需要取到无限接近于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_data
    
    
    def digits_after_decimal(num):
        '''
        返回变量 num 的小数点位数
        '''
        num = str(num)
        # 小数点所在位置
        digit_loc = num.find('.')
    
        # 小数点后的位数
        digits_decimal = len(num) - digit_loc - 1
        
        return digits_decimal
    
    def 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对应的Y
            y_sum = 0
            for 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)])/2
                y_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.1
        treatment_spell_data = extract_data(treatment_spell_path)
        # 画出表格2.2
        eruption_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 ) = 1 2 π exp ⁡ ( − x 2 2 ) 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 ) = 1 n h ∑ i = 1 n K ( x − X i h ) \hat{f}(x) = \frac{1}{nh} \sum_{i=1}^{n} K ( \frac{x-X_i}{h} ) f^(x)=nh1i=1nK(hxXi)
    在核密度估计中,一般把参数 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 random
    
    def 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_width
            y_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_width
            ytick_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.5
        for window_width in (0.2, 0.4, 0.8):
            fig = plt.figure()
            example_showing(window_width, title='简单案例')
            plt.show()
        
        # 产生一个双峰分布的数据及数据量为200个
        # 数据的峰值分别为-1和1数据使用正态分布,产生正态分布的方差为0.65
        np.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 data
        
        data = generate_bimodal()
        
        # 画出参考文献中的图2.6
        for 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.8
        fig = 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.9b
        for h in (20, 60):
            fig = plt.figure()
            plot_kernel_estimate(treatment_spell_data, kernel_function, window_width=h,
                                 size=(-200, 1000), title='间歇泉爆发持续时间 ')
            plt.show()
    

    最小近邻方法的核密度估计

    原理探讨与分析

    从上述方法的讨论与分析中,可以看到固定的窗口宽度或者说带宽宽度,会导致可视化效果不是有噪声,就是会把主体部分的细节覆盖掉。因此考虑能否根据数据与数据之间的距离,来动态的决定窗口长度呢?最小近邻方法便解决了这个问题:

    定义当前数据 x x x 与原始数据集 X i , i ∈ ( 1 , 2 , ⋯   , n ) X_i, i\in(1,2,\cdots, n) Xi,i(1,2,,n)的距离为 d i = ∣ x − X i ∣ d_i = |x- X_i| di=xXi,将数据按升序进行排序,并找到第 k 个最小距离 d k d_k dk。其中K为近邻个数是一个大于0的整数,从而可得经验密度函数如下:
    f ^ ( x ) = k 2 n d k \hat{f}(x) = \frac{k}{2n d_k} f^(x)=2ndkk
    这个公式其实很好理解,假设数据点的概率密度为 f ( x ) f(x) f(x),那么当带宽 r r r 较小时,我们期望有 2 r n f ( x ) 2 r n f(x) 2rnf(x) 个数据落在带宽 [ x − r , x + r ] [x-r, x+r] [xr,x+r]。过来说,若有 k 个数据落在区间 [ x − d k , x + d k ] [x-d_k, x+d_k] [xdk,x+dk] 中,即 k = 2 d k n f ( x ) k = 2d_k n f(x) k=2dknf(x)则可以估计概率密度函数为: f ^ ( x ) = k 2 n d k \hat{f}(x) = \frac{k}{2n d_k} f^(x)=2ndkk

    对于 x < X m i n x < X_{min} x<Xmin 的点,由于需要保证我们的窗口宽度不能太大,故需要将 d k = X ( k ) − x d_k = X_{(k)} - x dk=X(k)x;同样的,对于 x > X m a x x > X_{max} x>Xmax,需要将 d k = x − X ( n − k + 1 ) d_k = x - X_{(n-k+1)} dk=xX(nk+1),其中 X ( i ) , i ∈ ( 1 , 2 , ⋯   , n ) X_{(i)}, i\in(1,2,\cdots,n) X(i),i(1,2,,n) 为原始数据 X i , i ∈ ( 1 , 2 , ⋯   , n ) X_i, i\in(1,2,\cdots,n) Xi,i(1,2,,n) 的次序统计量。

    使用最小近邻方法进行密度估计的一个缺点是,最小近邻方法的核密度函数,不是一个概率密度函数,它对于X轴上的积分面积不是单位面积。因此在使用此种方法进行可视化时,会导致数据的拖尾非常严重。

    此外使用该方法进行密度估计时,在核密度 d k d_k dk 的中断点上,核密度是不可导的,这就会导致可视化的效果具有较多的毛刺,我们已见星权爆发持续时间的示例数据为例,采用最小近邻方法进行密度估计可视化效果如下所示:
    在这里插入图片描述

    为了缓解这个问题,我们可以采用正态分布的核密度函数与最小近邻方法相结合,其经验密度函数的估计公式如下:

    f ^ ( x ) = 1 n d k ∑ i = 1 n K ( x − X i d k ) \hat{f}(x)=\frac{1}{n d_{k}} \sum_{i=1}^{n} K\left(\frac{x-X_{i}}{d_{k}}\right) f^(x)=ndk1i=1nK(dkxXi)

    用示例数据的可视化效果如下,可以看到虽然所在的问题有所缓和,但依旧没有解决到本质上的问题,也即:

    1. 核密度函数不是概率密度函数
    2. 在某些点不可导

    在这里插入图片描述

    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_function
    
    
    def 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] - x
            elif 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,求出F
                y = 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.10
        
        nearest_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='间歇泉爆发持续时间 ')
    

    动态核密度方法

    原理探讨与分析

    受到最小近邻方法的启发,可以结合原始数据中数据与数据之间的间距,来动态的调整窗口的宽度,设 d j k d_jk djk 为数据集 X i X_i Xi 与剩余数据 X j X_j Xj 的第 [ k / 2 ] [k/2] [k/2] 个最小距离,从而有经验密度函数如下所示:
    f ^ ( t ) = 1 n ∑ j = 1 n 1 h d j , k K ( t − X j h d j , 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=1nhdj,k1K(hdj,ktXj)
    采用此种方法需要选择最小近邻个数 k 以及窗口长度(真心觉得解决了一个问题,又产生了另外一个问题),以精神治疗时长的数据为例我们采用动态和密度方法设置 k = 8 , h = 5 k=8, h=5 k=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_function
    
    def 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 = 0
            for 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_fn
                
            y = y_sum/n
            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'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.11
        variable_kernel_estimate(treatment_spell_data, window_width=5, 
                                 k=8, fn=kernel_function, 
                                 size=(-200,1000), title='自杀倾向治疗时间 ',
                                 interval=1)
    
    展开全文
  • 在后面的文章中,我们将采用参数化的广义线性混合模型,并展示如何切换到非参数化的随机效应表示,避免了正态分布的随机效应假设。 使用Dirichlet Process Mixture模型进行基本密度估计 提供了通过Dirichlet过程...

    原文链接:http://tecdat.cn/?p=23785

    原文出处:拓端数据部落公众号

    概述

    最近,我们使用贝叶斯非参数(BNP)混合模型进行马尔科夫链蒙特卡洛(MCMC)推断。

    在这篇文章中,我们通过展示如何使用具有不同内核的非参数混合模型进行密度估计。在后面的文章中,我们将采用参数化的广义线性混合模型,并展示如何切换到非参数化的随机效应表示,避免了正态分布的随机效应假设。

    使用Dirichlet Process Mixture模型进行基本密度估计

    提供了通过Dirichlet过程混合(DPM)模型进行非参数密度估计的机制(Ferguson, 1974; Lo, 1984; Escobar, 1994; Escobar and West, 1995)。对于一个独立和相同分布的样本 ,该模型的形式为

    这个模型实现是灵活的,运行任意核的混合。, 可以是共轭的,也可以是不共轭的(也是任意的)基度量 . 在共轭核/基数测量对的情况下,能够检测共轭的存在,并利用它来提高采样器的性能。

    为了说明这些能力,我们考虑对R中提供的Faithful火山数据集的喷发间隔时间的概率密度函数进行估计。

    data(faithful)

    观测值  对应于数据框架的第二列,而 .

    使用CRP表示法拟合高斯location-scale 分布混合分布

    模型说明

    我们首先考虑用混合正态分布的location-scaleDirichlet过程s来拟合转换后的数据

    其中 对应的是正态-逆伽马分布。这个模型可以解释为提供一个贝叶斯版本的核密度估计 用于使用高斯核和自适应带宽。在数据的原始尺度上,这可以转化为一个自适应的对数高斯核密度估计。

    引入辅助变量,表明混合的哪个成分产生了每个观测值,并对随机量进行积分,我们得到模型的CRP表示(Blackwell and MacQueen, 1973)。

    其中

    是向量中唯一值的数量,是第个唯一值在中出现的次数。这个说明清楚地表明,每个观测值都属于最多正态分布聚类中的任何一个,并且CRP分布与分区结构的先验分布相对应。

    这个模型的说明是这样的

        y[i] ~ dnorm(mu[i], var = s2[i])
        mu[i] <- muTilde[xi[i]]
        s2[i] <- s2Tilde[xi[i]]
      xi[1:n] ~ dCRP(alpha, size = n)
        muTilde[i] ~ dnorm(0, var = s2Tilde[i])
        s2Tilde[i] ~ dinvgamma(2, 1)
      alpha ~ dgamma(1, 1)
    

    请注意,在模型代码中,参数向量muTilde和s2Tilde的长度被设置为.我们这样做是因为目前的实现要求提前设置参数向量的长度,并且不允许它们的数量在迭代之间变化。因此,如果我们要确保算法总是按预期执行,我们需要在最坏的情况下工作,即有多少个成分就有多少个观测值的情况。但它的效率也有点低,无论是在内存需求方面(当 规模大时,需要维护大量未占用的成分)还是在计算负担方面(每次迭代都需要更新大量不需要后验推理的参数)。当我们在下面使用伽马分布的混合时,我们将展示一个能提高效率的计算捷径。

    还需要注意的是,的值控制着我们先验预期的成分数量,的值越大,对应于数据占据的成分数量越多。因此,通过指定一个先验值,我们为模型增加了灵活性。对Gamma先验的特殊选择允许使用数据增强方案从相应的全条件分布中有效取样。也可以选择其他的先验,在这种情况下,这个参数的默认采样是一个自适应的随机游走Metropolis-Hastings算法。

    运行MCMC算法

    下面的代码设置了数据和常数,初始化了参数,定义了模型对象,并建立和运行了MCMC算法。默认采样器是一个折叠的吉布斯采样器(Neal, 2000)。

    
    # 模型数据
    y = standlFaithful
    # 模型常量
    n = length(standlFaithful))
    # 参数初始化
    list(xi = sample(1:10, size=n, replace=TRUE),
    # 创建和编译模型
    Model(code,  data,  inits,  consts)
    ##定义模型...
    ##建立模型...
    ##设置数据和初始值...
    ##在模型上运行计算(随后的任何错误报告可能只是反映了模型变量的缺失值)... 
    ##检查模型的大小和尺寸......。
    ##模型构建完成。
    ## 编译完成。
    #MCMC的配置、创建和编译
    MCMC(conf)
    ## 编译......这可能需要一分钟
    ## 编译完成。
    
    

     

    我们可以从参数的后验分布中提取样本,并创建痕迹图、直方图和任何其他感兴趣的总结。例如,对于参数,我们有。
     

    # 参数的痕迹图
    ts.plot(samples[ , "alpha"], xlab = "iteration", ylab = expression(alpha))

    
    # 后验直方图
    hist(samples[ , "alpha"], xlab = expression(alpha), main = "", ylab = "Frequency")
    
    

    在这个模型下,对于一个新的观察,后验预测分布是最佳密度估计(在平方误差损失下)。这个估计的样本可以很容易地从我们的MCMC产生的样本中计算出来。

    
    # 参数的后验样本
     samples[, "alpha"]
    # 平均值的后验样本
     samples[, grep('muTilde', colnames(samples))] # 聚类平均数的后验样本。
    # 集群方差的后验样本
    samples[, grep('s2Tilde', colnames(samples))] # 聚类成员的后验样本。
    # 聚类成员关系的后验样本
    samples [, grep('xi', colnames(samples))] # 聚类成员的后验样本。
    
    hist(y, freq = FALSE,
         xlab = "标准化对数尺度上的等待时间")
    ##对标准化对数网格的密度进行点式估计
    
    

    然而,回顾一下,这是对等待时间的对数的密度估计。为了获得原始尺度上的密度,我们需要对内核进行适当的转换。 

    standlGrid*sd(lFaithful) + mean(lFaithful) # 对数尺度上的网格
    
    hist(faithful$waiting, freq = FALSE
      
    

    无论是哪种情况,都有明显的证据表明,数据中的等待时间有两个组成部分。

    生成混合分布的样本

    虽然混合分布的线性函数的后验分布的样本(比如上面的预测分布)可以直接从折叠采样器的实现中计算出来,但是对于非线性函数的推断需要我们首先从混合分布中生成样本。我们可以从随机度量中获得后验样本。需要注意的是,为了从 ,得到后验样本,我们需要监控所有参与其计算的随机变量,即成员变量xi,聚类参数muTilde和s2Tilde,以及浓度参数alpha。

    下面的代码从随机测量中生成后验样本。cMCMC对象包括模型和参数的后验样本。函数估计了一个截断水平,即truncG。后验样本是一个带列的矩阵,其中参数分布向量的维度(在本例中为)。

    outputG <- getSamplesDPmeasure(cmcmc)
     
    

    下面的代码使用随机测量的后验样本来计算的后验样本。请注意,这些样本是基于转换后的模型计算的,大于70的值对应于上述定义的网格上大于0.035的值。

    
         truncG <- outputG$trunc # G的截断水平
    
    
    
    probY70 <- rep(0, nrow(samples))  # P(y.tilde>70)的后验样本
    
    hist(probY70 )
    
    

    使用CRP表示法拟合伽马混合分布

    不限于在DPM模型中使用高斯核。就Old Faithful数据而言,除了我们在上一节中介绍的对数尺度上的高斯核的混合分布外,还有一种选择是数据原始尺度上的伽马混合分布。

    模型

    在这种情况下,模型的形式为

    其中对应于两个独立Gamma分布的乘积。下面的代码提供了该模型。

        y[i] ~ dgamma(shape = beta[i], scale = lambda[i])
        beta[i] <- betaTilde[xi[i]]
        lambda[i] <- lambdaTilde[xi[i]]
    

    请注意,在这种情况下,向量beta和lambda的长度为 。这样做是为了减少与采样算法有关的计算和存储负担。你可以把这种方法看作是对过程的截断,只不过它可以被认为是*精确的截断。事实上,在CRP表示法下,只要采样器的成分数严格低于采样器每次迭代的参数向量的长度,使用长度短于样本中观察值的参数向量就会生成一个合适的算法。

    运行MCMC算法

    下面的代码设置了模型数据和常数,初始化了参数,定义了模型对象,并建立和运行了Gamma混合分布的MCMC算法。请注意,在构建MCMC时,会产生一个关于聚类参数数量的警告信息。这是因为betaTilde和lambdaTilde的长度小于。另外,请注意,在执行过程中没有产生错误信息,这表明所需的集群数量未超过50个的上限。

    data <- list(y = waiting)
    Model(code, data = data)
     
    

    cModel <- compile
     
    

    samples <- runMCMC(cmcmc, niter = 7000, nburnin = 2000, setSeed = TRUE)
     
    

    在这种情况下,我们使用参数的后验样本来构建一个轨迹图,并估计的后验分布。
     

    # 参数的后验样本的跟踪图
    ts.plot(samples[ , 'alpha'], xlab = "iteration", ylab = expression(alpha))

    
    # 参数的后验样本的直方图 
    hist(samples[, 'alpha'])
    
    

    从混合分布中生成样本

    和以前一样,我们从后验分布中获得样本。
     

    outputG <- getSamplesDPmeasure(cmcmc)
     
    

    我们使用这些样本来创建一个数据密度的估计值,以及一个95%置信带。 

    
    for(iter in seq_len)) {
      density[iter, ] <- sapply(grid, function(x)
        sum( weightSamples[iter, ] * dgamma)))
    }
    
    
    hist(waiting, freq = FALSE
    

    我们再次看到,数据的密度是双峰的,看起来与我们之前得到的数据非常相似。

    使用stick-breaking 表示法拟合伽马DP混合分布

    模型

    Dirichlet过程混合物的另一种表示方法是使用随机分布的stick-breaking表示(Sethuraman, 1994)。


    引入辅助变量,表明哪个成分产生了每个观测值,上一节讨论的Gamma密度的混合物的相应模型的形式为

    其中 是两个独立Gamma分布的乘积。

    . 下面的代码提供了该模型说明。 

    
          y[i] ~ dgamma(shape = beta[i], scale = lambda[i])
          beta[i] <- betaStar[z[i]]
          lambda[i] <- lambdaStar[z[i]]
          z[i] ~ dcat(w[1:Trunc])
         # stick-breaking 
          v[i] ~ dbeta(1, alpha)
        w[1:Trunc] <- stick_breaking(v[1:(Trunc-1)]) # stick-breaking 权重
          betaStar[i] ~ dgamma(shape = 71, scale = 2)
    

    注意,截断水平已被设置为Trunc值,该值将在函数的常数参数中定义。

    运行MCMC算法

    下面的代码设置了模型数据和常量,初始化了参数,定义了模型对象,并建立和运行了Gamma混合分布的MCMC算法。当使用stick-breaking表示时,会指定一个分块Gibbs抽样器(Ishwaran, 2001; Ishwaran and James, 2002)。

    data <- list(y = waiting)
    consts <-length(waiting)
    betaStar = rgamma
                  lambdaStar = rgamma
                  v = rbeta
                  z = sample
                  alpha = 1
    
    

     
    

    compile(Model)
     
    

    MCMC(rModel, c("w", "betaStar", "lambdaStar", 'z', 'alpha'))
    comp(mcmc )
     
    

    MCMC(cmcmc, niter = 24000)
     
    

    使用stick-breaking近似法会自动提供随机分布的近似值,即 。下面的代码使用来自样本对象的后验样本计算后验样本,并从中计算出数据的密度估计。

    
      densitySamples[i, ] <- sapply(grid, function(x) sum(weightSamples  * dgamma(x, shape ,
                                        scale )))
    
    hist( waiting ylim=c(0,0.05),
     

    正如预期的那样,这个估计值看起来与我们通过CRP表示的过程获得的估计值相同。

    贝叶斯非参数化:非参数化随机效应

    我们将采用一个参数化的广义线性混合模型,并展示如何切换到非参数化的随机效应表示,避免了正态分布的随机效应假设。

    心肌梗死(MIs)的参数化meta分析

    我们将在对以前非常流行的糖尿病药物 "Avandia "的副作用进行meta分析的背景下,说明使用非参数混合模型对随机效应分布进行建模。我们分析的数据在引起对这种药物的安全性的严重质疑方面发挥了作用。问题是使用"Avandia "是否会增加心肌梗死(心脏病发作)的风险。,每项研究都有治疗和对照组。

    模型的制定

    我们首先进行基于标准的广义线性混合模型(GLMM)的meta分析。向量n和x分别包含对照组的患者总数和每项研究中对照组的心肌梗死患者人数。同样,向量m和y包含接受药物的病人的类似信息。该模型的形式为


    其中,随机效应、遵循共同的正态分布被合理地赋予非信息性先验。参数量化了对照组和治疗组之间的风险差异,而参数则量化了研究的具体变化。

    这个模型可以用以下代码指定。

    
            y[i] ~ dbin(size = m[i], prob = q[i]) # 药物MIs
            x[i] ~ dbin(size = n[i], prob = p[i]) # 控制MIs
            q[i] <- expit(theta + gamma[i]) # 药物的对数指数
            p[i] <- expit(gamma[i]) #对照组对数
            gamma[i] ~ dnorm(mu, var = tau2) # 研究效果
        theta ~ dflat() # 药物的影响
        # 随机效应超参数
        mu ~ dnorm(0, 10)
        tau2 ~ dinvgamma(2, 1)
    
    

    运行MCMC

    让我们来运行一个基本的MCMC。

    MCMC(codeParam,  data, inits,
                          constants, monitors = c("mu", "tau2", "theta", "gamma")
     
    

    par(mfrow = c(1, 4)
    
    hist(gammaMn)
    hist(samples[1000, gammaCols)

    结果表明,对照组和治疗组之间存在着整体的风险差异。但是正态性假设呢?我们的结论对该假设是否稳健?也许随机效应的分布是偏斜的。

    用于meta分析的基于DP的随机效应模型

    模型

    现在,我们对使用非参数分布。更具体地说,我们假设每个都是由位置尺度的正态混合分布产生的。


    这种模型引起了随机效应之间的聚类。与密度估计问题的情况一样,DP先验允许数据决定分量的数量,从最少的一个分量(即简化为参数模型)到最多的分量,即每个观测值有一个分量。如果数据支持这种行为,这允许随机效应的分布是多模态的,大大增加了其灵活性。这个模型可以用以下代码指定。

    
            y[i] ~ dbin(size = m[i], prob = q[i]) # 药物MIs
            x[i] ~ dbin(size = n[i], prob = p[i]) # MIs
            q[i] <- expit(theta + gamma[i]) # 药物的对数指数
            p[i] <- expit(gamma[i]) # 对照组对数值
            gamma[i] ~ dnorm(mu[i], var = tau2[i]) # 来自混合物的随机效应。
            mu[i]<- muTilde[xi[i]]                 # 来自聚类的随机效应的平均值 xi[i]
            tau2[i] <- tau2Tilde[xi[i]]           # 来自群组xi[i]的随机效应变量
        # 从基础测量中提取混合成分参数
            muTilde[i] ~ dnorm(mu0, var = var0)
            tau2Tilde[i] ~ dinvgamma(a0, b0)
        # 用于将研究报告聚类为混合成分的CRP
        xi[1:nStudies] ~ dCRP(alpha, size = nStudies)
        # 超参数
      
        theta ~ dflat() # 药物的影响
    
    

    运行MCMC

    以下代码对模型进行了编译,并对模型运行了一个压缩Gibbs抽样

    inits <- list(gamma = rnorm(nStudies))
    
    MCMC(code = BNP, data = data)
     
    

    
    hist(samplesBNP[, 'theta'], xlab = expression(theta), main = 'avandia的影响')
                  main = "随机效应分布")
                       main = "随机效应分布")
    
    # 推断出了多少个混合成分?
    xiRes <- samplesBNP[, xiCols].
    
    

    主要推论似乎对原始的参数化假设很稳健。这可能是由于没有太多证据表明随机效应分布中缺乏正态性。

    参考文献

    Blackwell, D. and MacQueen, J. 1973. Ferguson distributions via Polya urn schemes. The Annals of Statistics 1:353-355.

    Ferguson, T.S. 1974. Prior distribution on the spaces of probability measures. Annals of Statistics 2:615-629.

    Lo, A.Y. 1984. On a class of Bayesian nonparametric estimates I: Density estimates. The Annals of Statistics 12:351-357.


    最受欢迎的见解

    1.matlab使用贝叶斯优化的深度学习

    2.matlab贝叶斯隐马尔可夫hmm模型实现

    3.R语言Gibbs抽样的贝叶斯简单线性回归仿真

    4.R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归

    5.R语言中的Stan概率编程MCMC采样的贝叶斯模型

    6.Python用PyMC3实现贝叶斯线性回归模型

    7.R语言使用贝叶斯 层次模型进行空间数据分析

    8.R语言随机搜索变量选择SSVS估计贝叶斯向量自回归(BVAR)模型

    9.matlab贝叶斯隐马尔可夫hmm模型实现

    展开全文
  • 核密度估计的又是什么? KPCA中的的作用是什么? 深度学习CNN中的卷积和过滤有又什么关系? 核密度估计又是什么? 函数和激活函数什么关系? Linux内核叫Kernel? Jupyter Notebook的运行内核也叫...

    一切从好的问题开始

    SVM中为什么需要核函数?

    核密度估计的核又是什么?

    KPCA中的核的作用是什么?

    深度学习CNN中的卷积核?核和过滤有又什么关系?

    核密度估计又是什么?

    核函数和激活函数什么关系?

    Linux内核叫Kernel?

    Jupyter Notebook的运行内核也叫Kernel?

    展开全文
  • 参考url: ...密度评估器是一种利用D维数据集生成D维概率分布估计的算法,GMM算法用不同高斯分布的加权汇总来表示概率分布估计。核密度估计(kernel density estimation,KDE)算法将高斯混合理念扩展...

    参考url:

    https://jakevdp.github.io/PythonDataScienceHandbook/05.13-kernel-density-estimation.html

    密度评估器是一种利用D维数据集生成D维概率分布估计的算法,GMM算法用不同高斯分布的加权汇总来表示概率分布估计。核密度估计(kernel density estimation,KDE)算法将高斯混合理念扩展到了逻辑极限(logical extreme),它通过对每个点生成高斯分布的混合成分,获得本质上是无参数的密度评估器。

    1、KDE的由来:直方图

      密度估计评估器是一种寻找数据集生成概率分布模型的算法。

      一维数据的密度估计——直方图,是一个简单的密度评估器,直方图将数据分成若干区间,统计落入每个区间内的点的数量,然后用直观的方式将结果可视化。

      

      

       

      

       

    2、核密度估计的实际应用

      核密度估计的自由参数是核类型(kernel)参数,他可以指定每个点核密度分布的形状。

      核带宽(kernel bandwidth)参数控制每个点的核的大小

      核密度估计算法在sklearn.neighbors.KernelDensity评估器中实现,借助六个核中的任意一个核、二三十个距离量度就可以处理具有多个维度的KDE。

      由于KDE计算量非常大,因此Scikit-Learn评估器底层使用了一种基于树的算法,可以利用atol(绝对容错)和rtol(相对容错)参数来平衡计算时间与准确性,可以用Scikit-Learn的标准交叉检验工具来确定自由参数核带宽。

      

       通过交叉检验选择带宽

      在KDE中,带宽的选择不仅对找到合适的密度估计非常重要,也是在密度估计中控制偏差-方差平衡的关键:

      (1)带宽过窄将导致估计呈现高方差(即过拟合),而且每个点的出现或缺失都会引起很大的不同

      (2)带宽过宽将导致估计呈现高偏差(即欠拟合),而且带宽较大的核还会破坏数据结构

      机器学习中超参数的调优通常都是通过交叉检验完成的。

      

      

     

    展开全文
  • 核密度估计-KDE

    2021-09-23 23:30:29
    核密度估计的自由参数类型和带宽,前者指定每个点核密度分布的形状,后者指定每个点的大小。 一维数据密度估计——直方图,是一个简单的密度评估器,直方图将数据分成若干区间,统计落入每个区间内的点的...
  • Py之seaborn:数据可视seaborn库(二)的组合图可视密度图/核密度图分布可视、箱型图/散点图、小提琴图/散点图组合可视的简介、使用方法之最强攻略(建议收藏) 目录 二、组合图可视 1、密度图、...
  • 展开全部R语言实际上是函数的集合,用户可以使用base,stats等包中的基本函数,也可以自己编写函数完成一定62616964757a686964616fe4b893e5b19e31333363373835的功能。但是初学者往往认为编写R函数十分困难,或者...
  • 二维高斯分布概率密度函数数据集实战优化坐标轴与图像优化图像再次优化概率密度函数大家肯定都有听说过正态分布,其实正态分布只是概率密度分布的一种,正态分布的概率密度函数均值为μ ,标准差σ是高斯函数的一个...
  • 简介机器学习模型千变万化, 每种模型需要的数据形式都略有不同, 所以我们一般花在数据提取/清理/格式等操作的时间可能会占项目总时间的百分之七十(甚至更多)。不过, 借助别人已经写好的程序可以让我们的工作...
  • 非参数估计,即核密度估计(Kernel Density Estimation,KDE),不需要预先假设,从数据本身出发,来估计未知的密度函数。 一、估计过程 1、以每个点的数据+带宽(邻域)作为参数,用函数估计样本中每个数据点...
  • 在本篇博客中,我将与众位看官分享我自己写的关于用核密度函数加权的直方图的计算 欢饮批评指正!!! 核密度函数加权直方图在基于Mean shift的跟踪算法中经常用到。请客官查看我的博客: 样本均值漂移的原理Mean...
  • 首先介绍seaborn.distplot()函数 seaborn.distplot(a, bins=None, hist=True, kde=True, rug=False, fit=None, hist_kws=None, kde_kws=None, rug_kws=None, fit_kws=None, color=None, ...参数:a:Series、1
  • 二维高斯分布概率密度函数数据集实战优化坐标轴与图像优化图像再次优化概率密度函数大家肯定都有听说过正态分布,其实正态分布只是概率密度分布的一种,正态分布的概率密度函数均值为μ ,标准差σ是高斯函数的一个...
  • 开发人员可以修改表结构,可以随意修改其中的数据但是需要保证不影响其他开发同事。test: 测试环境开发可读写,开发人员可以通过工具修改表结构。online: 线上环境开发人员不允许直接在线上环境进行数据库操作,如果...
  • 超硬数据结构学霸笔记,考试面试吹牛就靠它

    万次阅读 多人点赞 2021-03-26 11:11:21
    上次发操作系统笔记,很快浏览上万,这次数据结构比上次硬核的多哦,同样的会发超硬代码,关注吧。
  • 这篇文章主要是写在论文之前。写者的毕业论文中没有涉及该部分实现,因此使用博客的方式记录案例或者功能的实现。 世界疫情实时展示和之前的中国疫情实时展示差不多,同样是使用...重点一如既往的是数据的获取和压入。
  • 函数功能:判定数据(或特征)的分布情况 调用方法:plt.hist(x, bins=10, range=None, normed=False, weights=None, cumulative=False, bottom=None, histtype='bar', align='mid', orientation='vertical', rwidth=...
  • 准静态 C V 法测量硅表面态密度分布及数据处理 钱敏1 刘蓓1 辛煜2 11 苏州大学 电子信息学院 微电子系 江苏 苏州 215021 21 苏州大学 物理科学与技术学院 江苏 苏州 215006 摘要 表面态问题的研究是半导体材料 器件...
  • 数据可视实战:词云展示这个词云是大数据可视的重要方式,就是把文本中的关键词进行展示,就是很高大上![在这里插入图片描述]...
  • 标准正态分布密度函数公式

    千次阅读 2021-02-01 03:17:02
    展开全部标准正态分布密度函数公式:正态曲线呈钟型62616964757a686964616fe58685e5aeb931333366306532,两头低,中间高,左右对称因其曲线呈钟形,因此人们又经常称之为钟形曲线。若随机变量X服从一个数学期望为μ...
  • 散点图散点图非常适合展示两个变量之间的关系,因为你可以直接看到数据的原始分布。 如下面第一张图所示的,你还可以通过对组进行简单地颜色编码来查看不同组数据的关系。想要可视三个变量之间的关系? 没问题! ...
  • 使用Matplotlib进行数据可视(一)

    多人点赞 2021-05-31 17:07:47
    本篇为《使用Python进行数据分析》中介绍Matplotlib库的基础使用方法的第一篇,主要内容为使用Matplotlib库的一些常用技巧,绘画简易的线形图、散点图、密度图和等高线图等用图,以及进行可视异常处理等。...
  • 优化算法的事就是在 模型表征空间中找到模型评估指标最好的模型。这里就引出了什么是模型的表征空间,以及什么是评估指标了。只有正确地应用表征空间以及评估指标,才可以更好地优化模型 譬如SVM的模型表征空间...
  • WP155_R0_数据中心空间和功率密度需求的计算数据中心空间和功率密度需求的计算第 155 号白皮书版本 0作者:Neil Rasmussen摘要采用每平方米多少瓦特(或每平方英尺多少瓦特)来衡量数据中心的功率密度是一种传统的方法...
  • 我们的读者可以通过可视互动或其他数据使用方式来探寻一个故事的背后发生了什么,因此,数据可视至关重要。数据可视的目的其实就是直观地展现数据,例如让花费数小时甚至更久才能归纳的数据量,转化成一眼就能...
  • 本文总结了在数据分析和可视中最有用的 50 个 Matplotlib 图表。这些图表列表允许您使用 python 的 matplotlib 和 seaborn 库选择要显示的可视对象。...
  • 一图胜千言,使用Python的matplotlib库,可以快速创建高质量的图形。...我们推出一个新的系列教程:Python数据可视,针对初级和中级用户,将理论和示例代码相结合,分别使用matplotlib, seaborn, plotly等工具实...
  • 直方图、密度图 import numpy as np import pandas as pd import matplotlib.pyplot as plt %matplotlib inline #plt.hist(x, bins=10, range=None, normed=False, weights=None, cumulative=False, bottom=None, ...
  • 近期在处理分布式云数据中心建设项目中,查阅和翻看了不少华为、华三等厂家的资料,发现虚拟应是促使云数据中心走向成熟应用的最大变因,特此记录下来,并向大家分享虚拟相关概念。 高度虚拟在迅速地改变...
  • Matlab数据可视

    千次阅读 多人点赞 2021-03-28 15:16:38
    Matlab数据可视内容:数据探索与可视综合实践直方图、箱线图、散点图:直方图箱线图散点图 内容:数据探索与可视综合实践 利用Matlab对UCI‐Aliva分类数据集进行初步探索和分析: (1)针对第A、B、C三个类别...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 31,381
精华内容 12,552
关键字:

非参数核密度可以做数据标准化吗