Skip to content

Instantly share code, notes, and snippets.

@seaside2mm
Forked from matsuken92/01_visualize_norm_dist.py
Last active March 19, 2022 01:46
Show Gist options
  • Save seaside2mm/48db46fcc6133bf6cc8d7884221153eb to your computer and use it in GitHub Desktop.
Save seaside2mm/48db46fcc6133bf6cc8d7884221153eb to your computer and use it in GitHub Desktop.
[统计分布可視化] #math
Error in user YAML: (<unknown>): could not find expected ':' while scanning a simple key at line 13 column 1
---
title:      【统计】python再认识最大似然估计
subtitle:   仰望星空,脚踏实地        
date:       2020-05-3          
author:     译者                   
header-img: https://images.unsplash.com/photo-1551288049-bebda4e38f71?ixlib=rb-1.2.1&ixid=eyJhcHBfaWQiOjEyMDd9&auto=format&fit=crop&w=3450&q=80
catalog: true    
mathjax: true                 
tags:                        
    - 统计 
    - python
Catogories: 数学
summary: 当你学习统计学或机器学习时,你会遇到 “似然函数 "的概念。 首先,有人说我看不懂,理解成是 "似是而非 "的 "可能”就足够了。
我想,如果你理解了概率和概率密度函数,就可以用数学的方式来处理这个概率,但我想尝试用图解的方式来做一个稍微直观的理解。

---

原文来源:https://qiita.com/kenmatsu4/items/b28d1b3b3d291d0cc698


当你学习统计学或机器学习时,你会遇到 “似然函数 "的概念。 首先,有人说我看不懂,理解成是 "似是而非 "的 "可能”就足够了。 我想,如果你理解了概率和概率密度函数,就可以用数学的方式来处理这个概率,但我想尝试用图解的方式来做一个稍微直观的理解。

你也可以在jupyter上找到完整的代码。

#以正态分布为例

正态分布的概率密度函数为 $$ {f(x)={1 \over \sqrt{2\pi\sigma^{2}}} \exp \left(-{1 \over 2}{(x-\mu)^2 \over \sigma^2} \right) } $$

在该图中,均值μ和标准差σ两个参数的值都是固定的(在上图中,均值μ=10,标准差σ=3),图中横轴上有$x$作为变量,图中的变量为:$x$是横轴上的变量。 作为输出的纵轴是概率密度$f(x)$。

似然函数的基本概念是:"在对数据进行抽样观察后,数据的原始参数来自于概率分布? " 因此,我认为可以考虑成一个反向的贝叶斯定理。 (事实上,可能性是贝叶斯定理的一个组成部分)。

(这里,"数据 "一词称为 "样本")。

现在考虑一个情况,我们有10个样本$x=(x_1,x_2,⋯,x_{10})$,我们知道它们遵循正态分布,但我们不确定两个参数,即均值μ和标准差σ的值是多少。

我们先来考虑 "同时分布10个样本中的这个值"。 我们还假设这10个样本都是IID(Independent Identical Distributions:从同一分布中独立抽取的样本)。 因为它是独立的,所以可以表示为各概率密度的乘积,并且 $$ {P(x_1, x_2,\cdots,x_{10}) = P(x_1)P(x_2)\cdots P(x_{10}) } $$ 这里,所有的$P(x_i)$被假定为正态分布,所以 $$ {P(x_1, x_2,\cdots,x_{10}) = f(x_1)f(x_2)\cdots f(x_{10}) } $$ 进一步展开, $$ {P(x_1, x_2,\cdots,x_{10}) = \prod_{i=1}^{10} {1 \over \sqrt{2\pi\sigma^{2}}} \exp \left(-{1 \over 2}{(x_i-\mu)^2 \over \sigma^2} \right) } $$ 我们现在有一个10个样本的概率密度函数。 但是等一下。 现在我们把样本作为一个确定的值,不再是一个不确定的概率值。 相反,不了解的是两个参数,即均值$μ$和标准差$σ$。 因此,我们可以把$x_i$视为常数,μ和σ是变量。即从样本值计算模型参数。

函数的形式完全相同, 将变量重新定义为μ和σ,其似然函数(likelihood funcion)定义为 $$ {L(\mu, \sigma) = \prod_{i=1}^{10} {1 \over \sqrt{2\pi\sigma^{2}}} \exp \left(-{1 \over 2}{(x_i-\mu)^2 \over \sigma^2} \right) } $$

图像理解

由于μ和σ是未知的,如果我们假设μ=0和σ=1,并画出一个图形,则

这将是一个看上去好像模型偏离很多。 这个时候,似然函数也很小。

由于概率(密度)是通过概率(密度)乘以样本数来计算的,所以0和1之间的数字必须乘以许多倍,结果是一个相当小的数字,几乎为零)。 经常用对数似然,因为乘法换成加法很容易,这种情况下,看上图题目中的数字就很容易理解了。

由于样本几乎完全位于零,所以$L(μ,σ)$也是一个相当小的概率(对数概率:约为-568)。

这一次,让我们试着让$μ=5$和$σ=4$。

虚线是每个样本对应的概率。 感觉比之前准确一些。 这一次,对数的可能性约为-20,这个数字要大得多。

动画模拟

我们来看看当$μ$变化时,似然函数是如何变化的,因为可以看到$μ=10$时,对数概率达到了最大值

当σ变化时,我们可以看出,当$σ=2.7$时,对数似然达到最大值。 原先生成的数据是$σ=3$,所以有一点误差,但数值相差无几,但也很接近。

num_frame = 30

min_x = -11
max_x = 21

x = np.linspace(min_x, max_x, 201)

def norm_dens(val, m, s):
    return (1/np.sqrt(2*np.pi*s**2))*np.exp(-0.5*(val-m)**2/s**2)

def animate(nframe):
    global num_frame
    plt.clf()

    m = 10
    s = nframe/float(num_frame) * 5
    y = norm_dens(x, m, s)

    L = np.prod([norm_dens(x_i, m, s) for x_i in data])
    l = np.log(L)

    plt.xlim(min_x, 16)
    plt.ylim(-0.01,.6)

    # 正态函数
    plt.plot(x,y)

    # 样本点
    plt.scatter(data, np.zeros_like(data), c="r", s=50)
    for d in data:
        plt.plot([d, d], [0, norm_dens(d, m, s)], "k--", lw=1)

    plt.title("sd:{0:.3f}, Likelihood:{1:.5f}, log Likelihood:{2:.5f}".format(s, L, l))

    #plt.show()

fig = plt.figure(figsize=(10,7))
anim = ani.FuncAnimation(fig, animate, frames=int(num_frame), blit=True)
anim.save('likelihood_s.gif', writer='imagemagick', fps=1, dpi=64)

最大似然估计

图中显示的是μ变化时对数似然的变化。 对似然函数的μ求导为0时,μ的值将是10。 这就是最大似然估计。 (我们暂时假设$\sigma$是固定的)。

这也是$\sigma$的变化概率图,可以看出,最大的概率是$\sigma=3$左右

最后,如果我们同时考察μ和σ,可以看出在μ大于10,σ稍小3的点上,似然值估计最大。

# 等高線
plt.figure(figsize=(8,5))
mu = np.linspace(5, 15, 200)
s = np.linspace(0, 5, 200)
MU, S = np.meshgrid(mu, s)

Z = np.array([(np.prod([norm_dens(x_i, a, b) for x_i in data])) for a, b in zip(MU.flatten(), S.flatten())])
plt.contour(MU, S, Z.reshape(MU.shape), cmap=cm.Blues)
plt.xlabel("mu")
plt.ylabel("s")
title date author header-img mathjax tags Catogories Summary
【统计】python标准偏差再讨论
2020-05-01 09:51:06 -0700
译者
true
统计
python
数学
对于想学习统计学的人来说,有一个非常重要但又很难理解的概念就是 "标准差"。 对于 "平均"大家很熟悉,你以为 "我懂了,我懂了",却突然出现了 "标准差"。

原文来源:https://qiita.com/kenmatsu4/items/e6c6acb289c02609e619

对于想学习统计学的人来说,有一个非常重要但又很难理解的概念就是 "标准差"。 对于 "平均"大家很熟悉,你以为 "我懂了,我懂了",却突然出现了 "标准差"。 $$ {\sigma = \sqrt{ {1 \over n} \sum_{i=1}^n(x_i - \bar{x})^2}

} $$ 有的人可能会不知所措,觉得不会做。

如果看图的图像,下面的红线的长度就是 "标准差"。 我们也会发现,为什么这个长度是标准差。

代码06 这边

在这篇文章中,我将试着为不擅长数学的人讲解一下什么是标准差。 如果你理解了公式,但不知道 "标准差 "是什么意思,我们将用一种方式来解释,帮助你直观地理解它。

在这篇文章中,n作为标准差的分母,在某些情况下,根据所要分析的情况,也可以用n-1来表示,但在这里,我们将用比较简单的符号n来表示。

1.平均

我们再来考虑一次 "平均"。 平均数有不同类型,如 "算术平均数"、"几何平均数"、"谐波平均数 "等,但所谓的 "平均数 "就是大家熟悉的 "算术平均数"。

所以你把所有的数据加起来,再除以这个数字。 在统计学中,这个 "均值 “被称为$\bar x$,定义如下 数据的数量假设为$n$。 $$ {\bar{x} = {1 \over n} \sum_{i=1}^{n} x_i } $$ 用图表表示如下,这是很直观的。

代码02 这边

2. 什么是 "偏差 "?

现在,下一步是对 "偏差 "概念的解释。Deviation, 偏差,是指每个数据与平均值之间的差值,如下图所示。 红线是每个数据的偏差。

代码03 这边

3. 平均偏差

现在,在说到标准差之前,我先介绍一下 "平均偏差",这很容易直观地理解。这是前文中介绍的 "偏差 "的平均数。 换句话说

ID 偏差
1 15
2 -18
3 4
4 -15
5 10
6 8
7 -4

但只有一个问题:如果你把它们全部加起来,就会得到0。 我想做的是上图红线长度的平均值,即

因此,考虑将长度的负的部分倒转,即绝对值之后再平均。

取这些 "偏差的绝对值 "的平均值,结果是10.57,也就是所谓的 "平均偏差"。衡量你离平均水平有多远的一个指标。本例中,与平均分81分相差10.57分左右。可以毫不夸张地说,这个概念几乎是标准偏差的想法。 只是计算方法稍有不同而已。所以,当你说 "标准差 "的时候,你可以把它想成是 =="它离平均数有多远的平均数"==。

公式表示: $$ \frac{1}{n}\sum_{i=1}^n|\bar x - x| $$ 写成图形时,它的表示方法是这样的,斜率反转,使数值在x为负值的区域内为正值。

4. 标准偏差

现在,终于,这篇文章的主角--"标准偏差"。 在 “平均偏差 "中,我们用绝对值将负数变为正数,而在 "标准偏差 "中,我们用平方。换句话说,思路是完全一样的,只是去掉负号的方式不同。

从下面可以看到,原点是均值点,相似之处是将负值改为正值。

在上一节中,与均值的差值用下图所示的线段长度来表示,

这次用平方面积来评价与均值的差异,如图所示。

把这些区域加在一起,取平均数。

公式表达为, $$ {\sigma = \sqrt{ {1 \over n} \sum_{i=1}^n(x_i - \bar{x})^2}} $$

5. 与方差联系

方差 说明数据的分散程度,可看做以面积为单位。 这时可以取根号,将面积变成它的长度。另外这也形象表示方差的sqaure的含义。

x = [96, 63, 85, 66, 91, 89, 77]
ave = np.average(x)
total = 0
for i in range(len(x)):
    total += (x[i] - ave)**2
print np.sqrt(total/len(x))

另外,改成 "平方 "而不是 "绝对值”的原因,这有点难以解释,但有一点是,用于平均偏差的绝对值在数学上很难处理,所以我们采用了容易处理的平方。 另外,我认为其中一个原因是,其中包含了几何中距离的概率。

一言以蔽之

  • 数学上的处理方式很方便。
  • 距离的概念有平方和。

6.例

##1:用图形理解

就像我在这篇文章的开头发的图。 这个数据的均值用红圈表示,逐一计算出每个数据的均值的偏差,再进一步计算出标准差,用红条的长度来表示。

所有数据的平均距离的平均值就是这个红条的长度。 标准差条的长度也会增加,因为数据是在一个大范围内发生的,所以标准差条的长度也会增加。

2:偏差值

【注】日本高考只有偏差值,没有分数。

偏差值也是以 "标准差 "为基础。计算每个人参加考试的平均值和标准偏差。

  • 平均分的人偏离度为50。
  • 得分比平均分高一个标准差的人的偏差为60。
  • 得分低于平均分一个标准差的人,其偏差为40。

按以下方法计算: $$ {偏差値 = { (分数 - 平均) \over 标准偏差} \times 10 + 50 } $$ 该值的计算方法为如果你的分数高于上图中的红条长度的平均分,说明你的偏离度为60!

最后

你觉得怎么样?通过对标准差的理解,是不是统计学很有意思!如果越来越多的人对数据分析感兴趣,我会很高兴的!

如果有不明白的地方,请在评论区说出来!

我还写过一篇文章(幻灯片),叫《统计学基础》,如果你喜欢,请你把它作为相关文章来参考。 http://qiita.com/kenmatsu4/items/5a59a7375140f29b31c2

title subtitle date author header-img catalog mathjax tags Catogories
【统计】python的中心极限定理再认识
仰望星空,脚踏实地
2020-04-29
译者
true
true
统计
python
数学

原文来源:https://qiita.com/kenmatsu4/items/351284ef430bcfd2c8ed


1.中心極限定理

如果你研究统计学,你一定会看到一个比较生硬的名字--中心极限定理。

根据维基百科的说法。

中心极限定理概率论中的一组定理。中心极限定理说明,在适当的条件下,大量相互独立随机变量的均值经适当标准化后依分布收敛正态分布。这组定理是数理统计学和误差分析的理论基础,指出了大量随机变量之和近似服从正态分布的条件。

但我不太明白...... 无论原始分布是什么,从它中抽取的样本均值都会接近正态分布。 样本方差也将接近于正态分布。 (准确地说,如果根据卡方分布,如果N大,可以用正态分布来近似)。

我想即使我们用语言解释,或者用公式证明(乘积率参数函数匹配之类的东西),也无法直观的理解它,所以本文的目的是想通过作图来尝试理解它。

更多理论内容,参考wiki

2. 准备作图

我们将用Python编写一个图形,准备过程如下。 本节提供了导入各种库和绘制图形的功能。

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import numpy.random as rd
import matplotlib.mlab as mlab
import scipy.stats as st

# 参数
n = 10000
sample_size = 10000

# 计算样本均值方差的函数
def sample_to_mean_var(sample):
    mean = np.mean(sample)
    var  = np.var(sample)
    return [mean, var]

# 画图
def plot_mean_var(stats, dist_name=""):
    mu = stats[:,0]
    var = stats[:,1]
    bins = 40

    # 均值图
    plt.figure(figsize=(7,5))
    plt.hist(mu, bins=bins, density=True, color="plum")
    plt.title("mu from %s distribution"%(dist_name))
    plt.show()

    # 方差图
    plt.figure(figsize=(7,5))
    plt.hist(var, bins=bins, color="lightblue", density=True)
    plt.title("var from %s distribution"%(dist_name))
    plt.show()

def plot_dist(data, bins, title =""):
    plt.figure(figsize=(7,5))
    plt.title(title)
    plt.hist(data, bins, color="lightgreen", density=True)
    plt.show()

3. 绘图

指数分布

首先,让我们尝试一下指数分布。 下图是参数λ为0.1、10000个样本的指数分布图。 它是一个完全不对称的分布,右边是长长的下摆。

lam = 0.1  
x = rd.exponential(1./lam, size=sample_size)
plot_dist(x, 100, "exponential dist")

下载

以这组10000个样本为一组,由这里计算出样本均值和样本方差。 在重复10,000次后,写出样本平均数和样本方差的直方图,我们得到如下结果.

# 由指数分布生成的样本
lam = 0.1
stats = np.array([sample_to_mean_var(rd.exponential(1./lam, size=sample_size)) for i in range(n)])
plot_mean_var(stats, dist_name="exponential")

我不知道,原来的分布相当偏斜,但样本均值和样本方差似乎是一个美丽的对称钟形。 中心极限定理是,样本均值和方差遵循正态分布。

下面,我也用其他的分布图来试试。

卡方分布

# 自由度5
df = 5
x = rd.chisquare(df, sample_size)
plot_dist(x, 50, "chi square dist")

df = 5   # 自由度

chi_stats = np.array([sample_to_mean_var(rd.chisquare(df, sample_size)) for i in range(n)])
plot_mean_var(chi_stats, dist_name="chi square")

双峰正态分布

我也会尝试一下奇形怪状的分布,比如两座山。

# 双峰正規分布
def generate_bimodal_norm():
    x = np.random.normal(0, 4, sample_size)
    y = np.random.normal(25, 8, sample_size)
    return np.append(x,y)

z = generate_bimodal_norm()
plot_dist(z, 70, "bi-modal normal dist")

binorm_stats = np.array([sample_to_mean_var(generate_bimodal_norm()) for i in range(n)])
plot_mean_var(binorm_stats, dist_name="bi-modal normal")

即使是这样的分布,样本均值或样本方差也是正态分布。 すごい,中心极限定理。

#4.结论 于是,我试着通过看图来直观地理解中心极限定理,虽然看公式和证明,似乎很难。 这也是正态分布在统计学中如此重要的原因:smile:

原文来源:https://qiita.com/kenmatsu4/items/74b8a9f696507af410a4


在统计学中,正态分布起着非常重要的作用。 表示正态分布(密度函数)的公式为 $$ {\Phi(x) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp \left( -\frac{(x-\mu)^2}{2\sigma^2} \right) } $$ 但这是一个非常复杂的公式......当方差$σ^2=1$,均值$μ=0$的标准正态分布时,看起来就简单多了。

$$ {\phi(x) = \frac{1}{\sqrt{2\pi}}\exp \left( -\frac{x^2}{2} \right) } $$ python作图则是如下钟型:

正态分布是平滑的、对称的。正态分布是用一个函数来表示连续概率,可以直观表现聚集在一点上的概率。 于是,用python看看二次函数的图像, $$ {f(x) = -x^2 } $$

越来越像了。 把它放在e的上面,做成一个钟的形状,python表示。

形状已经完全是正态分布。 这个正态分布的原始形式是$e^{-x^2}$。

此后,为了便于计算微分时的计算,$x$是$1/\sqrt{2}$。 换句话说,转换变量为$y = \sqrt{2}x$。 $$ {g(y) = \exp \left(-\frac{y^2}{2} \right) } $$ 对比图如下

它稍微向旁边扩散了一下。

为了让这个$f(x)$的面积为1(因为它是一个概率,所以所有可能发生的事件加起来都是100%),我们需要把积分成1。

由高斯积分(见维基百科)可知,$x$的整个范围内的积分值为 $$ {\int_{-\infty}^{\infty} e^{-x^2}dx = \sqrt{\pi} } $$ 所以,带入上面转换变量, $$ {\int_{-\infty}^{\infty} \exp \left( {-\frac{y^2}{2}}\right)dy = \sqrt{2\pi} } $$ 两边除以$2π$, $$ {\int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi}} \exp \left( {-\frac{y^2}{2}}\right)dy = 1 } $$ 我们得到了标准正态分布的公式:blush:

import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(-3,3, 500)

y1 = np.exp(-(x**2))
y2 = np.exp(-(x**2)/2)
y3 = (1/np.sqrt(2*np.pi))*np.exp(-(x**2)/2)

plt.xlim(-3,3)
plt.ylim(0,1.1)
         
plt.plot(x,y1,"b", label="exp(-x^2)")
plt.plot(x,y2,"g", label="exp(-(x^2)/2")
plt.plot(x,y3,"r", label="normal")
plt.legend()
plt.show()

这个公式是在$e^{-x^2}$的基础上,经过调整,如果把它积分面积就会是1。

title date author header-img catalog mathjax tags Catogories summary
【统计】python认识Q-Q图
2020-05-05 06:30:12 -0700
译者
true
true
统计
python
数学
我将写一篇文章,用一个动画图来解释Q图的意义。 我想试着解释一下这张图,因为这张图有点离奇,需要一些技巧来理解。在R中可以用qqnorm写出Q-Q图,但这是个黑匣子,不知道怎么操作,所以我自己用Python写的。

我将写一篇文章,用一个动画图来解释Q图的意义。 我想试着解释一下这张图,因为这张图有点离奇,需要一些技巧来理解。在R中可以用qqnorm写出Q-Q图,但这是个黑匣子,不知道怎么操作,所以我自己用Python写的。

这次要使用公寓租金数据。

Walk_min distance Price Type Area Direction Year
0 8 B 7900 1K 30.03 3
1 9 B 8500 1K 21.9 5
2 10 B 10800 1K 27.05 4
3 10 B 10800 1K 29.67 4
... ... ... ... ... ... ... ...
185 11 B 8600 1K 20.79 北東 0
186 8 B 7100 1K 22 西 17
187 9 B 18400 1LDK 54.68 西 10

您可以从这里的统计学考试二级官方教材中间的 "资料下载 "链接中下载。 解压下载的压缩文件,解压后在[第2章]-[文本]文件夹中解压出Mansion2.data,这就是我们要用到的数据。

而当你有了数据之后,首先要做的就是绘制图表,建立一个数据的图像。

价格柱状图看出来靠近在左侧,下摆右侧。 另外价格和尺寸之间似乎有一定的关联性。

既然这个Q-Q图将重点关注价格,那么我们将尝试在价格上再进一步解读该图。也就是说,"这个分布是否遵循正态分布?实际上,笔者根据这个数据得到的均值和标准差,应用正态分布的密度函数,如下图所示。

下面是绘制上述一系列图形的Python代码

什么是Q-Q图?

From wiki

是在统计学中,通过比较两个概率分布的分位数对这两个概率分布进行比较的概率图方法。首先选定分位数的对应概率区间集合,在此概率区间上,点(x,y)对应于第一个分布的一个分位数x和第二个分布在和x相同概率区间上相同的分位数。因此画出的是一条含参数的曲线,参数为概率区间的分割数[1]

如果被比较的两个分布比较相似,则其分位图近似地位于y = x上。如果两个分布线性相关,则分位图上的点近似地落在一条直线上,但并不一定是y = x。分位图同样可以用来估计一个分布的位置参数。

是的,这是一个Q-Q图的数据 "价格",也就是租金数据。 一眼看去,很难说清楚这张图显示的是什么。 为了给出一个课本上的解释,Q图是将得到的数据与理论上的分布进行比较,并考察它们的相似性。如果是相似的, 绘出的点在一条直线上。

那么上图应该如何解读呢? 上述图4被认为是对图2的改动。 换句话说,该图直观地显示了正态分布密度函数的理论分布与得到的租金数据是否是一条直线,与所得到的租金数据有多相似。

Q-Q图的结构

现在,这张图应该是类似于理论上的线性分布,但我将解释一下为什么是这样的,因为我认为有必要了解如何绘制这张图。

所以我们再用租金数据来说明一下。 其分布的形状是这样的。

从这里开始,用两个中间的图形来绘制Q-Q图。

首先要做的是将这些租房数据按递减的顺序逐一摆放,并打点,画出这个租房数据的图形。 一共有188个数据,均匀间隔在0到1之间。

作为第二张图,由于我们假设正态分布是理论上的正态分布,所以我们将写出并使用正态累积分布函数图。 这也是,这也代表了累积密度函数,188个点的累积密度函数与租金数据相等。

把这两个图形结合起来,我们可以写出一个Q-Q图的图形。 让我们用动画图来看看。

前两张图为此图右上和左下的中间过程图。 左上方是目标Q-Q图。 首先,右上角的租金数据图的横轴用分位数表示,左下角的正态累计分布函数的纵轴用分位数表示。 在右上角和左下角同时将这个量程从0滑到1。 用一条黑线表示。 黑线相交的点用红点表示。 这些红点同时被绘制出来,这就是Q-Q图。 点线代表的是它。 这个Q-Q图中的 "Q "代表 "分位数",我认为它的名字来源于右上图和左下图中的分位数是同时移动的。

按照正态分布的随机变量的Q-Q图

现在,你说如果数据和理论分布相同,那么Q-Q图就会是一条直线,所以我也想尝试一下。 这意味着使用遵循正态分布的随机变量。 下面是一个由188个随机数组成的直方图,这些随机数遵循正态分布。

为了进一步说明Q-Q图.....这绝对是一条直线,不是吗?

按指数分布的随机变量的Q-Q图

接下来是指数分布。

在这样的形状中,正态Q-Q图的右下角会凸起。

根据F分布的随机数的Q-Q图

F型分布,下摆稍长,向右转。

这在右下角也有一个凸Q-Q图。

根据β分布的随机数的Q-Q图

接下来,我将用左边的长下摆类型的分布和α=6, β=2的β分布写出Q-Q图。

这时,在左上角画出一个凸Q-Q图。

变成是α=0.5和β=0.5的β分布,两边的顶点都是顶点。 在本例中,我们画了一个Q-Q图,下半部分右下角向中间凸,上半部分左上角向中间凸。

本页描述的图形的完整Python代码可以在这里找到.

Error in user YAML: (<unknown>): could not find expected ':' while scanning a simple key at line 11 column 1
---
title: 【统计】核密度估计
date: 2020-05-05 12:14:50      
author:     译者                   
header-img: https://images.unsplash.com/photo-1551288049-bebda4e38f71?ixlib=rb-1.2.1&ixid=eyJhcHBfaWQiOjEyMDd9&auto=format&fit=crop&w=3450&q=80
catalog: true    
mathjax: true                 
tags:                        
    - 统计 
    - R
Catogories: 数学
Summary:对于大量一维数据的可视化,除了使用直方图(Histogram),还有一种更好的方法:核密度估计(Kernel Density Estimates,简称KDE)。
---

有一些数据,想“看看”它长什么样,我们一般会画直方图(Histogram)。现在你也可以用核密度估计。

什么是“核”

如果不了解背景,看到“核密度估计”这个概念基本上就是一脸懵逼。我们先说说这个核 (kernel) 是什么。

首先,“核”在不同的语境下的含义是不同的,例如在模式识别里,它的含义就和这里不同。在“非参数估计”的语境下,“核”是一个函数,用来提供权重。例如高斯函数 (Gaussian) 就是一个常用的核函数。

让我们举个例子,假设我们现在想买房,钱不够要找亲戚朋友借,我们用一个数组来表示 5 个亲戚的财产状况: [8, 2, 5, 6, 4]。我们是中间这个数 5。“核”可以类比 成朋友圈,但不同的亲戚朋友亲疏有别,在借钱的时候,关系好的朋友出力多,关系不好的朋友出力少,于是我们可以用权重来表示。总共能借到的钱是: 8*0.1 + 2*0.4 + 5 + 6*0.3 + 4*0.2 = 9.2

那么“核”的作用就是用来决定权重,例如高斯函数(即正态分布):

Kernel Exponential

如果还套用上面的例子的话,可以认为在 3 代血亲之外的亲戚就基本不会借钱给你了。

最后呢,一般要求核函数有下面两个性质:

  • 归一化:$\int_{- \infty}^{+ \infty}{K(u) du} = 1$
  • 对称性:对所有 $u$ 要求 $K(−u)=K(u)$

最后的最后: 一些常用的核

核密度估计

理解了“核”,核密度估计就容易理解了。

如果我们画直方图,其实目的是画出“概率密度函数”,而直方图本质上是认为频率等于概率。但这种假设不是必然的。核密度函数就是一种“平滑(smooth)”的手段。相当于是“我说我很牛逼你可能不信,但你可以听听我的朋友们是怎么评价我的,加权平均下就能更好地了解我了”。于是乎:

$(x_1,x_2,…,x_n)$ 是独立同分布的 n 个样本点,它的概率密度函数是 $f$,于是我们的估计: $$ \hat{f_h}(x) = \frac{1}{n}\sum_{i=1}^{n}{K_h(x - x_i)} = \frac{1}{nh}\sum_{i=1}^{n}{K(\frac{x-x_i}{h})} $$ 上面式子中 h 是人为指定的,代表“朋友圈”的大小,正式的叫法是“带宽”(bandwidth) 。而 $x−x_i$ 就是自己与朋友的亲疏程度,当然最后要正归化到 [-1, 1] 之间。

例子

核密度估计器与直方图密切相关,通过使用合适的核,可以被赋予平滑性或连续性等特性。为此,我们用这6个数据点来比较直方图和核密度估计器。

对于核密度估计,我们在每个数据点$x_i$上放置一个标准差为2.25的正态内核(用红色虚线表示)。对这些内核进行求和,得出核密度估计值(蓝色实心曲线)。与直方图的离散性相比,核密度估计的平滑性是显而易见的。对于连续随机变量来说,核密度估计更快地收敛到真实的PDF。

Comparison_of_1D_histogram_and_KDE

选择合适的带宽

选择不同的带宽,核密度估计的结果也大不相同,因此人们研究了一些算法来选择带宽。这方面对理解 KDE 本身没有什么太重要的意义,并且常见的算法在 scipy 里也已经都实现了,这里就不细说了,有兴趣的看看 wiki 吧。

截屏2020-05-05 12.24.31

R实现正态核估计

library(MASS) # Venables & Ripley, (2002), Modern Applied Statistics with S

#设定mu,sigma
set.seed(1)
mu <- c(0,0)
Sigmaa <- matrix(c(1,0,0,1),ncol=2)
Sigmab <- matrix(c(1,.75,.75,1),ncol=2)

#从两个均值一样,方差不同的分布中生成样本数据
dat1 <- mvrnorm(5000, mu = mu, Sigma = Sigmaa)
dat2 <- mvrnorm(5000, mu = mu, Sigma = Sigmab)

# 2维核密度估计
kde1 <- kde2d(dat1[,1],dat1[,2],n=100) 
kde2 <- kde2d(dat2[,1],dat2[,2],n=100)

# Perspective Plots
persp(kde1,phi=45,theta=30,shade=0.3,border=NA)
filled.contour(kde1) 
persp(kde2,phi=45,theta=30,shade=0.3,border=NA)
filled.contour(kde2)

截屏2020-05-05 12.40.18

截屏2020-05-05 12.40.40

截屏2020-05-05 12.40.56

参考

截屏2020-05-05 12.40.28

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment