MATLAB概率分布函数
 
mean Average or mean value of array
 
std Standard deviation
 
corrcoef Correlation coefficients
 
median Median value of array
 
max Largest elements in array
 
min Smallest elements in array
 
统计工具箱函数
 
表Ⅰ-1 概率密度函数
 
函数名对应分布的概率密度函数
 
betapdf 贝塔分布的概率密度函数
 
binopdf 二项分布的概率密度函数
 
chi2pdf 卡方分布的概率密度函数
 
exppdf 指数分布的概率密度函数
 
fpdf f分布的概率密度函数
 
gampdf 伽玛分布的概率密度函数
 
geopdf 几何分布的概率密度函数
 
hygepdf 超几何分布的概率密度函数normpdf 正态(高斯)分布的概率密度函数lognpdf 对数正态分布的概率密度函数
 
nbinpdf 负二项分布的概率密度函数
 
ncfpdf 非中心f分布的概率密度函数
 
nctpdf 非中心t分布的概率密度函数ncx2pdf 非中心卡方分布的概率密度函数poisspdf 泊松分布的概率密度函数
 
raylpdf 雷利分布的概率密度函数
 
tpdf 学生氏t分布的概率密度函数
 
unidpdf 离散均匀分布的概率密度函数
 
unifpdf 连续均匀分布的概率密度函数
 
weibpdf 威布尔分布的概率密度函数
 
表Ⅰ-2 累加分布函数
 
函数名对应分布的累加函数
 
betacdf 贝塔分布的累加函数
 
binocdf 二项分布的累加函数
 
chi2cdf 卡方分布的累加函数
 
expcdf 指数分布的累加函数
 
fcdf f分布的累加函数
 
gamcdf 伽玛分布的累加函数
 
geocdf 几何分布的累加函数
 
hygecdf 超几何分布的累加函数

1 概念
 
如果我们想知道某个随机变量$X$的分布$F$，这在一般情况下当然是无法准确知道的，但如果我们手上有它的一些独立同分布的样本，可不可以利用这些样本？一个很简单的办法就是，把这些样本的“频率”近似为随机变量的“概率”。
 
经验分布函数（empirical distribution function）：给每个点$1/n$的概率质量，得到CDF：
 $${\hat{F}}_n(x)=\frac{{\sum }_{i=1}^{n}I(X_i\le x)}{n}$$
 
2 性质
 
经验分布函数，有什么性质？它可以很好地近似真实的分布函数吗？我们给出如下几个定理。
 
定理：对于任意给定的$x$，有
 
$E({\hat{F}}_n(x))=F(x)$；
$V({\hat{F}}_n(x))=\frac{F(x)(1-F(x))}{n}\to 0$；
$\text{MSE}=\frac{F(x)(1-F(x))}{n}\to 0$；
${\hat{F}}_n(x)\stackrel{P}{⟶}F(x)$。
 
Glivenko-Cantelli定理：$X_1,\dots ,X_n\sim F$，那么
 $$\underset{x}{\sup }|{\hat{F}}_n(x)-F(x)|\stackrel{P}{⟶}0$$
 更准确地说，上式其实是几乎必然收敛的。
 
Dvoretzky-Kiefer-Wolfowitz (DKW) Inequity：$X_1,\dots ,X_n\sim F$，那么$\forall ϵ>0$，有
 $$P\left(\underset{x}{\sup }|{\hat{F}}_n(x)-F(x)|>ϵ\right)\le 2e^{-2n{ϵ}^2}$$
 
利用DKW不等式，可以构造出$F$的非参数的$1-\alpha$置信带：定义$L(x)=\max \left\{{\hat{F}}_n(x)-{ϵ}_n,0\right\}$，$U(x)=\max \left\{{\hat{F}}_n(x)+{ϵ}_n,0\right\}$，其中${ϵ}_n=\sqrt{\frac{1}{2n}\log (\frac{2}{\alpha })}$，那么有
 $$P[L(x)\le F(x)\le U(x),\forall x]\ge 1-\alpha$$
 
3 应用
 
经验分布函数有什么用？它可以用来计算一些statistical functional（统计泛函）。
 
假设要计算的statistical functional为$T(F)$，那么，可以利用经验分布函数，代替未知的分布函数，计算出$\theta =T(F)$的plug-in estimator（嵌入式估计量）：$\hat{\theta }=T({\hat{F}}_n)$。
 
如果存在某个$r(x)$使得$T(F)=\int r(x)dF(x)$，那么$T$就称为linear functional（线性泛函），这是因为这样的$T$必定满足$T(aF+bG)=aT(F)+bT(G)$。对于这样的linear functional $T(F)$，它的plug-in estimator可以写为：
 $$T({\hat{F}}_n)=\int r(x)d{\hat{F}}_n=\frac{1}{n}\sum _{i=1}^{n}r(X_i)$$

 
标准正态分布的分布函数 $\Phi(x)$ 可以说是"数据分析师"统计计算中非常重要的一个函数，基本上有正态分布的地方都或多或少会用上它。在一些特定的问题中，我们"数据分析师"需要大量多次地计算这个函数的取值，比如我经常需要算正态分布与另一个随机变量之和的分布，这时候就需要用到数值积分，而被积函数就包含
 
$\Phi(x)$。如果 $Z\sim N(0,1), X\sim f(x)$，$f$ 是 $X$ 的密度函数，那么 $Z+X$
 
的分布函数就是
 
我们"数据分析师"知道，$\Phi(x)$
 
没有简单的显式表达式，所以它需要用一定的数值方法进行计算。在大部分的科学计算软件中，计算的精度往往是第一位的，因此其算法一般会比较复杂。当这个函数需要被计算成千上万次的时候，速度可能就成为了一个瓶颈。
 
当然有问题就会有对策，一种常见的做法是略微放弃一些精度，以换取更简单的计算。在大部分实际应用中，一个合理的误差大小，例如
 
$10^{-7}$，一般就足够了。在这篇文章中，给大家介绍两种简单的方法，它们都比R中自带的 pnorm() 更快，且误差都控制在
 
$10^{-7}$ 的级别。
 
第一种办法来自于经典参考书 Abramowitz and Stegun: Handbook of
 
Mathematical Functions 的 公式
 
26.2.17 。其基本思想是把 $\Phi(x)$ 表达成正态密度函数 $\phi(x)$
 
和一个有理函数的乘积。这种办法可以保证误差小于 $7.5\times
 
10^{-8}$，一段C++实现可以在 这里 找到。(代码中的常数与书中的略有区别，是因为代码是针对误差函数
 
$\mathrm{erf}(x)$ 编写的，它与 $\Phi(x)$ 相差一些常数)
 
我们来对比一下这种方法与R中 pnorm() 的速度，并验证其精度。
 
library(Rcpp)
 
sourceCpp("test_as26217.cpp") x = seq(-6, 6, by = 1e-6) system.time(y 
 
可以看出，A&S 26.2.17
 
的速度大约是 pnorm() 的三倍，且误差也在预定的范围里，是对计算效率的一次巨大提升。
 
那么还有没有可能更快呢？答案是肯定的，而且你其实已经多次使用过这种方法了。怎么，不相信？看看下面这张图，你就明白了。
 

 
没错，这种更快的方法其实就是两个字：查表。它的基本想法是，我们预先计算好一系列的函数取值
 
$(x_i,\Phi(x_i))$，然后当我们需要计算某个点 $x_0$ 时，就找到离它最近的两个点 $x_k$ 和
 
$x_{k+1}$，再用线性插值的方法得到 $\Phi(x_0)$ 的近似取值：
 
什么？觉得这个方法太简单了？先别急，这里面还有不少学问。之前我们""说了，我们需要保证这种方法的误差不超过
 
$\epsilon=10^{-7}$，因此就需要合理地选择预先计算的点。由于
 
$\Phi(-x)=1-\Phi(x)$，我们暂且只需要考虑 $x$ 为正的情况。如果让 $x_i =
 
ih,i=0,1,\ldots,N$，那么对函数 $f$
 
进行线性插值的误差将不超过( 来源 )
 
其中 $\Vert f’’ \Vert_{\infty}$ 是函数二阶导绝对值的最大值。对于正态分布函数来说，它等于
 
$\phi(1)\approx 0.242$。于是令 $E(x)=10^{-7}$，我们就可以解出 $h\approx
 
0.001818$。最后，只要 $x_N>5.199$，即 $N\ge 2860$ 并另所有 $x>x_N$
 
的取值等于1，就可以保证整个实数域上 $\Phi(x)$ 的近似误差都不超过 $10^{-7}$。
 
这种简单方法的实现我放在了 Github
 
上 ，源程序和测试代码也可以在文章最后找到。下面给出它的表现：
 
library(Rcpp)
 
sourceCpp("test_fastncdf.cpp") x = seq(-6, 6, by = 1e-6) system.time(fasty 
 
与之前的结果相比，相当于速度是 pnorm() 的15倍！
 
我们似乎一直以为，在计算机和统计软件普及以后，一些传统的做法就会慢慢被淘汰，例如现在除了考试，或许大部分的时间我们都是在用软件而不是正态概率表。从教学与实际应用的角度来看，这种做法是 应该进行推广和普及的 ，但这也不妨碍我们从一些“旧知识”中汲取营养。关于这种大巧若拙的做法的故事还有很多，比如广为流传的 这一则 。在计算资源匮乏的年代，数据科学家"数据分析师"们想出了各种巧妙的办法来解决他们遇到的各种问题。现如今计算机的性能已经远不是当年可以媲迹，但前人的很多智慧却依然穿透了时间来为现在的我们提供帮助，不得不说这也是一种缘分吧。http://cda.pinggu.org/view/17447.html

 
图 1 显示的是正态(或高斯)分布。它是一条连续的贝尔曲线，期望两边的值是相等的，可以理解为期望就是平均值。它是一个概率分布，因此曲线下方的面积是1。正态分布是由两个参数完全定义的：期望和标准差，它们是衡量期望两边的值如何分布的一种方式。
 

 
图 1 正态分布
 
期望和标准差分别是用希腊字母 μ 和 σ 来表示的，变量 x 有 n 个样本，这些是由下面的公式定义的：
 

 
因此，期望就是值的和除以值的个数一换句话说，也就是平均值。可以通过值和期望的差值的平方子和除以 n-1，然后对结果开方来得到标准差。对于不同的期望和标准差的值，正态分布的相对宽度和高度分布曲线的变化是相当大的。但是，分布值总是如图 1 所示。这意味着，如果知道一个符合正态分布的变量的期望和标准差，例如在大量人口中个体的身高，就可以知道 95% 的人身高不超过期望的 2σ。标准正态分布的期望为 0，标准差为 1。
 
uniform_distribution 模板定义了可以产生随机浮点值的分布对象类型，默认是 double 类型。默认构造函数创建的是标准正态分布，因此期望是 0，方差是 1.0:
 
std::normal_distribution<> dist; // mu: 0 sigma: 1
 
下面展示了如何创建一个有特定值和标准差的正态分布：
 
double mu {50.0}, sigma {10.0};
 
std::normal_distribution<> norm {mu, sigma};
 
这里定义了一个生成 double 值的分布对象，期望为 50.0，标准差是 10.0。为了生成值，可以将一个随机数生成器传给 norm 函数对象。例如：
 
std::random_device rd;
 
std::default_random_engine rng {rd()};
 
std::cout << "Normally distributed values: "<< norm (rng) << " " << norm (rng) << std::endl; // 39.6153 45.5608
 
可以通过调用对象的成员函数 mean() 和 stddev() 来获取它的期望值和标准差：
 
std::cout<
 
通过调用无参数的成员函数 param()，可以得到一个封装了这两个值的 param_type 对象。为了设置期望或标准差，需要将一个 param_type 对象传给成员函数 Pamm()。分布类有用来获取期望和标准差的成员，param_type 对象拥有和它们的名字相同的成员函数。下 面是一个示例：
 
using Params = std::normal_distribution<>::param_type; // Type alias for readability
 
double mu {50.0}, sigma {10.0};
 
std::normal_distribution<> norm {mu, sigma};// Create distribution
 
auto params = norm.param(); // Get mean and standard deviation
 
norm.param(Params {params.mean(),params.stddev() + 5.0}); // Modify params
 
std::cout << "mu: "<< norm.mean() << " sigma: " << norm.stddev ()<< std::endl; // mu: 50 sigma: 15
 
这里调用无参数的 param() 来获取包含期望和方差的 param_type 对象。在第二个 param() 调用中，通过传入一个 Pams 对象将标准差增加了 5.0。
 
可以通过传入一个 param_type 对象作为一个分布对象调用的第二个参数来临时设置期望和标准差：
 
using Params = std::normal_distribution<>::param_type; // Type alias for readability
 
std::random_device rd;
 
std::default_random_engine rng {rd()};
 
std::normal_distribution<> norm {50.0, 10.0}; // Create distribution
 
Params new_p {100.0, 30.0};// mu=100 sigma=30
 
std::cout << norm(rng, new_p) << std::endl; // Generate value with new_p: 100.925
 
std::cout << norm,mean() << " " << norm.stddev()<< std::endl;// 50 10
 
new_p 定义的期望和标准差只会被应用到它被作为第二个参数传入的 norm 的执行中。原始的期望和标准差会被应用到随后的没有第二个参数的 norm 调用中。
 
成员函数 min() 和 max() 返回的是分布可以产生的最小值和最大值。对于分布来说，这并不是特别有用。因为返回值的类型可以这样表示最大值和最小值：
 
std::cout << "min: " << norm.min () << " max: " << norm.max ()<< std::endl; // min:4.94066e-324 max: 1.7 9769e+308