二项分布和Beta分布

转载 2013年12月06日 00:10:53

http://hyry.dip.jp/tech/slice/slice.html/42


本文通过实例介绍二项分布和Beta分布的含义,并使用pymc对抛硬币进行模拟实验,从而获得Beta分布。

In [15]:
%pylab inline
import pylab as pl
import numpy as np
from scipy import stats
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.kernel.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.

在概率论和统计学中,二项分布是n个独立的[是/非]试验中成功的次数的离散概率分布,其中每次试验的成功概率为p。举两个例子就很容易理解二项分布的含义了:

  • 抛一次硬币出现正面的概率是0.5(p),抛10(n)次硬币,出现k次正面的概率。

  • 掷一次骰子出现六点的概率是1/6,投掷6次骰子出现k次六点的概率。

在上面的两个例子中,每次抛硬币或者掷骰子都和上次的结果无关,所以每次实验都是独立的。二项分布是一个离散分布,k的取值范围为从0到n,只有n+1种可能的结果。

scipy.stats.binom为二项分布,下面用它计算抛十次硬币,出现k次正面的概率分布。

In [16]:
n = 10
k = np.arange(n+1)
pcoin = stats.binom.pmf(k, n, 0.5)
pcoin
Out[16]:
array([ 0.00097656,  0.00976563,  0.04394531,  0.1171875 ,  0.20507813,
        0.24609375,  0.20507813,  0.1171875 ,  0.04394531,  0.00976563,
        0.00097656])
In [17]:
pl.stem(k, pcoin, basefmt="k-")
pl.margins(0.1)

下面是投掷6次骰子,出现6点的概率分布。

In [18]:
n = 6
k = np.arange(n+1)
pdice = stats.binom.pmf(k, n, 1.0/6)
pdice
Out[18]:
array([  3.34897977e-01,   4.01877572e-01,   2.00938786e-01,
         5.35836763e-02,   8.03755144e-03,   6.43004115e-04,
         2.14334705e-05])
In [19]:
pl.stem(k, pdice, basefmt="k-")
pl.margins(0.1)

对于硬币或者骰子这样的简单实验,我们事先能很准确地掌握系统成功的概率。然而通常情况下,系统成功的概率是未知的。为了测试系统的成功概率p,我们做n次试验,统计成功的次数k,于是很直观地就可以计算出p=k/n。然而由于系统成功的概率是未知的,这个公式计算出的p只是系统成功概率的最佳估计。也就是说实际上p也可能为其它的值,只是为其它的值的概率较小。

例如有某种特殊的硬币,我们事先完全无法确定它出现正面的概率。然后抛10次硬币,出现5次正面,于是我们认为硬币出现正面的概率最可能是0.5。但是即使硬币出现正面的概率为0.4,也会出现抛10次出现5次正面的情况。因此我们并不能完全确定硬币出现正面的概率就是0.5,所以p也是一个随机变量,它符合Beta分布。

Beta分布是一个连续分布,由于它描述概率p的分布,因此其取值范围为0到1。 Beta分布有αβ两个参数,其中α为成功次数加1,β为失败次数加1。

连续分布用概率密度函数描述,下面绘制实验10次,成功4次和5次时,系统成功概率p的分布情况。可以看出k=5时,曲线的峰值在p=0.5处,而k=4时,曲线的峰值在p=0.4处。

In [20]:
n = 10
k = 5
p = np.linspace(0, 1, 100)
pbeta = stats.beta.pdf(p, k+1, n-k+1)
plot(p, pbeta, label="k=5", lw=2)

k = 4
pbeta = stats.beta.pdf(p, k+1, n-k+1)
plot(p, pbeta, label="k=4", lw=2)
xlabel("$p$")
legend(loc="best");

下面绘制n=10,k=4n=20,k=8的概率分布。可以看出峰值都在p=0.4处,但是n=20的山峰更陡峭。也就是说随着实验次数的增加,p取其它值的可能就越小,对p的估计就更有信心,因此山峰也就更陡峭了。

In [30]:
n = 10
k = 4
p = np.linspace(0, 1, 100)
pbeta = stats.beta.pdf(p, k+1, n-k+1)
plot(p, pbeta, label="n=10", lw=2)

n = 20
k = 8
pbeta = stats.beta.pdf(p, k+1, n-k+1)
plot(p, pbeta, label="n=20", lw=2)
xlabel("$p$")
legend(loc="best");

假设我们的知识库中没有Beta分布,如何通过模拟实验找出p的概率分布呢?pymc是一个用于统计估计的库,它可以通过 先验概率和 观测值 模拟出 后验概率 的分布。下面先解释一下这两个概率:

  • 先验概率:在贝叶斯统计中,某一不确定量p的先验概率分布是在考虑"观测数据"前,能表达p不确定性的概率分布。

  • 后验概率:在考虑相关证据或数据后所得到的不确定量p的概率分布。

拿前面抛硬币的实验来说,如果在做实验之前能确信硬币出现正面的概率大概在0.5附近的话,那么它的先验概率就是一个以0.5为中心的山峰波形。而如果是某种特殊的硬币,我们对其出现正面的概率完全不了解,那么它的先验概率就是一个从0到1的平均分布。为了估计这个特殊硬币出现正面的概率,我们做了20次实验,其中出现了8次正面。通过这个实验,硬币出现正面的可能性的后验概率就如上图中的绿色曲线所示。

pymc库可以通过先验概率和观测值模拟出后验概率的分布,这对于一些复杂的系统的估计是很有用的。下面我们看看如何用pymc来对这个特殊硬币出现正面的可能性进行估计:

  • 首先pcoin是这个特殊硬币出现正面的概率,由于我们没有任何先验知识,因此它的先验概率是一个从0到1的平均分布(Uniform)。

  • 假设我们做了20次实验,其中8次为正面。根据前面的介绍可知,出现正面的次数符合二项分布(Binomial),并且这个二项分布的概率ppcoin。这个通过value参数指定了实验的结果。因此experiment虽然是一个二项分布,但是它已经不能取其它值了。

In [32]:
import pymc
pcoin = pymc.Uniform("pcoin", 0, 1)
experiment = pymc.Binomial("experiment", 20, pcoin, value=8)

接下来通过MCMC对象模拟pcoin的后验概率。MCMC是Markov chain Monte Carlo(马尔科夫蒙特卡洛)的缩写,它是一种用马尔可夫链从随机分布取样的算法。通过调用MCMC对象的sample(),可以对pcoin的后验概率分布进行取样。这里30000为取样次数,5000表示不保存头5000次取样值。这时因为MCMC算法通常有一个收敛过程,我们希望只保留收敛之后的取样值。

In [33]:
mc = pymc.MCMC([pcoin])
mc.sample(30000, 5000)
 
[****************100%******************]  30000 of 30000 complete

通过MCMC对象trace()可以获得某个不确定量的取样值。下面的程序获得pcoin的25000次取样值,并用hist()显示其分布情况。由结果可知pcoin的分布与前面介绍的Beta分布一致。

In [31]:
pcoin_trace = mc.trace("pcoin")[:]
hist(pcoin_trace, normed=True, bins=30);
plot(p, pbeta, "r", label="n=20", lw=2)
Out[31]:
[<matplotlib.lines.Line2D at 0x5182190>]
In [34]:
pcoin_trace.shape
Out[34]:
(25000,)

相关文章推荐

伯努利分布、二项分布、多项分布、Beta分布、Dirichlet分布

1. 伯努利分布 伯努利分布(Bernoulli distribution)又名两点分布或0-1分布,介绍伯努利分布前首先需要引入伯努利试验(Bernoulli trial)。 伯努利试验...

带你理解beta分布

相信大家学过统计学的都对 正态分布 二项分布 均匀分布 等等很熟悉了,但是却鲜少有人去介绍beta分布的。用一句话来说,beta分布可以看作一个概率的概率分布,当你不知道一个东西的具体概率是多少时,它...

Beta分布从入门到精通

最近一直有点小忙,但是不知道在瞎忙什么,终于有时间把Beta分布的整理弄完。 下面的内容,夹杂着英文和中文,呵呵~ Beta Distribution Beta Distribution...

概率分布之Beta分布与Dirichlet分布

Beta分布与Dirichlet分布的定义域均为[0,1],在实际使用中,通常将两者作为概率的分布,Beta分布描述的是单变量分布,Dirichlet分布描述的是多变量分布,因此,Beta分布可作为二...
  • jteng
  • jteng
  • 2017年03月04日 21:53
  • 1193

Beta 分布的应用

从随机变量到顺序统计量考虑如下的游戏:有一个魔盒(随机数生成器),上有一个按钮,每按一下按钮,就均匀地输出一个 U∼[0,1]U\sim[0,1]之间的随机数,现在按上下,得到10个随机数,第7大的数...

直观理解Beta分布

像正态分布、二项分布和均匀分布这样的一些分布,在统计学习当中往往会结合一些现实世界中的实际应用来解释,因此对于统计学的初学者来说也很容易清晰地理解这些分布。但是我发现Beta分布就很少会用这样可以凭直...

各种分布的python例子

#-*- encoding:utf-8 -*- import numpy as np from scipy import stats import matplotlib.pyplot as plt #...

在数据处理中常用的5种概率分布的python实现

1.二项分布 from scipy.stats import binom, norm, beta, expon import numpy as np import matplotlib.pyplot...

关于Beta分布、二项分布与Dirichlet分布、多项分布的关系

From :

认识beta和dirichlet分布

认识Beta/Dirichlet分布 本文主要是对rickjin的《LDA数学八卦》以及PRML一书中关于Beta分布和Dirichlet分布的知识整理。 我们知道Gamma函数  ...
  • liuhex
  • liuhex
  • 2017年02月16日 18:00
  • 291
内容举报
返回顶部
收藏助手
不良信息举报
您举报文章:二项分布和Beta分布
举报原因:
原因补充:

(最多只允许输入30个字)