精华内容
下载资源
问答
  • Python实现Parzen window classifier
    2021-04-15 16:29:11

    不指定近邻个数的近邻分类器 

    import numpy as np
    from sklearn import datasets
    import pandas as pd
    from sklearn.metrics.pairwise import rbf_kernel
    from sklearn.model_selection import StratifiedKFold
    from sklearn.metrics import accuracy_score
    from sklearn.neighbors import KNeighborsClassifier
    
    class PWC():
        def __init__(self,train_X,train_y,kerel = "rbf"):
            self.kernel = kerel
            self.X = train_X
            self.y = tr
    更多相关内容
  • 基于Matlab的Parzen Window函数。简单实用。
  • kernel density estimation是在概率论中用来估计未知的密度函数,属于非参数检验方法之一,由Rosenblatt (1955)和Emanuel Parzen(1962)提出,又名Parzen窗(Parzen window) 本文翻译自英国雷丁大学(Reading ...

    kernel density estimation是在概率论中用来估计未知的密度函数,属于非参数检验方法之一,由Rosenblatt (1955)和Emanuel Parzen(1962)提出,又名Parzen窗(Parzen window)
    本文翻译自英国雷丁大学(Reading University)Xia Hong老师的讲义材料

    [统计寓于生活.jpg)]

    1 概率密度函数

    连续概率函数 p ( x ) p(x) p(x)的数学定义满足以下特性:

    1. x x x介于 a a a b b b两点之间的概率为 P ( a < x < b ) = ∫ a b p ( x ) d x P(a<x<b)=\int_a^bp(x)dx P(a<x<b)=abp(x)dx
    2. x是非负实数
    3. 概率函数的积分: ∫ − ∞ ∞ p ( x ) d x = 1 \int_{-\infty}^\infty p(x)dx=1 p(x)dx=1

    最常见的概率函数是高斯函数(Gaussian Function,又称为正态分布): p ( x ) = 1 2 π σ exp ⁡ ( − ( x − c ) 2 2 σ 2 ) p(x)=\frac{1}{\sqrt{2\pi\sigma}}\exp(-\frac{(x-c)^2}{2\sigma^2}) p(x)=2πσ 1exp(2σ2(xc)2),这里 c c c是均值, σ \sigma σ是标准差

    Guassian-Function.

    拓展到对于向量 x \mathbf{x} x,非负函数 p ( x ) p(\mathbf{x}) p(x)有以下特性:

    1. x \mathbf{x} x在区域 R R R里的概率为 P = ∫ R p ( x ) d x P=\int_R p(\mathbf{x})d\mathbf{x} P=Rp(x)dx
    2. 概率函数的积分为 ∫ p ( x ) d x = 1 \int p(\mathbf{x})d\mathbf{x}=1 p(x)dx=1

    2D-Gaussian-

    2 密度估计

    密度估计:给定的一系列数量为 n n n的样本 x 1 , ⋅ ⋅ ⋅ , x n \mathbf{x}_1, \cdot\cdot\cdot,\mathbf{x}_n x1,,xn,可以估计密度函数 p ( x ) p(\mathbf{x}) p(x),从而根据任意新样本 x \mathbf{x} x可以得到输出 p ( x ) p(\mathbf{x}) p(x)
    大部分未知密度函数估计方法的基本思想都很简单,主要是依赖于样本落在区域 R R R的概率 P P P,即有 P = ∫ R p ( x ) d x P=\int_R p(\mathbf{x})d\mathbf{x} P=Rp(x)dx
    假设区域 R R R很小, P ( x ) P(\mathbf{x}) P(x)在区域内波动很小,上式可以写做 P = ∫ R p ( x ) d x ≈ p ( x ) ∫ R d x = p ( x ) V P=\int_R p(\mathbf{x})d\mathbf{x}\approx p(\mathbf{x})\int_R d\mathbf{x}=p(\mathbf{x})V P=Rp(x)dxp(x)Rdx=p(x)V,这里 V V V是区域 R R R的“量”(二维即为面积)
    从另一方面看,假设 n n n个样本 x 1 , ⋅ ⋅ ⋅ , x n \mathbf{x}_1, \cdot\cdot\cdot,\mathbf{x}_n x1,,xn都是独立且服从概率密度函数 p ( x ) p(\mathbf{x}) p(x),且 n n n个样本中有 k k k个落在区域 R R R里面,则有 P = k / n P=k/n P=k/n,因此 p ( x ) p(\mathbf{x}) p(x)的估计式为 p ( x ) = k / n V p(\mathbf{x})=\frac{k/n}{V} p(x)=Vk/n

    3 Parzen窗密度估计

    考虑 R R R是中心在 x \mathbf{x} x的超立方体(例如二维平面),令 h h h为超立方体的边缘长度,所以对于二维平面有有 V = h 2 V=h^2 V=h2,对于三维立体有 V = h 3 V=h^3 V=h3

    Parzen-Window

    引入 ϕ ( x i − x h ) = { 1 ∣ x i k − x k ∣ h < = 1 / 2 , k = 1 , 2 0 o t h e r w i s e \phi(\frac{\mathbf{x}_i-\mathbf{x}}{h})=\left\{ \begin{aligned} 1\quad& \frac{|x_{ik}-x_{k}|}{h}<=1/2, k=1,2 \\ 0\quad& otherwise \end{aligned} \right. ϕ(hxix)=10hxikxk<=1/2,k=1,2otherwise
    Parzen概率密度公式(二维)为 p ( x ) = k / n V = 1 n ∑ i = 1 n 1 h 2 ϕ ( x i − x h ) p(\mathbf{x})=\frac{k/n}{V}=\frac{1}{n}\sum_{i=1}^n {\frac{1}{h^2}\phi(\frac{\mathbf{x}_i-\mathbf{x}}{h})} p(x)=Vk/n=n1i=1nh21ϕ(hxix) ϕ ( x i − x h ) \phi(\frac{\mathbf{x}_i-\mathbf{x}}{h}) ϕ(hxix)即为窗函数
    我们归纳这个思想并拓展到其他Parzen窗密度估计法中
    例如,如果使用高斯函数,对于一维有: p ( x ) = 1 n ∑ i = 1 n 1 2 π σ exp ⁡ ( − ( x i − x ) 2 2 σ 2 ) p(x)=\frac{1}{n}\sum_{i=1}^{n}{\frac{1}{\sqrt{2\pi\sigma}}\exp(-\frac{(x_i-x)^2}{2\sigma^2})} p(x)=n1i=1n2πσ 1exp(2σ2(xix)2),这是对 n n n个将数据点作为中心的高斯函数的简单求平均,公式中的 σ \sigma σ需要再做确定


    例子
    给定一个系列的5个数据点 x 1 = 2 x_1=2 x1=2 x 2 = 2.5 x_2=2.5 x2=2.5 x 3 = 3 x_3=3 x3=3 x 4 = 1 x_4=1 x4=1 x 5 = 6 x_5=6 x5=6,参数 σ = 1 \sigma=1 σ=1,中心 x = 3 x=3 x=3的高斯函数作为窗函数,求出Parzen概率密度估计(pdf)
    解答
    1 2 π exp ⁡ ( − ( x 1 − x ) 2 2 ) = 1 2 π exp ⁡ ( − ( 2 − 3 ) 2 2 ) = 0.2420 \frac{1}{\sqrt{2\pi}}\exp(-\frac{(x_1-x)^2}{2})=\frac{1}{\sqrt{2\pi}}\exp(-\frac{(2-3)^2}{2})=0.2420 2π 1exp(2(x1x)2)=2π 1exp(2(23)2)=0.2420
    1 2 π exp ⁡ ( − ( x 2 − x ) 2 2 ) = 1 2 π exp ⁡ ( − ( 2.5 − 3 ) 2 2 ) = 0.3521 \frac{1}{\sqrt{2\pi}}\exp(-\frac{(x_2-x)^2}{2})=\frac{1}{\sqrt{2\pi}}\exp(-\frac{(2.5-3)^2}{2})=0.3521 2π 1exp(2(x2x)2)=2π 1exp(2(2.53)2)=0.3521
    1 2 π exp ⁡ ( − ( x 3 − x ) 2 2 ) = 1 2 π exp ⁡ ( − ( 3 − 3 ) 2 2 ) = 0.3989 \frac{1}{\sqrt{2\pi}}\exp(-\frac{(x_3-x)^2}{2})=\frac{1}{\sqrt{2\pi}}\exp(-\frac{(3-3)^2}{2})=0.3989 2π 1exp(2(x3x)2)=2π 1exp(2(33)2)=0.3989
    1 2 π exp ⁡ ( − ( x 4 − x ) 2 2 ) = 1 2 π exp ⁡ ( − ( 1 − 3 ) 2 2 ) = 0.0540 \frac{1}{\sqrt{2\pi}}\exp(-\frac{(x_4-x)^2}{2})=\frac{1}{\sqrt{2\pi}}\exp(-\frac{(1-3)^2}{2})=0.0540 2π 1exp(2(x4x)2)=2π 1exp(2(13)2)=0.0540
    1 2 π exp ⁡ ( − ( x 5 − x ) 2 2 ) = 1 2 π exp ⁡ ( − ( 6 − 3 ) 2 2 ) = 0.0044 \frac{1}{\sqrt{2\pi}}\exp(-\frac{(x_5-x)^2}{2})=\frac{1}{\sqrt{2\pi}}\exp(-\frac{(6-3)^2}{2})=0.0044 2π 1exp(2(x5x)2)=2π 1exp(2(63)2)=0.0044
    因此, p ( x = 3 ) = ( 0.2420 + 0.3521 + 0.3989 + 0.0540 + 0.0044 ) / 5 = 0.2103 p(x=3)=(0.2420 + 0.3521 + 0.3989+0.0540 + 0.0044)/5 = 0.2103 p(x=3)=(0.2420+0.3521+0.3989+0.0540+0.0044)/5=0.2103


    下面用图形化语言表示Parzen窗,每个数据点密度函数(虚线)对于最终的概率密度函数(实线)有相同的贡献度

    lines-of-data-points

    line-of-pdf.

    展开全文
  • Parzen Window and Likelihood

    千次阅读 2018-01-30 23:05:04
    Kernel density estimation via the Parzen-Rosenblatt window method Generative Adversarial Nets code Nice explanation 对Ian提到的采用Gassian Parzen Window 去拟合G产生的样本并估计对数似然性(即给予观察...

    Kernel density estimation via the Parzen-Rosenblatt window method

    Generative Adversarial Nets code

    Nice explanation

    对Ian提到的采用Gassian Parzen Window 去拟合G产生的样本并估计对数似然性(即给予观察样本,判断模型正确的可能性)如何实现的原理很感兴趣?
    于是找到源代码parzen_ll.py来看,
    但是苦于并没有学过theano及pylearn2所以看起来不是明白,本来想尝试安装py2learn来测试,结果看了官网的安装要求,请卸载python,清除所有组件。这个代价有点大
    ,毕竟我是要做Tensorflow的。因此只能破罐破摔来,硬着头皮看代码,有些东西只能靠google,看了一个大概。首先来说一下Parzen Window

    2 Defining the Region Rn

    Parzen windown 是一种应用广泛的非参数估计方法,从采样的样本 p(xn) p ( x n ) 中估计概率密度函数 p(x) p ( x )
    这种方法最基本的思想是给定一个特定的区域(Window)对落入其中的样本进行计数,可以得到样本落入该区域的概率大约为:

    p(x)=knn p ( x ) = k n n

    从数学的角度来求取在区域R里k个观察样本的概率,我们可以考虑一个Binomial(二项式分布):

    pk=[nk]pk(1p)nk p k = [ n k ] ⋅ p k ⋅ ( 1 − p ) n − k

    然后在二项式分布假设的情况下, 所以我们可以求得其均值为:

    E[k]=np E [ k ] = n ⋅ p

    在连续性假设条件下:

    p(x)=Rf(x)dx=f(x)V p ( x ) = ∫ R f ( x ) d x = f ( x ) ⋅ V

    其中v为区域R的体积,我们可以对公式进行变形,可以得到:

    f(x)=p(x)V=k/nV f ( x ) = p ( x ) V = k / n V

    上面简单的公式可以让计算给定点 x x 的概率密度,通过计算有多少点落入确定性区域。

    确定区域的两种不同的方法

    一种方法是固定体积(volume)如下图,原理是取一个固定大小的区域R,观察有多少样本落入其中。

    这里写图片描述
    另外一种方法是样本数目K固定;k近邻算法就是这个原理,即是,对于一定大小的样本总数,我们取能框住k个样本的区域R
    这里写图片描述

    3D- hypercubes

    这里写图片描述

    The window function

    如上面所示,我们可以可视化区域R,也就很容易计算出多少点在区域内,数学表达上是:

    ϕ(u)=[10|uj|1/2;j=1,2,...,dotherwise ϕ ( u ) = [ 1 | u j | ≤ 1 / 2 ; j = 1 , 2 , . . . , d 0 o t h e r w i s e

    扩展概念(立方体核,hyercube window)我们可以得到下面的概念:
    通过下面的公式可以获得落入该区域的样本个数Kn

    kn=i=1nϕ(xxi,hn) k n = ∑ i = 1 n ϕ ( x − x i , h n )

    式中:

    u=(xxihn) u = ( x − x i h n )

    也就是:

    ϕ(xxi,hn)=[10xxihn1/2;j=1,2,...,dotherwise ϕ ( x − x i , h n ) = [ 1 | x − x i h n | ≤ 1 / 2 ; j = 1 , 2 , . . . , d 0 o t h e r w i s e

    进而我们可以的得到Rn的概率

    而将将简单的立方体核函数替换掉,采用高斯核:

    ϕ(xxi,hn)=1(2π)d/2|Σ|1/2exp((xxi)TΣ1(xxi)2h2) ϕ ( x − x i , h n ) = 1 ( 2 π ) d / 2 | Σ | 1 / 2 exp ⁡ ( − ( x − x i ) T Σ − 1 ( x − x i ) 2 h 2 )

    效果图

    Parzen window estimation

    基于上述公式,我们可以得到:

    pn(x)=1ni=1n1hdϕ[xxi,hn] p n ( x ) = 1 n ∑ i = 1 n 1 h d ϕ [ x − x i , h n ]

    式中有:

    Vn=hd V n = h d

    k=i=1nϕ[xxi,hn] k = ∑ i = 1 n ϕ [ x − x i , h n ]

    Parzen-window技术的关键性参数

    这有两个关键性参数:
    - 核函数(kernel function)
    - 核宽度(window width)

    核密度估计需要满足如下条件:
    1) 有限值非负密度函数:

    f(y,h)dy=1 ∫ f ( y , h ) d y = 1

    2) 对于核宽度的要求

    limnnhd(n)= lim n → ∞ ⁡ n h d ( n ) = ∞

    limnhd=0 lim n → ∞ ⁡ h d = 0

    limnkn= lim n → ∞ ⁡ k n = ∞

    一式原因在于:
    limnpn(x)=p(x) lim n → ∞ ⁡ p n ( x ) = p ( x )

    limnknnV=p(x) lim n → ∞ ⁡ k n n ⋅ V = p ( x )

    构造Parzen window估计函数

    def parzen_estimation(mu, sigma, mode='gauss'):
        """
        Implementation of a parzen-window estimation
        Keyword arguments:
            x: A "nxd"-dimentional numpy array, which each sample is
                      stored in a separate row (=training example)
            mu: point x for density estimation, "dx1"-dimensional numpy array
            sigma: window width
        Return the density estimate p(x)
        """
    
        def log_mean_exp(a):
            max_ = a.max(axis=1)
            return max_ + np.log(np.exp(a - np.expand_dims(max_, axis=0)).mean(1))
    
        def gaussian_window(x, mu, sigma):
            a = (np.expand_dims(x, axis=1) - np.expand_dims(mu, axis=0)) / sigma
            b = np.sum(- 0.5 * (a ** 2), axis=-1)
            E = log_mean_exp(b)
            Z = mu.shape[1] * np.log(sigma * np.sqrt(np.pi * 2))
            return np.exp(E - Z)
    
        def hypercube_kernel(x, mu, h):
            n, d = mu.shape
            a = (np.expand_dims(x, axis=1) - np.expand_dims(mu, axis=0)) / h
            b = np.all(np.less(np.abs(a), 1/2), axis=-1)
            kn = np.sum(b.astype(int), axis=-1)
            return kn / (n * h**d)
    
        if mode is 'gauss':
            return lambda x: gaussian_window(x, mu, sigma)
        elif mode is 'hypercube':
            return lambda x: hypercube_kernel(x, mu, h=sigma)

    将Parzen Window 方法应用于高斯数据集

    我们构造如下数据集:

    f(x)=1(2π)d/2|Σ|1/2exp((xxi)TΣ1(xxi)2) f ( x ) = 1 ( 2 π ) d / 2 | Σ | 1 / 2 exp ⁡ ( − ( x − x i ) T Σ − 1 ( x − x i ) 2 )

    设定:
    u=[00]Σ=[1001] u = [ 0 0 ] Σ = [ 1 0 0 1 ]

    import numpy as np
    
    # Generate 10,000 random 2D-patterns
    mu_vec = np.array([0,0])
    cov_mat = np.array([[1,0],[0,1]])
    x_2Dgauss = np.random.multivariate_normal(mu_vec, cov_mat, 10000)
    
    print(x_2Dgauss.shape)

    这里写图片描述

    不同核函数对于pdf的估计效果不同

    这里写图片描述

    可以看到高斯核函数的概率密度估计过渡更加自然,而超维度的高斯估计则虽然相比较精度较高一些但是毛刺比较严重。

    模型[0, 0]
    真实值0.15915494309189535
    hypercube(h=0.3)0.15444444
    Gauss(h=0.3)0.14109129
    import matplotlib.pyplot as plt
        from mpl_toolkits.mplot3d import Axes3D
        from matplotlib import cm
        from matplotlib.ticker import LinearLocator, FormatStrFormatter
    
        # Make data
        mu_vec = np.array([0, 0])
        cov_mat = np.array([[1, 0], [0, 1]])
        x_2Dgauss = np.random.multivariate_normal(mu_vec, cov_mat, 10000)
    
        # Pdf estimation
        def pdf_multivaraible_gauss(x, mu, cov):
            part1 = 1 / ((2 * np.pi) ** (len(mu)/2) * (np.linalg.det(cov)**(1/2)))
            part2 = (-1/2) * (x-mu).T.dot(np.linalg.inv(cov)).dot((x-mu))
            return float(part1 * np.exp(part2))
    
        ph = parzen_estimation(x_2Dgauss, 0.3, mode='hypercube')
        pg = parzen_estimation(x_2Dgauss, 0.3, mode='gauss')
        print(pdf_multivaraible_gauss(np.array([[0], [0]]), np.array([[0], [0]]), cov_mat))
        print(ph([[0, 0]]))
        print(pg([[0, 0]]))
    
        x = np.linspace(-5, 5, 100)
        x, y = np.meshgrid(x, x)
        zg = []
        zt = []
        zh = []
        for i, j in zip(x.ravel(), y.ravel()):
            zg.append(pg([[i, j]]))
            zh.append(ph([[i, j]]))
            zt.append(pdf_multivaraible_gauss(np.array([[i], [j]]), np.array([[0], [1]]), cov_mat))
        zg = np.asarray(zg).reshape(100, 100)
        zh = np.asarray(zh).reshape(100, 100)
        zt = np.asarray(zt).reshape(100, 100)
    
        # Plot the surface
        fig = plt.figure(figsize=(18, 6))
        ax1 = fig.add_subplot(131, projection='3d')
        surf = ax1.plot_surface(x, y, zg, rstride=1, cstride=1,
                                cmap=cm.coolwarm, antialiased=False)
        # Customize the z axis
        ax1.zaxis.set_major_locator(LinearLocator(10))
        ax1.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
        ax1.set(zlim=[0, 0.2], xlabel='x', ylabel='y', zlabel='p(x,y)',
                title="Gauss Bivariate Gaussian density")
    
        # Re-Plot
        ax2 = fig.add_subplot(132, projection='3d')
        ax2.plot_surface(x, y, zt, rstride=1, cstride=1, cmap=cm.coolwarm, antialiased=False)
    
        # Customize the z axis
        ax2.zaxis.set_major_locator(LinearLocator(10))
        ax2.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
        ax2.set(zlim=[0, 0.2], xlabel='x', ylabel='y', zlabel='p(x,y)',
                title="True Bivariate Gaussian density")
    
        # Re-Plot
        ax3 = fig.add_subplot(133, projection='3d')
        ax3.plot_surface(x, y, zh, rstride=1, cstride=1, cmap=cm.coolwarm, antialiased=False)
    
        # Customize the z axis
        ax3.zaxis.set_major_locator(LinearLocator(10))
        ax3.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
        ax3.set(zlim=[0, 0.2], xlabel='x', ylabel='y', zlabel='p(x,y)',
                title="Hypercube Bivariate Gaussian density")
    
        # Add a color bar which maps values to colors.
        fig.savefig("./Gauss_kernel_{}.png".format('Bivariate_Gaussian'))
        plt.show()

    不同核宽度对pdf估计的影响

    这里写图片描述
    这里写图片描述

    可以看到当核宽度非常小(0.01)时, pdf已经不可控, 而当h=0.5, 估计效果又有所下降。因此如何能根据数据样本来确定核宽度大小,是一个非常重要的问题。

        import matplotlib.pyplot as plt
        from mpl_toolkits.mplot3d import Axes3D
        from matplotlib import cm
        from matplotlib.ticker import LinearLocator, FormatStrFormatter
    
        # Make data
        mu_vec = np.array([0, 0])
        cov_mat = np.array([[1, 0], [0, 1]])
        x_2Dgauss = np.random.multivariate_normal(mu_vec, cov_mat, 10000)
    
        # Pdf estimation
        def pdf_multivaraible_gauss(x, mu, cov):
            part1 = 1 / ((2 * np.pi) ** (len(mu)/2) * (np.linalg.det(cov)**(1/2)))
            part2 = (-1/2) * (x-mu).T.dot(np.linalg.inv(cov)).dot((x-mu))
            return float(part1 * np.exp(part2))
    
        pg1 = parzen_estimation(x_2Dgauss, 0.01, mode='gauss')
        pg2 = parzen_estimation(x_2Dgauss, 0.1, mode='gauss')
        pg3 = parzen_estimation(x_2Dgauss, 0.5, mode='gauss')
        x = np.linspace(-5, 5, 100)
        x, y = np.meshgrid(x, x)
        zg = []
        zh = []
        zt = []
        for i, j in zip(x.ravel(), y.ravel()):
            zg.append(pg1([[i, j]]))
            zh.append(pg2([[i, j]]))
            zt.append(pg3([[i, j]]))
        zg = np.asarray(zg).reshape(100, 100)
        zh = np.asarray(zh).reshape(100, 100)
        zt = np.asarray(zt).reshape(100, 100)
    
        # Plot the surface
        fig = plt.figure(figsize=(18, 6))
        ax1 = fig.add_subplot(131, projection='3d')
        surf = ax1.plot_surface(x, y, zg, rstride=1, cstride=1,
                                cmap=cm.coolwarm, antialiased=False)
        # Customize the z axis
        ax1.zaxis.set_major_locator(LinearLocator(10))
        ax1.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
        ax1.set(zlim=[0, 0.2], xlabel='x', ylabel='y', zlabel='p(x,y)',
                title="Bivariate Gaussian density-0.01")
    
        # Re-Plot
        ax2 = fig.add_subplot(132, projection='3d')
        ax2.plot_surface(x, y, zh, rstride=1, cstride=1, cmap=cm.coolwarm, antialiased=False)
    
        # Customize the z axis
        ax2.zaxis.set_major_locator(LinearLocator(10))
        ax2.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
        ax2.set(zlim=[0, 0.2], xlabel='x', ylabel='y', zlabel='p(x,y)',
                title="Bivariate Gaussian density-0.1")
    
        # Re-Plot
        ax3 = fig.add_subplot(133, projection='3d')
        ax3.plot_surface(x, y, zt, rstride=1, cstride=1, cmap=cm.coolwarm, antialiased=False)
    
        # Customize the z axis
        ax3.zaxis.set_major_locator(LinearLocator(10))
        ax3.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
        ax3.set(zlim=[0, 0.2], xlabel='x', ylabel='y', zlabel='p(x,y)',
                title="Bivariate Gaussian density-0.5")
    
        # Add a color bar which maps values to colors.
        fig.savefig("./Gauss_kernel_{}.png".format('Gauss_width'))
        plt.show()

    已知真实分布情况下

    我们可以计算出最好的核宽度为,通过误差图来观察
    这里写图片描述

    可以看到 h = 0.21529369369369367

    当然大多数情况下并不知道真实的概率分布,那可以通过最大似然法来选择核宽度

    首先介绍似然性

    Likelihood

    可能性其实就是所有采的样本在你所建立的模型中概率的乘积,让这个乘积最大就是最大似然法,而为了计算方便,采用对数的方法,将
    乘积运算转换为加法运算。

    对parzen_ll的理解

    终于到了主题:

    def theano_parzen(mu, sigma):
        """
        Credit: Yann N. Dauphin
        """
    
        x = T.matrix()
        mu = theano.shared(mu)
        a = ( x.dimshuffle(0, 'x', 1) - mu.dimshuffle('x', 0, 1) ) / sigma
        E = log_mean_exp(-0.5*(a**2).sum(2))
        Z = mu.shape[1] * T.log(sigma * numpy.sqrt(numpy.pi * 2))
    
        return theano.function([x], E - Z)

    这个函数用于创造一个parzen的计算器,其中x的大小为2D向量,通过dimshuffule来变为Ax1xB 这样的向量,mu为图像的所有像素的值,sigma对于所有维都是固定的均为sigma。通过创造一个log_mean_exp函数计算平均值

    def log_mean_exp(a):
        """
        Credit: Yann N. Dauphin
        """
    
        max_ = a.max(1)
    
        return max_ + T.log(T.exp(a - max_.dimshuffle(0, 'x')).mean(1))

    定义好Parzen函数之后,通过get_nll函数获得似然值,其中的操作是每次输入一个batch_size,最后对各个batch_size求取平均值

    def get_nll(x, parzen, batch_size=10):
        """
        Credit: Yann N. Dauphin
        """
    
        inds = range(x.shape[0])
        n_batches = int(numpy.ceil(float(len(inds)) / batch_size))
    
        times = []
        nlls = []
        for i in range(n_batches):
            begin = time.time()
            nll = parzen(x[inds[i::n_batches]])
            end = time.time()
            times.append(end-begin)
            nlls.extend(nll)
    
            if i % 10 == 0:
                print i, numpy.mean(times), numpy.mean(nlls)
    
        return numpy.array(nlls)

    其中可以追查mu的来源,mu ==samples ,于是在下面我们可以发现出mu的蛛丝马迹,reshape函数将三个通道的图片扁平为一维向量,而整个空间定格为整个图片的像素点数。

    samples = model.generator.sample(args.num_samples).eval()
        output_space = model.generator.mlp.get_output_space()
        if 'Conv2D' in str(output_space):
            samples = output_space.convert(samples, output_space.axes, ('b', 0, 1, 'c'))
            samples = samples.reshape((samples.shape[0], numpy.prod(samples.shape[1:])))
        del model
        gc.collect()

    值得一提的是为了获得较好的估计值,他采用验证集来确定最大的似然值对应的sigma,验证集采用的是50000-60000的MNIST字体。

    def cross_validate_sigma(samples, data, sigmas, batch_size):
    
        lls = []
        for sigma in sigmas:
            print sigma
            parzen = theano_parzen(samples, sigma)
            tmp = get_nll(data, parzen, batch_size = batch_size)
            lls.append(numpy.asarray(tmp).mean())
            del parzen
            gc.collect()
    
        ind = numpy.argmax(lls)
        return sigmas[ind]

    Parzen Window的公式大致对的上,明确了输入(对于MNIST)是60000x1x784的一个批次,大约10个,主要可能考虑到对于高达784的多变量高斯分布计算量必然很大。
    然后考虑到likelihood的计算公式4

    将Parzen Window估计概率并应用于求解最大似然法,目前来说,对求解似然性于不是最好的办法,但是也没有更好的办法。其基本思路大致是,产生的样本通过高斯核Parzen窗口法计算出一个概率模型Pg(784维的高斯分布),然后估计测试样本的概率,从而计算出在该分布下的对数似然性,方差参数采用交叉验证来获得。

    最大似然法

    import matplotlib.pyplot as plt
        from mpl_toolkits.mplot3d import Axes3D
        from matplotlib import cm
        from matplotlib.ticker import LinearLocator, FormatStrFormatter
    
        # Make data
        np.random.seed(2017)
        mu_vec = np.array([0, 0])
        cov_mat = np.array([[1, 0], [0, 1]])
        x_2Dgauss = np.random.multivariate_normal(mu_vec, cov_mat, 10000)
        x_test = np.random.multivariate_normal(mu_vec, cov_mat, 1000)
    
        pg = parzen_estimation(x_2Dgauss, 0.01, mode='gauss')
        sigma, nll = cross_validate_sigma(x_2Dgauss, x_test, np.linspace(1e-4, 1, 20), batch_size=100)
        fig = plt.figure(figsize=(6, 6))
        ax = fig.gca()
        ax.plot(np.linspace(1e-4, 1, 20), nll,  '-*r')
        ax.set(xlabel='log likelihood', ylabel='\sigma value', title='NLL method ')
        print(sigma)
        fig.savefig('./{}.png'.format('NLL_method'))
        plt.show()
        # Pdf estimation
        def pdf_multivaraible_gauss(x, mu, cov):
            part1 = 1 / ((2 * np.pi) ** (len(mu)/2) * (np.linalg.det(cov)**(1/2)))
            part2 = (-1/2) * (x-mu).T.dot(np.linalg.inv(cov)).dot((x-mu))
            return float(part1 * np.exp(part2))
    
    
        tr = pdf_multivaraible_gauss(np.array([[0], [0]]), np.array([[0], [0]]), cov_mat)
        errs = []
        sigmas = np.linspace(1e-4, 1, 1000)
        for sigma in sigmas:
            pg = parzen_estimation(x_2Dgauss, sigma, mode='gauss')
            err = np.abs(pg(np.array([[0, 0]]))-tr)
            errs.append(err)
        fig = plt.figure(figsize=(6, 6))
        ax = fig.gca()
        ax.plot(sigmas, errs, '-*b')
        ax.set(xlabel='AE', ylabel='\sigma value', title='AE method ')
        ind = np.argmin(errs)
        print(sigmas[ind])
        fig.savefig('./{}.png'.format('AE_method'))
        plt.show()

    最大似然法

    可以看到在h = 0.21060526315789474, 附近似然性取得最大值

    方法最优值
    最大似然法0.21060526315789474
    误差法0.21529369369369367

    可以看到极大似然法可以有效估计核宽度。

    Parzen window 技术的优缺点:

    Parzen window技术作为一种非参数方法,有着和其他方法一样的缺点,就是每次估计需要整个样本参与运算,而并非是参数估计,只是带入参数计算。
    几个挑战:
    1) 数据样本的大小, Parzen window的计算复杂度为 n2d n 2 d , 其中 n n 为训练样本数目, 而d为样本维度。

    2) 合适的核宽度的选择, 一般来说:

    hn1n h n ∝ 1 n

    hn1logn h n ∝ 1 l o g ( n )

    3) 核函数的选择

    4) 应用于两个方面:Bayes估计与互信息的求取

    展开全文
  • Parzen window

    2019-08-05 19:34:14
    同时可以对窗口函数做一定的泛化,就有其他的Parzen window密度估计方法。 例如在1-D的情况下使用Gaussian函数: 这种方法就相当于将n个点为中心的高斯函数计算平均。其中标准差 需要预先设定。   例子:...

    转载自:http://blog.sina.com.cn/s/blog_679e13290101cpr1.html

    主要参考资料: http://www.personal.rdg.ac.uk/~sis01xh/teaching/CY2D2/Pattern2.pdf

    在数学上一个连续概率密度函数p(x)的需满足以下的条件:
    1、x在a和b之间的概率为:
    Parzen <wbr>window概率密度估计


    2、对所有的x,p(x)非负
    3、p(x)的积分值为1
    Parzen <wbr>window概率密度估计


    最经常使用的概率密度函数就是高斯函数(正态分布)
    Parzen <wbr>window概率密度估计

    Parzen <wbr>window概率密度估计
    将一维的情况扩展到多维,现在的 x就是一个向量,p( x)也需要满足下列条件:
    1、在一个区域 Rx的概率为
    Parzen <wbr>window概率密度估计
    2、概率密度函数的积分值为1
    Parzen <wbr>window概率密度估计

    Parzen <wbr>window概率密度估计
    密度估计
    给点n个数据样本x1,x2,....,xn,我们可以估计概率密度函数p(x),对于新的样本x就可以计算出相应的p(x).这个过程就是密度估计。
    密度估计的基础是:一个向量x落入到区域R的概率为
    Parzen <wbr>window概率密度估计
    假设R非常小,所以p(x)的变化也很小,上面的公式就改写为:
    Parzen <wbr>window概率密度估计

    其中V是R的“体积”
     
    另一方面,假设x1,...,xn是根据密度函数p(x)独立取的n个样本点,其中有k个样本点落入到区域R中,关于R的概率就为:
    Parzen <wbr>window概率密度估计
    这样就可以得到一个p(x)的估计函数:
    Parzen <wbr>window概率密度估计
     
    Parzen window密度估计
    假设R是以x为中心的超立方体,h为这个超立方体的边长,在2-D的方形中有V=h*h,3-D的立方体中有V=h^3。
    Parzen <wbr>window概率密度估计
    Parzen <wbr>window概率密度估计
    给定上面的公式,表示的是Xi是否落在方形中。
    Parzen概率密度估计公式的表示如下:
    Parzen <wbr>window概率密度估计
    其中 Parzen <wbr>window概率密度估计被称作窗口函数(window function)。
    同时可以对窗口函数做一定的泛化,就有其他的Parzen window密度估计方法。
    例如在1-D的情况下使用Gaussian函数:
    Parzen <wbr>window概率密度估计
    这种方法就相当于将n个点为中心的高斯函数计算平均。其中标准差 Parzen <wbr>window概率密度估计需要预先设定。
     
    例子:
    给定五个点:x1=2, x2=2.5, x3=3, x4=1, x5=6, 计算x=3位置的Parzen概率密度函数,采用 Parzen <wbr>window概率密度估计的高斯函数作为window function。
    计算过程如下:
    Parzen <wbr>window概率密度估计

    Parzen <wbr>window概率密度估计

    Parzen <wbr>window概率密度估计 
    采用图形的方式进行显示,并假设上面的5个点对整个密度函数做出相等的贡献:
    Parzen <wbr>window概率密度估计
    采用Parzen Window对这个五个点估计得到的概率密度函数为:
    Parzen <wbr>window概率密度估计



    展开全文
  • Parzen_parzen_parzenwindowmatlab_parzenwindow.zip
  • Parzen_parzen_parzenwindowmatlab_parzenwindow_源码.zip
  • import numpy as np from sklearn.base import BaseEstimator, ClassifierMixin from sklearn.utils import check_random_state, check_array from sklearn.metrics.pairwise import rbf_kernel ...
  • 计算过程如下: 采用图形的方式进行显示,并假设上面的5个点对整个密度函数做出相等的贡献: 采用Parzen Window对这个五个点估计得到的概率密度函数为: 补充: 自适应带宽选取 参数与非参数密度估计
  • Parzen window 概率密度估计

    千次阅读 2018-01-25 17:09:27
    非参数估计:已知样本所属的类别,但未知总体概率密度函数的形式,要求我们直接推断概率密度函数本身。非参数估计的方法主要有:直方图法、核方法。...window概率密度估计" title="Parzen window概率密度估计" sty
  • 这种假设尤其适用于样本定义域离散的情况,而当定义域连续情况下,则需要进一步设计求解规则,比如说Parzen window和K-nearest。   Parzen Windows 是一种通过窗口来估计窗口中心点概率密度的方法。假设窗口内的点...
  • https://blog.csdn.net/unixtch/article/details/78556499
  • PARZEN窗和K近邻算法的python实现。 现实生活中常常会有这样的问题:缺乏足够的先验知识,因此难以人工标注类别或进行人工类别标注的成本太高。很自然地,我们希望计算机能代我们完成这些工作,或至少提供一些帮助...
  • 背景 我们有非常多的数据,想求出该数据的概率密度函数。我们不知道这些数据服从哪一种分布,所以采用概率密度的无参数估计。 kernel density estimator 前提假设 假设采集到的数据 x∈RD\mathbf{x} \in \mathbb{R}^...
  • [1] O. Chapelle, "Active Learning for Parzen Window Classifier", Proceedings of the Tenth International Workshop Artificial Intelligence and Statistics, 2005.
  •  Parzen窗口密度估计是一种基本数据插值技术,给定一些随机样本x,PWDE估计由这些样本驱动的PDF。PWDE叠加放置于观察值上核函数,基于这样的方式,每个观察值对于估计PDF都有贡献。利用PWDE来估计PDF的公式如下,P...
  • matlab中提供了核平滑密度估计函数... %红色为参考 PRML上核密度估计(parzen窗密度估计)的实现: %% Demonstrate a non-parametric (parzen) density estimator in 1D %% % This file is from pmtk3.googlecode....
  • 在基于熵的音频相似度度量中,用到Parzen窗法对所提取的MFCC参数进行概率密度函数估计,其MATLAB实现如下:function p=Parzen(xi,x,h1,f)%xi为样本,x为概率密度函数的自变量的取值,%h1为样本数为1时的窗宽,f为窗...

空空如也

空空如也

1 2 3 4 5 ... 15
收藏数 285
精华内容 114
关键字:

parzen window