gp参数估计 python python非参数估计

目录

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

• 优点
• 缺点

• 2、非参有什么用？
• 3、直方图估计

• 3.1 思想
• 3.2 优缺点

• 优点
• 缺点

• 3.3 代码实现

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

• matplotlib hist
• seaborn distplot

• 方式二：我的代码实现

• 定义直方图类
• 测试
• 尝试不同的区间宽度h

• 总结
• 参考资料：

开新坑，写一下非参课程的实验。以下是我的简易理解，不会写的太深入，可能会有错误疏漏，欢迎各位指出。若想深入了解非参的话，还是建议大家找找教材和论文进行阅读。

在开始讲直方图估计之前，先讲一些非参的基本概念。

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()

总结

这次尝试按着书上的公式去实现了直方图估计，代码写的可能有点冗余，不够精简，但基本上达到了我的期望水平，能够按自己意愿任意修改起始点$x_0$和窗宽$h$，完成实验作业。代码还有优化空间，但继续优化有点重复造轮子的感觉，真要实际工作中用还是调包快，我就不继续优化了~

参考资料：

[1] 【知乎】参数检验与非参数检验的区别与联系,各自的优缺点,如何选择？ [2] Härdle W, Müller M, Sperlich S, et al. Nonparametric and semiparametric models[M]. Berlin: Springer, 2004.