【非参】python实现直方图估计

Kevin Davis

已于 2023-02-03 16:51:48 修改

阅读量1.1k

点赞数 1

分类专栏：非参数模型文章标签： python 人工智能

于 2022-12-15 22:35:19 首次发布

本文链接：https://blog.csdn.net/weixin_44590417/article/details/128212830

版权

非参数模型专栏收录该内容

3 篇文章 0 订阅

订阅专栏

1、非参数模型和参数模型的区别是什么？

参数模型假设总体的分布类型已知，分布的参数未知，抽取样本数据来估计分布的参数，并用假设检验来检验未知参数的合理性。

然而实际生活中的很多数据它的分布类型我们是不知道的，因而参数模型的适用范围较窄，但统计检验的效力较强。

非参数模型则不需要对数据的总体分布进行假设，理论上适用于任何数据，因此适用范围非常广。

优点

适用范围广，较为稳健。

缺点

维数诅咒或维数灾难，对于高维度的数据（即变量数目较多），非参数模型计算非常缓慢。
在模型假设（即分布假设）正确的前提下，参数模型的表现更佳。

2、非参有什么用？

这个问题问得比较功利，但却可以帮助大家快速了解非参的作用（~~没用还学个锤子~~）。根据我的浅薄理解，目前非参主要有以下用途：

（1）估计分布函数 $F(x)=P(X\leq x)$ 。
（2）估计概率密度函数 $f (x)$ 。
（2）回归拟合及预测。
（3）做检验。

由于非参的检验效力相对较弱，常见的用途主要为前三条。拟合及预测的作用可以参考线性回归，根据给定的数据X来预测Y。非参还引入了连接函数（link function），增加了非线性，提高了模型的估计能力，但也带来了可解释性降低的问题。

（1）和（2）的作用相同，可以获取数据的分布情况信息。如何得知数据的分布情况呢？最简单的方法就是画个图，看看数据的分布情况，看数据在哪个区间的分布比较密集，哪个区间的分布比较稀疏。直方图估计主要的用途就是估计数据的分布。

3、直方图估计

3.1 思想

直方图估计的思想非常简单，选定一个初始点 $x_0$ ，将数据划分成等长的若干区间 $B_j = [x_0+(j-1)h, x_0 + h), j\in Z,$ 然后计算落在每个区间上的样本点个数（即频数） $n_j$ ，除以样本容量 $n$ 得到频率，再除以区间长度 $h$ 即可得到每个区间的估计概率密度函数 $f_j = \frac{n_j}{nh}.$
用数学公式表达为：
$\hat{f}_h(x) = \frac{1}{nh}\sum_{i=1}^{n} \sum_j I(X_i \in B_j) I(x\in B_j),$
其中，
$I(X_i \in B_J) = \begin{cases} 1 & if \quad X_i \in B_j, \\ 0 & otherwise. \end{cases}$
这里 $I(X_i \in B_J)$ 是示性函数，表示如果 $X_i$ 在 $B_j$ 区间内则为1，否则为0。这个 $\hat{f}_h(x)$ 的计算公式看起来比较复杂，其表达的意思是：先判断 $x$ 在哪个区间，然后加总求和这个区间内的所有样本，从而得到这个区间的样本频数，最后套用上文公式计算得到概率密度函数 $\hat{f}_h(x)$ 。

另一种更为简单的表达方式为：
$\hat{f}_h(m_j) = \frac{1}{nh} \# \{X_i \in [m_j - \frac{h}{2}, m_j + \frac{h}{2}) \}.$

$\#\{\}$ 号表示是满足 $\{\}$ 内表达式的样本个数。 $m_j$ 为第 $j$ 个区间的中心点。该方法先求解出落在以 $m_j$ 为中心点的区间的样本个数，然后除以 $nh$ 得到概率密度函数 $\hat{f}_h(m_j)$ 。

3.2 优缺点

优点

简单直接，容易理解
直方图能够非常方便地观察数据的分布情况

缺点

直方图的形状依赖于窗宽 $h$ 和原点 $x_0$
直方图估计的密度函数不是连续的，在跳跃处是不可导的
将连续数据离散化，每一个 $x\in [m_j - \frac{h}{2}, m_j + \frac{h}{2})$ 的估计 $f$ 都是 $\hat{f}_h (m_j)$

3.3 代码实现

方式一：懒人版，直接调包

matplotlib hist

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

np.random.seed(0)

n = 200
x = np.random.normal(0, 0.1, n)

plt.figure(figsize=(10, 7), dpi=80)
plt.rc('font', family='Times New Roman', weight = 'medium', size=14)  #设置英文字体
plt.hist(x,edgecolor='white')
plt.xlabel('x')
plt.ylabel('y')
plt.title(f'matplotlib hist')

plt.show()

在这里插入图片描述

seaborn distplot

plt.figure(figsize=(10, 7), dpi=80)
sns.distplot(a=x, kde=True)  #kde = True 增加密度线
plt.xlabel('x')
plt.ylabel('y')
plt.title(f'seaborn distplot')
plt.show()

在这里插入图片描述

方式二：我的代码实现

定义直方图类

class Histogram():
    def __init__(self, verbose=False):
        self.verbose = verbose #是否输出计算过程，默认不输出，default: False

    def check(self, x, lower, upper):
        '''
        示性函数，判断x是否落在区间[lower, upper) 里，是返回1，否返回0。
        '''
        if x >= lower and x < upper:
            return 1
        else:
            return 0

    def find_left(self, x, m_j, h):
        # 向左搜索
        center_left = []  #记录左侧的中点，记得要排序
        fh_list = []  #记录每一个箱子的密度估计值f_h

        lower = m_j - h / 2
        upper = m_j + h / 2

        while lower > x.min() - h:
            f = np.array([self.check(i, lower, upper) for i in x])  #示性函数，判断x的每一个值是否落在当前区间
            fh = f.sum() / (n * h)
            if self.verbose:
                print(f'Left: m_j = {m_j}, h = {h}, lower = {lower}, upper = {upper}, f.sum = {f.sum()}, fh = {fh}')
            center_left.append(m_j)
            fh_list.append(fh)
            
            m_j -= h
            lower = m_j - h / 2
            upper = m_j + h / 2

        center_left = [center_left[-i] for i in range(1, len(center_left)+1)]  #翻转
        fh_list = [fh_list[-i] for i in range(1, len(fh_list)+1)]

        return np.array(center_left), np.array(fh_list)

    def find_right(self, x, m_j, h):
        # 向右搜索
        center_right = []  #记录右侧的中点，不需要排序
        fh_list = []  #记录每一个箱子的密度估计值f_h

        lower = m_j - h / 2
        upper = m_j + h / 2

        while upper < x.max() + h:
            f = np.array([self.check(i, lower, upper) for i in x])  #示性函数，判断x的每一个值是否落在当前区间
            fh = f.sum() / (n * h)
            if self.verbose:
                print(f'Right: m_j = {m_j}, h = {h}, lower = {lower}, upper = {upper}, f.sum = {f.sum()}, fh = {fh}')
            center_right.append(m_j)
            fh_list.append(fh)
            
            m_j += h
            lower = m_j - h / 2
            upper = m_j + h / 2

        return np.array(center_right), np.array(fh_list)

    def fit(self, x, m_j, h):
        '''
        直方图拟合函数。
        Fit the data with hitogram nonparametric method.

        input:
        x: the data, must be a vector, ndarray, [n, ]. 数据，必须为一个向量，numpy数组。
        m_j: the inital center points, float or int. 初始数据中心。
        h: the window length, float or int. 窗宽参数。

        output:
        centers: each center of bins, ndarray, [n, ]. 每一个箱子的中心点m_j。
        fh_list: the estimated density corrensponding to each center, ndarray, [n,]. 每一个中心点对应的频率密度估计值。
        '''

        assert x.ndim == 1

        self.h = h
        self.x0 = m_j
        if m_j < x.min():
            centers, fh_list = self.find_right(x, m_j, h)
        elif m_j > x.max():
            centers, fh_list = self.find_left(x, m_j, h)
        else:
            c_left, fh_left = self.find_left(x, m_j, h)
            c_right, fh_right = self.find_right(x, m_j, h)
            centers = np.concatenate((c_left, c_right[1:]))
            fh_list = np.concatenate((fh_left, fh_right[1:]))
        
        self.centers = centers
        self.fh_list = fh_list

        return centers, fh_list

    def plot(self, saveflag=False):
        # 画图函数
        plt.figure(figsize=(10, 7), dpi=80)
        plt.rc('font', family='Times New Roman', weight = 'medium', size=14)  #设置英文字体
        # plt.plot(self.centers, self.fh_list, c='gold', marker='o')  
        plt.bar(self.centers, self.fh_list, width=self.h, edgecolor='white')  
        plt.xlabel('x')
        plt.ylabel('y')
        plt.title(f'h = {self.h}, x0 = {self.x0}')

        if saveflag:
            plt.savefig('histogram.png')
        else:
            plt.show()

测试

x0 = 0
h = 0.05
estimator = Histogram(verbose=True)
centers, fh_list = estimator.fit(x, x0, h)
estimator.plot()

在这里插入图片描述

尝试不同的区间宽度h

x0 = 0
h_list = [0.007, 0.05, 0.02, 0.1]
estimator = Histogram(verbose=False)
centers_set = []
fh_set = []
for h in h_list:
    centers, fh_list = estimator.fit(x, x0, h)
    centers_set.append(centers)
    fh_set.append(fh_list)

plt.figure(figsize=(16, 12), dpi=80)
plt.rc('font', family='Times New Roman', weight = 'medium', size=14)  #设置英文字体

for i in range(len(h_list)):
    ax = plt.subplot(2, 2, i+1)
    ax.bar(centers_set[i], fh_set[i], width=h_list[i], edgecolor='white')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_title(f'h = {h_list[i]}, x0 = {x0}')

plt.show()

在这里插入图片描述