devbox
知新_RL
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

机器学习中的概率分布_分布问题属于机器学习中什么问题

作者：知新_RL | 2024-04-17 07:15:30

赞

踩

分布问题属于机器学习中什么问题

1. 概述

在机器学习中，采用概率视角来理解数据特征，比较常用的有多种概率分布函数，本文对其作基本的总结。具体而言，主要有均匀分布、高斯分布、伯努利分布、泊松分布、二项分布、多项式分布、贝塔分布、贝塔-二项分布、负二项分布、狄里克雷分布，伽马分布等。

2. 概率分布函数

按随机变量的所属类型即离散型和连续型，概率分布分为离散型随机变量概率分布和连续型变量概率分布。

2.1 离散型概率分布

2.1.1  二项分布(Binomial Distribution)

 想象投掷一枚硬币n次，令变量X\in \{0,1,2,...,n \}表示头朝上出现的次数，如果定义参数p为单次投掷硬币头朝上的概率，则称X服从二项分布，表示为 X\sim Bin(n,p),相应的概率质量函数pmf（连续变量称pdf）定义如下：

                                Bin(k,n|p)=Pr(X=k))\triangleq \binom{n}{k}p^k (1-p)^{n-k}

其中系数表达式定义如下

                                                \binom{n}{k} \triangleq \frac{n!}{(n-k)!k!}

相应的系数就是通常所说的二项式系数，表示从可供选项数目n中选出k个选项的总方案数。二项分布的展示效果及代码如下。

  1.  
  2. import numpy as np
  3. from matplotlib import pyplot as plt
  4.  
  5. import operator as op
  6. from functools import reduce
  7.  
  8. def const(n, r):
  9.     r = min(r, n-r)
  10.     numer = reduce(op.mul, range(n, n-r, -1), 1)
  11.     denom = reduce(op.mul, range(1, r+1), 1)
  12.     return numer / denom
  13.  
  14. def binomial(n, p):
  15.     q = 1 - p
  16.     #y = [const(n, k) * (p ** k) * (q ** (n-k)) for k in range(n)]
  17.     xlist = []
  18.     ylist = []
  19.     for k in range(n):
  20.         x = k
  21.         y = const(n, k) * (p ** k) * (q ** (n-k))
  22.         xlist.append(x)
  23.         ylist.append(y)
  24.     return xlist, ylist, np.mean(ylist), np.std(ylist)
  25.  
  26. for ls in [(0.5, 20), (0.7, 40), (0.5, 40)]:
  27.     p, n_experiment = ls[0], ls[1]
  28.     #x = np.arange(n_experiment)
  29.     x, y, u, s = binomial(n_experiment, p)
  30.     plt.scatter(x, y, label=r'$\mu=%.2f,\ \sigma=%.2f$' % (u, s))
  31.  
  32. plt.legend()
  33. #plt.savefig('./graph/binomial.png')
  34. plt.show()

2.1.2 伯努利分布(Bernoulli Distribution)

想象投掷一枚硬币，得到的结果是正面（头像）朝上和反面朝上，用变量X来定义这个事件，则X\in \{0,1 \},则称这个事件X服从伯努利分布,表示为 X\sim Ber(p),其中参数p为头朝上的概率，相应的概率质量函数pmf（连续变量称pdf）定义如下：

                                Ber(x|p)=p^{\mathbb{I}(x=1)}(1-p)^{\mathbb{I}(x=0)}

 上式可简化表示为：

                                  Ber(x|p)=\begin{cases} p& \text{ if } x=1 \\ 1-p& \text{ if } x=0 \end{cases}

由上可知，伯努利分布是二项分布一个特例，即n=1。伯努利分布的展示效果及代码如下。

  1.  
  2. import random
  3. import numpy as np
  4. from matplotlib import pyplot as plt
  5.  
  6. def bernoulli(p, k):
  7.     return p if k else 1 - p
  8.  
  9. n_experiment = 100
  10. p = 0.6
  11. x = np.arange(n_experiment)
  12. y = []
  13. for _ in range(n_experiment):
  14.     pick = bernoulli(p, k=bool(random.getrandbits(1)))
  15.     y.append(pick)
  16.  
  17. u, s = np.mean(y), np.std(y)
  18. plt.scatter(x, y, label=r'$\mu=%.2f,\ \sigma=%.2f$' % (u, s), marker="o", c= "red")
  19. plt.legend()
  20. #plt.savefig('./graph/bernoulli.png')
  21. plt.show()

2.1.3 多项分布(Multinomial Distribution)

多项分布是二项分布的推广。二项分布做n次伯努利试验，规定了每次试验的结果只有两个比如0,1，常用于模拟双面硬币的投掷结果。如果现在还是做n次试验，但每次试验的结果可以有K个（比如投掷骰子），且K个结果发生的概率互斥且和为1，则发生其中一个结果X次的概率就是多项分布。令\mathbf x=(x_1,x_2,...,x_K)表示随机变量，对应的是上述投掷结果，其中x_j(1\leqslant j\leq K)表示的是第j面朝上的次数，则相应的pmf如下：

                                        Mu( \mathbf x|n,\mathbf p)=\begin{pmatrix} n\\ x_1,x2,...x_K \end{pmatrix}\prod_{j=1}^{K}p_j^{x_j}

其中p_j表示第j面朝上的概率，系数表达式定义如下

                                        \begin{pmatrix} n\\ x_1,x2,...x_K \end{pmatrix}=\frac{n!}{x_1!x_2!\begin{matrix}\cdot\cdot \cdot x_K! \end{matrix}}

上述系数就是通常所说的多项式系数。具体而言，多项式系数的推导可理解为：

一共有 \small \binom{n}{x_1}种方法确定第一个子集，第一个子集确定后，从剩下的\small n-x_1中确定第二个子集的可选方法\small \binom{n-x_1}{x_2},再从剩下的\small \small n-x_1-x_2中确定第三个子集\small \binom{n-x_1-x_2}{x_3}，以此类推，第\small K个子集中有可选方法\small \binom{n-x_1-x_2-...-x_{K-1}}{x_K}，n=\sum_{k=1}^{K}x_k   是投掷的总次数。

                                                                   ​​​​​​​        ​​​​​​​           \frac{n!}{x_1!(n-x_1)!}\cdot \frac{(n-x_1)!}{x_2!(n-x_1-x_2)!}\cdot \frac{(n-x_1-x_2)!}{x_3!(n-x_1-x_2-x_3)!}\cdot \cdot \cdot \frac{(n-x_1-\cdot \cdot \cdot -x_{K-1})!}{x_K!(n-x_1-\cdot \cdot \cdot -x_K)!}

经过消去约简后得到：

                                            \frac{n!}{x_1!x_2!\begin{matrix}\cdot\cdot \cdot x_K! \end{matrix}}

多项分布的展示效果及代码如下。

  1. import numpy as np
  2. from matplotlib import pyplot as plt
  3.  
  4. import operator as op
  5. from functools import reduce
  6.  
  7. def factorial(n):
  8.     return reduce(op.mul, range(1, n + 1), 1)
  9.  
  10. def const(n, a, b, c):
  11.     """
  12.         return n! / a! b! c!, where a+b+c == n
  13.     """
  14.     assert  a + b + c == n
  15.  
  16.     numer = factorial(n)
  17.     denom = factorial(a) * factorial(b) * factorial(c)
  18.     return numer / denom
  19.  
  20. def multinomial(n):
  21.     """
  22.     :param x : list, sum(x) should be `n`
  23.     :param n : number of trial
  24.     :param p: list, sum(p) should be `1`
  25.     """
  26.     # get all a,b,c where a+b+c == n, a<b<c
  27.     ls = []
  28.     for i in range(1, n + 1):
  29.         for j in range(i, n + 1):
  30.             for k in range(j, n + 1):
  31.                 if i + j + k == n:
  32.                     ls.append([i, j, k])
  33.  
  34.     y = [const(n, l[0], l[1], l[2]) for l in ls]
  35.     x = np.arange(len(y))
  36.     return x, y, np.mean(y), np.std(y)
  37.  
  38. for n_experiment in [20, 21, 22]:
  39.     x, y, u, s = multinomial(n_experiment)
  40.     plt.scatter(x, y, label=r'$trial=%d$' % (n_experiment))
  41.  
  42. plt.legend()
  43. #plt.savefig('./graph/multinomial.png')
  44. plt.show()

2.1.4 类别分布(Categorial Distribution)

 基于上述的分析， 随机变量X\in \{0,1,2,...,n \}，如果取n=1,则对于投掷K面的骰子（特别注意：此时可使用指示函数来表示每个面的取值结果\mathbb{I}(x=K)，相当于每个面只有0或1的可能,共K个面，不等于常见的骰子取1-6），得到的\mathbf x就是由0,1组成的向量,特别地，如果投掷得到第k面，对应的就是第k个向量。在这种情况下，我们可以将x想象为有K种状态的分类随机变量，那么\mathbf x就是它的dummy encoding。也就是

                                                        \mathbf x=[\mathbb{I}(x=1),...,\mathbb{I}(x=K) ]

如果K=3，那么我们就将1,2,3三种状态编码为（1，0，0），（0，1，0），（0，0，1），这就是现在所谓的one-hot encoding，此时可以理解为K面只有一面“hot”,对应的pmf可简化为：

                                                        Mu( \mathbf x|1, \mathbf p)=\prod_{j=1}^{K}{p_j^{ \mathbb{I}(x_j=1)}}

上述的类别分布常用于对不同的类别进行编码，Gustavo Lacerda也建议称之为multinoulli  distribution 与 multinomial distribution进行对应。上述四种最基本的离散分布对比：

        

类别分布的展示效果及代码实现如下。

  1. import random
  2. import numpy as np
  3. from matplotlib import pyplot as plt
  4.  
  5. def categorical(p, k):
  6.     return p[k]
  7.  
  8. n_experiment = 100
  9. p = [0.3, 0.1, 0.6]
  10. x = np.arange(n_experiment)
  11. y = []
  12. for _ in range(n_experiment):
  13.     pick = categorical(p, k=random.randint(0, len(p) - 1))
  14.     y.append(pick)
  15.  
  16. u, s = np.mean(y), np.std(y)
  17. plt.scatter(x, y, c = "red", label=r'$\mu=%.2f,\ \sigma=%.2f$' % (u, s))
  18. plt.legend()
  19. #plt.savefig('./graph/categorical.png')
  20. plt.show()

2.1.5 负二项分布(Negavite Binomial Distribution)

负二项分布的定义不唯一，我们以上述伯努利试验不断重复多次为例，定义在单次伯努利试验中的成功概率为p,也就是前述的\theta,失败的概率是(1 - p )。我们观察这个序列，直到发生预定义的成功次数r, 那么我们看到的失败次数X将具有负二项（或Pascal）分布，记为：X\sim NB(r, p)。
我们可以类似地使用负二项分布来模拟某台机器在它发生故障之前的工作天数（r = 1）。对应的pmf如下。

                                f(k; r, p) \equiv \Pr(X = k) = \binom{k+r-1}{r-1} (1-p)^kp^r

其中r是成功的次数，k是失败的次数，p是成功的概率。在这个公式中，均值与方差分别如下：

                                                        mean =(1-p)r/{p}

                                                        var =(1-p)r/{p}^2

表达式右侧括号中的项就是二项式系数，并且定义如下

                \binom{k+r-1}{r-1} = \frac{(k+r-1)!}{(r-1)!\,(k)!} = \frac{(k+r-1)(k+r-2)\dotsm(r)}{(k)!}

从k+r-1 个样本而不是k+r 个样本中选择了k个失败，因为根据定义， k+r 个样本中的最后一个试验是成功的。上述系数展开式最右边一项继续展开为：

                                (k+r−1)⋯(r)(k)!=(−1)k(−r)(−r−1)(−r−2)⋯(−r−k+1)(k)!=(−1)k(−rk)

这也就是所谓的“负二项”的来历。

通过最后一个表达式和二项式系数，对于每个0\leq p< 1和q=1-p。

               p^{-r} = (1-q)^{-r} = \sum_{k=0}^\infty \binom{-r}{k}(-q)^k = \sum_{k=0}^\infty \binom{-r}{k}(-1)^k q^k = \sum_{k=0}^\infty \binom{k+r-1}{k}q^k

因此概率质量函数的项满足求和为1。

                                                \sum_{k=0}^\infty \binom{k+r-1}{k}(1-p)^kp^r = p^{-r}p^r = 1

​要理解概率质量函数的上述定义，需要注意每个特定序列的 r 次成功和k 次失败的概率是p^r(1-p)^k，因为k  +  r次试验的结果是独立发生的。由于第r次成功总是最后出现，所以仍然需要 从剩余的k+r-1  次试验中选择失败的 k次试验。上述二项式系数，给出了k+r-1 的精确解释。负二项分布的展示效果及代码如下。

  1. import numpy as np
  2. from matplotlib import pyplot as plt
  3.  
  4. import operator as op
  5. from functools import reduce
  6.  
  7. def factorial(n):
  8.     return reduce(op.mul, range(1, n + 1), 1)
  9.  
  10.  
  11. def const(k, r):
  12.     numer = factorial(k+r-1)
  13.     denom = factorial(r-1) * factorial(k)
  14.     return numer / denom
  15.  
  16.  
  17. def negative_binomial(ns, rs, p):
  18.     q = 1 - p
  19.     xlist = []
  20.     ylist = []
  21.     ks = ns - rs # total failed times 
  22.     for x in range(ks):
  23.         r = ns - x # success times r
  24.         y = const(x, r) * (q ** x) * (p ** r) 
  25.         xlist.append(x)
  26.         ylist.append(y)
  27.     return xlist, ylist, np.mean(ylist), np.std(ylist)
  28.  
  29.  
  30. def negative_binomial_total(ns, rs, p):
  31.     q = 1 - p
  32.     xlist = list(range(rs,ns))
  33.     ylist = []
  34.     
  35.     for x in xlist:
  36.         k = x - rs # failed times k
  37.         y = const(k, rs) * (q ** k) * (p ** rs) 
  38.         ylist.append(y)
  39.     return xlist, ylist, np.mean(ylist), np.std(ylist)
  40.  
  41. fig,axes = plt.subplots(2,1, figsize=(11,7)) #get 1x1 subplots
  42. axes[0].set_xlabel("failed times(k)", fontsize="12") #set the coordinate axis-x label
  43. axes[0].set_ylabel("Probability", fontsize="12") #set the coordinate axis-y label
  44. axes[0].set_title(u"Negative Binomial Distribution",fontsize=12)
  45.  
  46. axes[1].set_xlabel("total times(n)", fontsize="12") #set the coordinate axis-x label
  47. axes[1].set_ylabel("Probability", fontsize="12") #set the coordinate axis-y label
  48. #axes[1].set_title(u"Negative Binomial Distribution",fontsize=12)
  49.  
  50. for ls in [(0.1, 3, 30), (0.5, 14, 30), (0.8, 23, 30), (0.3, 9, 30)]:
  51.     p, n_sucess, n_experiment = ls[0], ls[1], ls[2]
  52.     x, y, u, s = negative_binomial(n_experiment, n_sucess, p)
  53.     xt, yt, ut, st = negative_binomial_total(n_experiment, n_sucess, p)
  54.     axes[0].scatter(x, y, label=r'$\mu=%.2f,\ \sigma=%.2f$' % (u, s))
  55.     axes[1].scatter(xt, yt, label=r'$\mu=%.2f,\ \sigma=%.2f$' % (ut, st))
  56.     
  57.  
  58. axes[0].legend()
  59. axes[1].legend()
  60. #plt.savefig('./graph/binomial.png')
  61. plt.show()

2.1.6 泊松分布(Poisson Distribution)

泊松分布由法国数学家 Siméon Denis Poisson 于1837年于他的概率论著作中首次发表。泊松分布适合于描述单位时间内随机事件发生的次数的概率分布。如某一服务设施在一定时间内受到的服务请求的次数，电话交换机接到呼叫的次数、汽车站台的候客人数、机器出现的故障数、自然灾害发生的次数、DNA序列的变异数、放射性原子核的衰变数、激光的光子数分布等等。泊松分布通常用于对小概率事件进行建模。如果令X\in \{0,1,2,...\},有一个泊松分布，那么定义相应的分布参数\lambda>0,对应的表示为X\sim Poi (\lambda)，得到的pmf如下

                                        Poi (x|\lambda) \equiv \Pr(X = k) = e^{-\lambda}\frac{\lambda^k}{k!}

直观理解可参阅 知乎--马同学-如何直观的理解泊松分布。泊松分布的展示效果及代码如下。

  • 当 k 小于 λ 时，k 每增加 1 累积概率增速会逐渐加大
  • 而当 k 大于 λ 后，累积概率的降速也会逐渐加大。

 

  1. import numpy as np
  2. from matplotlib import pyplot as plt
  3.  
  4. import operator as op
  5. from functools import reduce
  6.  
  7. def poisson(x, lamb):
  8.     xlist = []
  9.     ylist = []
  10.     for k in x:
  11.         xlist.append(k)
  12.         numer = np.exp(-lamb)*(lamb**k)
  13.         denom = reduce(op.mul, range(1, int(k+1)), 1)
  14.         y = numer/denom
  15.         ylist.append(y)
  16.     return xlist, ylist, np.mean(ylist), np.std(ylist)
  17.  
  18. for lamb in [0.5, 1, 1.5, 2, 5, 8]:
  19.  
  20.     x = np.arange(0, 20, 1, dtype=np.int32)
  21.     x, y, u, s = poisson(x, lamb=lamb)
  22.     plt.plot(x, y, '-o', label=r'$\mu=%.2f,\ \sigma=%.2f,'
  23.                          r'\ \lambda=%.2f$' % (u, s, lamb))
  24. plt.legend()
  25. #plt.savefig('./graph/exponential.png')
  26. plt.show()

2.1.7 几何分布(Geometric Distribution)

对于上述伯努利试验，获得一次成功所需的伯努利试验次数X的概率分布，在集合上表示为X \in \{1,2,...\}。几何分布给出了第一次成功发生需要k个独立试验的概率，每个试验的成功概率为p，失败的概率为1-p。如果每次试验的成功概率为p，则第k次试验（共k次试验中）第一次成功的概率为:

                                        Geo(x)=Pr(X=k)=(1-p)^{k-1}p

其中 k\geqslant 1, 0< p\leqslant 1。几何分布的展示效果及代码如下。

  1. import numpy as np
  2. from matplotlib import pyplot as plt
  3.  
  4. def geometric(n, p):
  5.     q = 1 - p
  6.     xlist = []
  7.     ylist = []
  8.     for k in range(n):
  9.         x = k
  10.         y = q ** (k-1) * p
  11.         xlist.append(x)
  12.         ylist.append(y)
  13.     return xlist, ylist, np.mean(ylist), np.std(ylist)
  14.  
  15. for ls in [(0.5, 10), (0.7, 15), (0.3, 20)]:
  16.     p, n_experiment = ls[0], ls[1]
  17.     #x = np.arange(n_experiment)
  18.     x, y, u, s = geometric(n_experiment, p)
  19.     plt.plot(x, y, '-o', label=r'$\mu=%.2f,\ \sigma=%.2f$' % (u, s))
  20.  
  21. plt.legend()
  22. #plt.savefig('./graph/binomial.png')
  23. plt.show()

2.1.8 超几何分布(Hypergeometric Distribution)

超几何分布可视为几何分布的一种扩展，描述了从有限N个物件（其中包含M个指定种类的物件）中抽出n个物件，采取不放回方式，成功抽出该指定种类的物件的次数。超几何分布中的参数是M,N,n，超几何分布记作X\sim H(N,M,n) 。假定一个箱子中共有100颗球，70颗为白球，30颗为红球，从中不放回地拿出10个球，那么共取出 k 个红球的概率(pmf)为：

                HyperGeo(x)=Pr(X=k)=\frac{C_M^kC_{N-M}^{n-k}}{C_N^n}

                C_M^k= \frac{M!}{(M-k)!k!},C_{N-M}^{n-k}= \frac{(N-M)!}{(N-M-n+k)!(n-k)!},C_N^n= \frac{N!}{(N-n)!n!}

超几何分布的展示效果及代码实现如下。

  1. import numpy as np
  2. from matplotlib import pyplot as plt
  3.  
  4. import operator as op
  5. from functools import reduce
  6.  
  7. def const(ns, ks):
  8.     ks = min(ks, ns-ks)
  9.     numer = reduce(op.mul, range(ns, ns-ks, -1), 1)
  10.     denom = reduce(op.mul, range(1, ks+1), 1)
  11.     return numer / denom
  12.  
  13. def hyper_geometric(N, M, n):
  14.     xlist = []
  15.     ylist = []
  16.     for k in range(n):
  17.         x = k
  18.         y = const(M, k) * const(N-M, n-k) / const(N, n) 
  19.         xlist.append(x)
  20.         ylist.append(y)
  21.     return xlist, ylist, np.mean(ylist), np.std(ylist)
  22.  
  23. for ls in [(10, 30, 100), (15, 60, 100), (20, 80, 100)]:
  24.     ns, Ms, Ns = ls[0], ls[1], ls[2]
  25.     #x = np.arange(n_experiment)
  26.     x, y, u, s = hyper_geometric(Ns, Ms, ns)
  27.     plt.plot(x, y, '-o', label=r'$\mu=%.2f,\ \sigma=%.2f$' % (u, s))
  28.  
  29. plt.legend()
  30. #plt.savefig('./graph/binomial.png')

2.2 连续型概率分布

2.2.1 均匀分布(Uniform Distribution)

还是以投掷骰子为例，比如结果是1到6。得到任何一个结果的概率是相等的，这就是均匀分布的基础。与伯努利分布不同，均匀分布的所有可能结果各自的个数也是相等的。如果变量X服从均匀分布的，则pdf可表示为：

                                                Unif(x|a,b)=\frac{1}{b-a}\mathbb{I}(a\leq x \leq b)

均匀分布也是常见的用于模拟等概率事件的分布函数，均匀分布效果展示及代码实现如下。

  1. import numpy as np
  2. from matplotlib import pyplot as plt
  3.  
  4. def uniform(x, a, b):
  5.  
  6.     y = [1 / (b - a) if a <= val and val <= b
  7.                     else 0 for val in x]
  8.  
  9.     return x, y, np.mean(y), np.std(y)
  10.  
  11. x = np.arange(-100, 100) # define range of x
  12. for ls in [(-50, 50),(-15, 20), (20, 40), (10, 20)]:
  13.     a, b = ls[0], ls[1]
  14.     x, y, u, s = uniform(x, a, b)
  15.     plt.plot(x, y, label=r'mean=%.2f, variance=%.2f' % (u, s))
  16.  
  17. plt.legend()
  18. #plt.savefig('./graph/uniform.png')
  19. plt.show()

2.2.2 高斯分布(Gaussian Distribution)

根据Kevin P. Murphy的总结，现在广泛提及的高斯分布这一提法存在误导性，并非由高斯首先提出。按照时间先后，高斯分布由数学家Moivre首先发现，而此时Laplace还未出生，然后Laplace 揭示了该分布的一些特殊性质（此时高斯还只有六岁），此后在1800~1810年随着高斯将高斯分布结合天体运动与最小二乘法分析误差的估计得到广泛推广应用，高斯分布（正态分布）从首次出现到最终确立，其时间简史为：

1705年， J. Bernoulli 的著作《猜测的艺术》问世，提出伯努利大数定律；

1730-1733年，棣莫弗Moivre从二项分布的渐进公式得到正态密度函数（高斯分布的雏形）首次提出中心极限定理；

1780年， P. S. M. Laplace建立中心极限定理的一般形成；

1805年，勒让德Legendre 发明最小二乘法；

1809年，高斯引入正态误差理论，不但补充了最小二乘法，而且首次导出高斯分布；

1811年，拉普拉斯利用中心极限定理论证高斯分布；

1837年，海根G.Hagen 提出元误差学说，自此之后，逐步正式确立误差服从高斯分布。

高斯分布的pdf表示如下

                                                \mathcal N(x\mid \mu,\sigma^2) = \frac{1}{ \sqrt{2\pi \sigma^2 }} \exp \left[-\frac{(x - \mu )^2}{2\sigma ^2} \right]

其中\mu,\sigma^2分别表示均值和方差。高斯分布得到广泛使用可总结有以下几个原因：

1) 高斯分布的定义只有两个参数，且这两个参数即均值与方差，抓取了一个概率分布最基本的属性特征。

2)根据中心极限定理，任意一个独立随机变量的和都趋于高斯分布，这使得它能够较好的模拟残差或噪声的分布特征。

3)在满足给定均值与均方差的约束前提下，高斯分布基于的假设最少，不确定性最大则具有最大的熵值。

4)高斯分布的数学形式相对简单，易于积分或微分运算，在其它分布函数中得到广泛的扩展应用。

高斯分布的展示效果及代码如下。

  1. import numpy as np
  2. from matplotlib import pyplot as plt
  3.  
  4. def gaussian(x, n):
  5.     u = x.mean()
  6.     s = x.std()
  7.  
  8.     # divide [x.min(), x.max()] by n
  9.     x = np.linspace(x.min(), x.max(), n)
  10.  
  11.     a = ((x - u) ** 2) / (2 * (s ** 2))
  12.     y = 1 / (s * np.sqrt(2 * np.pi)) * np.exp(-a)
  13.  
  14.     return x, y, x.mean(), x.std()
  15.  
  16. x = np.arange(-150, 150) # define range of x
  17. x, y, u, s = gaussian(x, 10000)
  18.  
  19. plt.plot(x, y, label=r'$\mu=%.2f,\ \sigma=%.2f$' % (u, s))
  20. plt.legend()
  21. #plt.savefig('./graph/gaussian.png')
  22. plt.show()

2.2.3 正态分布(Normal Distribution)

正态分布是高斯分布标准化后的结果，此时均值为0，均方差为1。相应的的pdf如下：

                                                        \mathcal N(x \mid 0,1) = \frac{1}{ \sqrt{2\pi }} \exp \left[-\frac{x ^2}{2} \right]

正态分布的展示效果及代码如下。

  1.  
  2. import numpy as np
  3. from matplotlib import pyplot as plt
  4.  
  5. def normal(x, n):
  6.     u = x.mean()
  7.     s = x.std()
  8.  
  9.     # normalization
  10.     x = (x - u) / s
  11.  
  12.     # divide [x.min(), x.max()] by n
  13.     x = np.linspace(x.min(), x.max(), n)
  14.  
  15.     a = ((x - 0) ** 2) / (2 * (1 ** 2))
  16.     y = 1 / (s * np.sqrt(2 * np.pi)) * np.exp(-a)
  17.  
  18.     return x, y, x.mean(), x.std()
  19.  
  20. x = np.arange(-150, 150) # define range of x
  21. x, y, u, s = normal(x, 10000)
  22.  
  23. plt.plot(x, y, label=r'$\mu=%.2f,\ \sigma=%.2f$' % (u, s))
  24. plt.legend()
  25. #plt.savefig('./graph/normal.png')
  26. plt.show()

2.2.4 狄拉克分布(Dirac Distribution)

对于上述高斯分布当 \sigma^2\rightarrow 0，概率分布中的所有质量都集中在一个点上，使得概率分布图变得无限“瘦长”，此时

                                                        \lim_{\sigma^2\rightarrow 0}\mathcal N(x|\mu,\sigma^2) =\delta(x-\mu)

这里的 \delta 称之为狄拉克-德尔塔函数（Dirac delta function）,也称之为狄拉克分布，也称之为退化概率密度函数（degenerate pdf），亦可定义为

                                                        \delta(x)=\begin{cases} \infty& \text{ if } x=0 \\ 0 &\text{ if } x\neq 0 \end{cases}

积分形式如下：

                                                                \int_{-\infty }^{\infty}\delta(x)dx=1

上述\delta函数的一个重要的属性是具有筛选特性(sifting property)，即能够从一个求和或积分式中筛选出单项：

                                                            \int_{-\infty }^{\infty}f(x)\delta(x-\mu)dx=f(\mu)

  狄拉克-德尔塔函数的展示效果及代码实现如下。  

  1. import numpy as np
  2. import matplotlib.pyplot as plt
  3.  
  4. def dirac_delta(x,sig):
  5.     res = np.zeros_like(x)
  6.     res[(-(1/(2*sig))<=x) & (x<=(1/(2*sig)))] = 1
  7.     return res
  8.  
  9. x = np.linspace(-5,5,1000)
  10.  
  11. for sigma in np.arange(1,5,0.1):
  12.     plt.cla()
  13.     plt.grid()
  14.     plt.title('Dirac Delta function',size=12)
  15.     plt.xlabel('x',size=10)
  16.     plt.ylabel("$\delta(x)$",size=12)
  17.     plt.ylim(0,1)
  18.     plt.plot(x,dirac_delta(x,sigma),color='red')
  19.     plt.pause(0.5)
  20.  
  21. plt.show()

2.2.5 经验分布(Empirical Distribution)

       假设我们有一组  N个样本的数据集 D={x(1),x(2),...,x(N)}，源自于概率分布p(X)，其中X∈R，我们可以使用一组delta函数或“spike”来近似这个pdf，基于这些样本数据可得到：    

                                                p^N(x)=1N∑n=1Nδx(i)(x)

上述这个分布就称为数据集D 的经验分布(empirical distribution), 相应的得到其cdf如下：

                                 F^N(x)=1N∑n=1NI(x(i)≤x)=1N∑n=1Nux(i)(x)​​​​​​​

其中，uy(x) 称为阶跃函数，具体形式如下：

                                                uy(x)={1 if x≥y0 if x<y

经验分布的展示效果及代码实现如下。

        ​​​​​​​

  1. import numpy as np
  2. import statsmodels.api as sm # recommended import according to the docs
  3. import matplotlib.pyplot as plt
  4.  
  5. sample = np.random.uniform(0, 1, 5)
  6. ecdf = sm.distributions.ECDF(sample)
  7.  
  8. x = np.linspace(min(sample), max(sample))
  9. y = ecdf(x)
  10. plt.step(x, y)
  11. plt.show()

 2.2.6 学生分布(Student Distribution)

学生分布由英国统计学家威廉·西利·戈塞特 Gosset 在爱尔兰吉尼斯啤酒厂工作期间提出并发表。为了满足吉尼斯董事关于研究人员发表科学研究不提及啤酒、吉尼斯或人员的姓氏的要求，Gosset 以他的笔名“学生”予以首次发表。学生分布的pdf如下

                                \mathcal \tau (x|\mu,\sigma^2,\nu) = \frac{\Gamma(\frac{\nu+1}{2})} {\sqrt{\nu\pi \sigma^2}\,\Gamma(\frac{\nu}{2})} \left[1+\frac{1}{\nu}\frac{(x - \mu )^2}{\sigma ^2} \right]^{-\frac{(\nu +1)}{2}}

标准化后即\mu=0,\sigma=1,则上述简化为：

                                        f(x) = \frac{\Gamma(\frac{\nu+1}{2})} {\sqrt{\nu\pi}\,\Gamma(\frac{\nu}{2})} \left(1+\frac{x^2}\nu \right)^{-(\nu+1)/2}

高斯分布对数据中的离群点比较敏感，比如具有明显长尾效应的数据点采用高斯分布建模，则得到的分布易受到噪声影响，主要原因在于log-probability 以中心位置仅仅只是二次衰减。而学生分布可视为高斯分布的一个扩展，在远离中心位置快速衰减(阶次高于2)，几乎不受离群点的影响。学生分布的展示效果及代码实现如下。

  1. import numpy as np
  2. from matplotlib import pyplot as plt
  3.  
  4. def gamma_function(n):
  5.     cal = 1
  6.     for i in range(2, n):
  7.         cal *= i
  8.     return cal
  9.  
  10. def student_t(x, freedom, n):
  11.  
  12.     # divide [x.min(), x.max()] by n
  13.     x = np.linspace(x.min(), x.max(), n)
  14.  
  15.     c = gamma_function((freedom + 1) // 2) \
  16.         / np.sqrt(freedom * np.pi) * gamma_function(freedom // 2)
  17.     y = c * (1 + x**2 / freedom) ** (-((freedom + 1) / 2))
  18.  
  19.     return x, y, np.mean(y), np.std(y)
  20.  
  21. for freedom in [1, 2, 3, 5]:
  22.  
  23.     x = np.arange(-10, 10) # define range of x
  24.     x, y, _, _ = student_t(x, freedom=freedom, n=10000)
  25.     plt.plot(x, y, label=r'$v=%d$' % (freedom))
  26.  
  27. plt.legend()
  28. plt.savefig('./graph/student_t.png')
  29. plt.show()

2.2.7 柯西分布(Cauchy Distribution)

柯西分布又称为洛伦兹分布( Cauchy distribution/Lorentiz  distribution)，是学生分布的一个特例，即相应的自由度参数\nu =1。

                                        Cau (x|\mu,\sigma^2) = \frac{1}{ {\pi \sigma }} \left[1+\frac{(x - \mu )^2}{\sigma ^2} \right]^{-1}

柯西分布的展示效果及代码实现如下。

  1. import numpy as np
  2. import math
  3. from matplotlib import pyplot as plt
  4.  
  5. def cauchy(x, mu, sigma):
  6.     c = 1/(math.pi * sigma)
  7.     y = (1 + ((x -mu)/sigma)**2)**(-1) * c
  8.     return x, y, np.mean(y), np.std(y)
  9.  
  10. for ls in [(0, 1), (1, 2), (2, 0.5), (2, 5),(3, 2) ]:
  11.     a, b = ls[0], ls[1]
  12.  
  13.     x = np.arange(-4, 20, 0.01, dtype=np.float)
  14.     x, y, u, s = cauchy(x, mu = a, sigma =b)
  15.     plt.plot(x, y, label=r'$\mu=%.2f,\ \sigma=%.2f$'% (a, b))
  16. plt.legend()
  17. #plt.savefig('./graph/gamma.png')
  18. plt.show()

2.2.8 拉普拉斯分布(Laplace Distribution)

拉普拉斯分布又称为双侧指数分布，它的概率密度函数pdf表示如下：

                                                Lap(x|\mu,b)=\frac{1}{2b}exp(-\frac{|x-\mu|}{b})

其中参数\mu,b分别表示位置参数和比例参数（或缩放参数）,且b>0。相应地主要参数有：

                                                mean = \mu,mode=\mu,var=2b^2

拉普拉斯分布在0位置相比高斯分布有更多的概率密度设置，这个特性使得它增强稀疏的作用。对拉普拉斯取负对数可以得到L_1正则化的形式，而不同于高斯分布得到的L_2正则化形式，这有利于得到稀疏解。laplace分布的展示效果及代码实现如下。

  1. import numpy as np
  2. from matplotlib import pyplot as plt
  3.  
  4. def laplace(x, mu, b):
  5.     
  6.     y =  np.exp(-abs(x-mu)/b)/(2*b)
  7.     return x, y, np.mean(y), np.std(y)
  8.  
  9. for ls in [(0, 1), (1, 5), (3, 1), (2, 2)]:
  10.     a1, a2 = ls[0], ls[1]
  11.  
  12.     x = np.arange(-10, 10, 0.01, dtype=np.float)
  13.     x, y, u, s = laplace(x, mu=a1, b=a2)
  14.     plt.plot(x, y, label=r'$\mu=%d,\ b=%d$' % (a1, a2))
  15. plt.legend()
  16. #plt.savefig('./graph/gamma.png')
  17. plt.show()

2.2.9 伽马分布(Gamma Distribution)

gamma 分布是一类非常重要的基础分布，因为很多常用的分布函数可视作伽马分布的一个特例，它的分布函数pdf主要有两个参数形状参数\alpha 和尺度参数 \beta 来定义，具体如下

                                                Ga (x|\alpha,\beta) =\frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{\alpha-1}e^{-\beta x}, x> 0

                                                   \Gamma (z)=\int_{0 }^{\infty}u^{z-1}e^{-u}du

上式中，\Gamma(z)称之为gamma函数，根据定义可知伽马函数具有如下性质：

                Γ(z+1)=∫0∞uze−udu=[−uze−u]0∞+∫0∞zuz−1e−udu=limu→∞(−uze−u)−(−0ze−0)+z∫0∞uz−1e−udu=z∫0∞uz−1e−udu=zΓ(z)

gamma分布具有如下性质

                                                mean =\frac{\alpha}{\beta},mode=\frac{\alpha-1}{\beta},var =\frac{\alpha}{\beta^2}

当\alpha = n时，伽马分布可视为n个指数分布的和。

伽马分布实际上将指数分布与泊松分布联系在一起。

1.指数分布解决的问题是，要发生一个随机事件需要等待多长的时间

2.伽马分布解决的问题是，要发生n个随机事件需要等待多长时间（每个事件发生概率符合指数分布）

3.泊松分布解决的问题是在特定的时间内，发生n次随机事件的概率（小概率事件）

伽马分布的展示效果及代码如下。

  1. import numpy as np
  2. from matplotlib import pyplot as plt
  3.  
  4. def gamma_function(n):
  5.     cal = 1
  6.     for i in range(2, n):
  7.         cal *= i
  8.     return cal
  9.  
  10. def gamma(x, a, b):
  11.     c = (b ** a) / gamma_function(a)
  12.     y = c * (x ** (a - 1)) * np.exp(-b * x)
  13.     return x, y, np.mean(y), np.std(y)
  14.  
  15. for ls in [(1, 1), (2, 1), (3, 1), (2, 2)]:
  16.     a, b = ls[0], ls[1]
  17.  
  18.     x = np.arange(0, 20, 0.01, dtype=np.float)
  19.     x, y, u, s = gamma(x, a=a, b=b)
  20.     plt.plot(x, y, label=r'$\mu=%.2f,\ \sigma=%.2f,'
  21.                          r'\ \alpha=%d,\ \beta=%d$' % (u, s, a, b))
  22. plt.legend()
  23. #plt.savefig('./graph/gamma.png')
  24. plt.show()

2.2.10 指数分布(Exponential Distribution)

指数分布是gamma分布的一个特例，即它的pdf定义如下

                               Expon(x|\lambda)= \begin{cases} Ga(1,\lambda)= \frac{\lambda e^{-\lambda x}}{\Gamma(1)}=\lambda e^{-\lambda x}& \text{ if } x\geq 0 \\ 0 & \text{ if } x< 0 \end{cases}

 指数分布描述了泊松过程中事件之间的距离或时间的分布。这里，\lambda称之为rate参数，代表单位时间内的事件数量（比如每2h接到3个投诉电话或每10h设备发生一次故障），x 代表时间。指数分布与泊松分布存在一定的关联性，因为指数分布也是泊松过程的体现。泊松分布处理固定时间段内事件发生的次数(对次数建模)，属于离散性概率分布，而指数分布通常涉及到某个特定事件发生之前的时间量(对时间建模), 属于连续概率分布。

假设某个人通常每2小时接到 1 个电话。现在计算一个电话在下一小时内打来的概率？

解决方案： 假设每2小时 1 个电话。因此，预计每小时接到半个电话。所以可以取  λ = 0.5

因此，计算过程如下：

                        p(0\leq X\leq 1)=0.5exp(-0.5)=0.303265,

因此，下一小时内电话到达的概率为 0.303265.

指数分布的展示效果及代码如下。

  1.  
  2. import numpy as np
  3. from matplotlib import pyplot as plt
  4.  
  5. def exponential(x, lamb):
  6.     y = lamb * np.exp(-lamb * x)
  7.     return x, y, np.mean(y), np.std(y)
  8.  
  9. for lamb in [1, 2, 3]:
  10.  
  11.     x = np.arange(0, 6, 0.01, dtype=np.float)
  12.     x, y, u, s = exponential(x, lamb=lamb)
  13.     plt.plot(x, y, label=r'$\mu=%.2f,\ \sigma=%.2f,'
  14.                          r'\ \lambda=%d$' % (u, s, lamb))
  15. plt.legend()
  16. #plt.savefig('./graph/exponential.png')
  17. plt.show()

2.2.11 卡方分布(Chi-Squared Distribution)

k个独立的标准正态分布变量的平方和服从自由度为k的卡方分布。卡方分布是伽马分布的一个特例，相应的pdf定义如下：

                                      \chi(x|k)\triangleq Ga(x|\frac{k}{2},\frac{1}{2}) =\frac{\frac{1}{2}^{k/2}}{\Gamma(k/2)}x^{\frac{k}{2}-1}e^{-\frac{1}{2} x}

比如定义k个随机变量{Z_1,...,Z_k}是相互独立，符合标准正态分布的随机变量（数学期望为0、方差为1） Z_i \sim N(0,1),那么 S\sim \chi_k ^{2}，其中S的定义如下

                                                 S =\sum_{i=1}^{k}Z_i^{2}

卡方分布常用于对样本估计的假设检验与置信区间的计算。卡方分布的展示效果及代码如下。

  1. import numpy as np
  2. from matplotlib import pyplot as plt
  3.  
  4. def gamma_function(n):
  5.     cal = 1
  6.     for i in range(2, n):
  7.         cal *= i
  8.     return cal
  9.  
  10. def chi_squared(x, k):
  11.     
  12.     c = 1 / (2 ** (k/2) * gamma_function(k//2))
  13.     y = c * (x ** (k/2 - 1)) * np.exp(-x /2)
  14.  
  15.     return x, y, np.mean(y), np.std(y)
  16.  
  17. for k in [2, 3, 4, 6]:
  18.     x = np.arange(0, 10, 0.01, dtype=np.float)
  19.     x, y, _, _ = chi_squared(x, k)
  20.     plt.plot(x, y, label=r'$k=%d$' % (k))
  21.  
  22. plt.legend()
  23. #plt.savefig('./graph/chi-squared.png')
  24. plt.show()

2.2.12 贝塔分布(Beta Distribution)

贝塔分布是关于概率的概率分布，由于贝塔分布是对概率建模，所以它的定义域 x 的取值在0-1之间，相应的pdf对应如下:

                                        Beta(x|\alpha,\beta) =\frac{1}{B(\alpha,\beta)}x^{\alpha-1}(1-x)^{\beta-1}

其中系数项分母定义如下：

                ​​​​​​​        ​​​​​​​                                B(\alpha,\beta) \triangleq \frac{\Gamma(\alpha) \Gamma(\beta) }{\Gamma(\alpha+\beta)}

        ​​​​​​​        ​​​​​​​        mean =\frac{a}{\alpha+\beta} ,mode=\frac{\alpha-1}{\alpha+\beta-2},var =\frac{\alpha \beta}{(\alpha+\beta)^2(\alpha+\beta+1)}

观察贝塔分布的前置项，可以视为满足概率分布的归一化参数，直观的对比贝塔与二项分布：

 二项分布和 beta 分布的区别在于前者对次数建模，模拟成功的次数 (x)；而后者对概率建模，模拟成功的概率 (p)。 换句话说，概率是二项分布的参数，在 Beta分布中，概率是一个随机变量。

关于参数\alpha,\beta的解释：

      可以将α-1 视为成功的次数，将β -1 视为失败的次数，就像二项式中的x和n-x项一样。 可以选择 α 和 β 参数，如果认为成功的概率非常高，假设为 90%，则将 α 设置为 90，将β 设置为 10。反之，如果不这么认为，β 为 90，α 为 10。 随着α变大（更成功的事件），大部分概率分布将向右移动，而β的增加使分布向左移动（更多失败）。 此外，如果α和β都增加，分布将变窄，因为更加确定。

贝塔分布与伽马分布一样，也是一种非常常用的基础分布，对于贝塔分布的一个特例，即\alpha = \beta =1,我们就能得到均匀分布。贝塔分布的展示效果及代码实现如下。

 

  1. import numpy as np
  2. from matplotlib import pyplot as plt
  3.  
  4. def gamma_function(n):
  5.     cal = 1
  6.     for i in range(2, n):
  7.         cal *= i
  8.     return cal
  9.  
  10. def beta(x, a, b):
  11.  
  12.     gamma = gamma_function(a + b) / \
  13.             (gamma_function(a) * gamma_function(b))
  14.     y = gamma * (x ** (a - 1)) * ((1 - x) ** (b - 1))
  15.     return x, y, np.mean(y), np.std(y)
  16.  
  17. for ls in [(1, 3), (5, 1), (2, 2), (2, 5)]:
  18.     a, b = ls[0], ls[1]
  19.  
  20.     # x in [0, 1], trial is 1/0.001 = 1000
  21.     x = np.arange(0, 1, 0.001, dtype=np.float)
  22.     x, y, u, s = beta(x, a=a, b=b)
  23.     plt.plot(x, y, label=r'$\mu=%.2f,\ \sigma=%.2f,'
  24.                          r'\ \alpha=%d,\ \beta=%d$' % (u, s, a, b))
  25. plt.legend()
  26. #plt.savefig('./graph/beta.png')
  27. plt.show()

2.2.13 贝塔-二项分布(Beta-Binomial Distribution)

作为一种复合分布，beta-binomial是二项分布的共轭分布，可以认为二项分布中的参数p是从 beta 分布中随机抽取的,即

        ​​​​​​​        ​​​​​​​                        ​​​​​​​        f(x\mid p)=\binom{n}{k} p^x (1-p)^{n-x}

上式中p不是一个常量,而是作为变量且服从如下贝塔分布

        ​​​​​​​        ​​​​​​​                         f(p)=\frac{\Gamma (\alpha+\beta)}{\Gamma (\alpha)\Gamma(\beta)}p^{\alpha-1}(1-p)^{\beta-1},0< x<1,\alpha> 0,\beta> 0,

由此得到贝塔-二项分布的pdf如下

     ​​​​​​​                        f(x∣n,α,β)=∫01Bin(x|n,p)Beta(p∣α,β)dp=(nx)1B(α,β)∫01px+α−1(1−p)n−x+β−1dp=(nx)B(x+α,n−x+β)B(α,β)

使用beta 函数的性质，上式可以写成如下

        ​​​​​​​        f(x\mid n,\alpha,\beta) = \frac{\Gamma(n+1)}{\Gamma(x+1)\Gamma(n-x+1)} \frac{\Gamma(x+\alpha)\Gamma(n-x+\beta)}{\Gamma(n+\alpha+\beta)} \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}

贝塔-二项分布的展示效果及代码实现如下。

  1. import numpy as np
  2. from matplotlib import pyplot as plt
  3.  
  4. def gamma_function(n):
  5.     cal = 1
  6.     for i in range(2, n):
  7.         cal *= i
  8.     return cal
  9.  
  10. def beta_binom(n, a, b):
  11.     
  12.     multipliers3 = gamma_function(a + b) / \
  13.             (gamma_function(a) * gamma_function(b))
  14.             
  15.     xlist = list(range(1,n))
  16.     ylist = []
  17.     
  18.     for x in xlist:
  19.     
  20.         multipliers1 = gamma_function(n + 1) / \
  21.                 (gamma_function(x+1) * gamma_function(n-x+1))
  22.                 
  23.         multipliers2 = gamma_function(x + a) * gamma_function(n-x+b) / \
  24.                 (gamma_function(n+a+b))
  25.             
  26.         y = multipliers1 *multipliers2 * multipliers3
  27.         ylist.append(y)
  28.         
  29.     return xlist, ylist, np.mean(ylist), np.std(ylist)
  30.  
  31. for ls in [(1, 3, 15), (5, 1, 15), (2, 2, 15), (2, 5, 15)]:
  32.     a, b, n = ls[0], ls[1], ls[2]
  33.  
  34.     # x in [0, 1], trial is 1/0.001 = 1000
  35.     xn = np.arange(0, 100, 1, dtype=np.int32)
  36.     x, y, u, s = beta_binom(n=n, a=a, b=b)
  37.     plt.plot(x, y, label=r'$\mu=%.2f,\ \sigma=%.2f,'
  38.                          r'\ \alpha=%d,\ \beta=%d$' % (u, s, a, b))
  39. plt.legend()
  40. #plt.savefig('./graph/beta.png')
  41. plt.show()

2.2.14 帕累托分布(Pareto Distribution)

帕累托分布常用于模拟具有明显长尾特性的分布，对应的pmf如下。

                                Pareto(x|k,m) = km^kx^{-(k+1)}\mathbb{I}(x\geqslant m)

        ​​​​​​​        ​​​​​​​     mean =\frac{km}{k-1} \ if \ {k>1},mode=m,var =\frac{m^2k}{(k-1)^2(k-2)} \ if \ {k>2}

帕累托分布的展示效果及代码如下。

  1. import numpy as np
  2. import matplotlib.pyplot as plt
  3. k, m = 6., 2.  # distribution width、specify mode   
  4. s = (np.random.pareto(k, 1000) + 1) * m
  5. # draw histgram
  6. count, bins, _ = plt.hist(s, 100, density=True)
  7.  
  8. # draw line graph  
  9. fit = k * m ** k / (bins ** (k + 1))
  10. plt.plot(bins, max(count) * fit / max(fit), linewidth=2, color='r')
  11. plt.show()

2.2.15 狄里克雷分布(Dirichlet Distribution)

Dirichlet分布可以看做是分布之上的分布，是（多项）分布的概率的概率，是一类在实数域以正单纯形（standard simplex）为支撑集（support）的高维连续概率分布，是Beta分布在高维情形的推广，在LDA文本数据挖掘中得到了广泛应用。Dirichlet分布的pdf如下

                                        S_K =\{\mathbf x: 0\leq x \leq 1, \sum_{k=1}^{K}x_k=1\}

                                        \(Dir(\mathbf x| \boldsymbol \alpha) =\frac{1}{B(\boldsymbol \alpha)}\coprod_{k=1}^{K}x_k^{\alpha_k-1}\mathbb{I}(\mathbf x \in S_K)\\)​​​​​​​

其中，

                        {B(\boldsymbol \alpha)}=\frac{\coprod_{k=1}^{K}\Gamma (\alpha_k)}{\Gamma (\alpha_0)}=\frac{\Gamma (\alpha_1)\Gamma (\alpha_2)...\Gamma (\alpha_K)}{\Gamma (\alpha_1+\alpha_2+...+\alpha_K)}

当K=2时，上式退化为贝塔分布的归一化因子，其中  \alpha_0 = \alpha_1+\alpha_2+...+\alpha_K 。

Dirichlet分布的展示效果及代码如下。

 

  1. from random import randint
  2. import numpy as np
  3. from matplotlib import pyplot as plt
  4.  
  5. def normalization(x, s):
  6.     """
  7.     :return: normalizated list, where sum(x) == s
  8.     """
  9.     return [(i * s) / sum(x) for i in x]
  10.  
  11. def sampling():
  12.     return normalization([randint(1, 100),
  13.             randint(1, 100), randint(1, 100)], s=1)
  14.  
  15. def gamma_function(n):
  16.     cal = 1
  17.     for i in range(2, n):
  18.         cal *= i
  19.     return cal
  20.  
  21. def beta_function(alpha):
  22.     """
  23.     :param alpha: list, len(alpha) is k
  24.     :return:
  25.     """
  26.     numerator = 1
  27.     for a in alpha:
  28.         numerator *= gamma_function(a)
  29.     denominator = gamma_function(sum(alpha))
  30.     return numerator / denominator
  31.  
  32. def dirichlet(x, a, n):
  33.     """
  34.     :param x: list of [x[1,...,K], x[1,...,K], ...], shape is (n_trial, K)
  35.     :param a: list of coefficient, a_i > 0
  36.     :param n: number of trial
  37.     :return:
  38.     """
  39.     c = (1 / beta_function(a))
  40.     y = [c * (xn[0] ** (a[0] - 1)) * (xn[1] ** (a[1] - 1))
  41.          * (xn[2] ** (a[2] - 1)) for xn in x]
  42.     x = np.arange(n)
  43.     return x, y, np.mean(y), np.std(y)
  44.  
  45. n_experiment = 1200
  46. for ls in [(6, 2, 2), (3, 7, 5), (6, 2, 6), (2, 3, 4)]:
  47.     alpha = list(ls)
  48.  
  49.     # random samping [x[1,...,K], x[1,...,K], ...], shape is (n_trial, K)
  50.     # each sum of row should be one.
  51.     x = [sampling() for _ in range(1, n_experiment + 1)]
  52.  
  53.     x, y, u, s = dirichlet(x, alpha, n=n_experiment)
  54.     plt.plot(x, y, label=r'$\alpha=(%d,%d,%d)$' % (ls[0], ls[1], ls[2]))
  55.  
  56. plt.legend()
  57. #plt.savefig('./graph/dirichlet.png')
  58. plt.show()

2.2.16 多变量分布(Multivariate Distribution)

Multivariate Gaussian

        ​​​​​​​        ​​​​​​​        ​​​​​​​        \large \mathcal N(\mathbf x\mid \boldsymbol \mu,\boldsymbol \Sigma) = \frac{1}{(2\pi)^{D/2} \left |\boldsymbol \Sigma \right |^{1/2}} exp\left [ -\frac{1}{2}({\mathbf x}-{\boldsymbol\mu})^\mathrm{T}{\boldsymbol\Sigma}^{-1}({\mathbf x}-{\boldsymbol\mu}) \right ]​​​​​​​

Multivariate Student      

                    \large \mathcal \tau (\mathbf x|\boldsymbol \mu,\boldsymbol \Sigma,\nu) = \frac{\Gamma\left[(\nu+D)/2\right]}{\Gamma(\nu/2)\nu^{D/2}\pi^{D/2}\left|{\boldsymbol\Sigma}\right|^{1/2}}\left[1+\frac{1}{\nu}({\mathbf x}-{\boldsymbol\mu})^T{\boldsymbol\Sigma}^{-1}({\mathbf x}-{\boldsymbol\mu})\right]^{-(\nu+D)/2}

2.2.17 其它分布

其它分布还有Erlang分布、逆伽马（Inverse Gamma）分布、逆高斯（Inverse Gaussian）分布，对应分布函数可以视为上述分布的一种特殊形式或变体，具体的含义不再详述。

        ​​​​​​​        ​​​​​​​        ​​​​​​​        ​​​​​​​        Erlang(x \mid \lambda)=Ga(x \mid 2, \lambda)

        ​​​​​​​        ​​​​​​​        ​​​​​​​        ​​​​​​​        IGa(x \mid a,b) \triangleq \frac{b^{a}}{\Gamma(a)}x^{-(a+1)}e^{-b/x}

        ​​​​​​​        ​​​​​​​                        IGau(x \mid \mu,\mu^2) \triangleq \frac{\mu}{\sqrt{2 \pi x^3}} \exp\biggl(-\frac{(x-\mu)^2}{2x}\biggr)

3.概率分布关系及总结

4.参考文献

1. MLAPP

2.PRML

3.常见12种概率分布

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/知新_RL/article/detail/438873
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号