Beta分布从入门到精通

最近一直有点小忙，但是不知道在瞎忙什么，终于有时间把Beta分布的整理弄完。

下面的内容，夹杂着英文和中文，呵呵～

Beta Distribution

Beta Distribution Definition:

The Beta distribution is a special case of the Dirichlet distribution, and is related to the Gamma distribution. It has the probability distribution function:

这里，因为Beta分数是二项分布的参数p的概率分布, 所以x（即p）的取值范围为0 <= x <= 1

where the normalisation, B, is thebeta function, Beta function could also be expressed by Gamma function:

Gamma函数在实数域可以表示为：

Gamma函数在整数域可以表示为：

Γ(n)=(n?1)!

Gamma函数有以下性质：

因为Beta函数可以表示为Gamma函数，所以Beta分布还可以表示为：

0 <= x <= 1

Beta分布可以理解为二项分布的参数p的分布，所以，这里重新定义Beta分布：

Beta分布的期望：

Beta分布的方差：

Beta分布的众数 mode:

Beta分布的偏度 Skewness：

Beta分布的峰度 Kurtosis:

Beta Distribution Examples

Beta分布可以说是一个百变星君，根据参数a,b的不同，可以呈现出多种完全不同的概率分布图.

生成Beta分布的代码：

from scipy.stats import beta
import matplotlib.pyplot as plt
import numpy as np

a, b = 2, 1

mean, var, skew, kurt = beta.stats(a, b, moments='mvsk')

x = np.linspace(0, 1, 100)
plt.plot(x, beta.pdf(x, a, b), 'r-', lw=5, alpha=0.6, label='beta pdf')
plt.show()

然后，根据调整代码中的a,b的取值，可以得到不同的Beta分布：

a, b = 2, 1:

a, b = 2, 2

a, b = 8, 2

a, b = 0.01, 20

a, b = 1, 1

这样一个一个的绘制，是不是太逊了，画在一起：

代码：

from scipy.stats import beta
import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(0, 1, 100)

a_array = [1, 2, 4, 8]
b_array = [1, 2, 4, 8]

fig, axarr = plt.subplots(len(a_array), len(b_array))

for i, a in enumerate(a_array):
  for j, b in enumerate(b_array):
    axarr[i, j].plot(x, beta.pdf(x, a, b), 'r', lw=1, alpha=0.6, label='a='+str(a)+',b='+str(b))
    axarr[i, j].legend(frameon=False)

plt.show()

将所有的Beta分布绘制在一个图上：

代码：

from scipy.stats import beta
import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(0, 1, 100)

a_array = [1, 2, 4, 8]
b_array = [1, 2, 4, 8]

for i, a in enumerate(a_array):
  for j, b in enumerate(b_array):
    plt.plot(x, beta.pdf(x, a, b), lw=1, alpha=0.6, label='a='+str(a)+',b='+str(b))

plt.legend(frameon=False)
plt.show()

Beta Mean

由公式可以得到，Beta分布的均值，也可以通过采样的方法，在一个Beta分布中，采样，计算均值。

代码：

import numpy as np
import numpy.random as nprnd
import scipy.stats as spstat
import scipy.special as ssp
import itertools as itt

import matplotlib.pyplot as plt
import pylab as pl

N = (np.arange(200) + 3) ** 2 * 20

betamean = np.zeros_like(N, dtype=np.float64)

for idx, i in enumerate(N):
    betamean[idx] = np.mean(nprnd.beta(2, 1, i))

plt.plot(N, betamean, color='steelblue', lw=2)
plt.xscale('log')
plt.show()

print spstat.beta(2, 1).mean()

print spstat.beta(2, 1).mean(), 2.0 / (2 + 1)
print spstat.beta(2, 1).var(), 2 * 1.0 / (2 + 1 + 1) / (2 + 1) ** 2

运行结果：

这里可以看到，随着采样点的增加，样本点的均值也就更加的收敛，更加的接近?, ? 是一个通过公式计算得到的。这样，这个图片的结果也符合大数定理，随着采样点的增加，只要样本点无限大，那么最终的均值就会无限的接近?.

Conjugate Prior

A conjugate prior,p(p), of a likelihood, p(x|p), is a distribution that results in a posterior distribution, p(p|x)with the same functional form as the prior
and a parameterisation that incorporates the observationx.

这句话，猛的一读，晕头转向，但是，仔细读上三五遍，基本上就理解了什么叫“共轭先验”。

基本上说，一个参数的共轭先验p(p)是这样的一个分布：在这个分布的基础上加上观测样本能够得到一个与先验分布具有相同的函数形式的后验概率分布p(p|x)，并且这个后验概率分布p(p|x)融合了观测样本x。也就是说共轭先验p(p)和后验概率分布p(p|x)具有相当的函数形式。

说点人话吧。。。

Beta分布是二项分布的参数p的共轭先验，也就是说，二项分布的参数p的共轭先验是一个Beta分布，其中，Beta分布中的两个参数a,b可以看作两个二项分布的参数p的先验知识，可以称为伪计数，例如
a, b = 2, 1, 这就意味着，二项分布的参数p的先验知识为：在三次实验中，a出现两次，b出现1次，也可以理解为发生了2次，没有发生的有1次。

后验概率也符合Beta分布：

Beta(p|a, b) + count(m1, m2) = Beta(p| a+m1, b+m2)

在二项分布的参数的先验分布的基础上，加上观测数据，就可以得到二项分布的参数p的后验概率分布也符合Beta分布。这里，
m1, m2 分别表示对应于 x=1 和 x=0在观测数据中出现的次数。

话说，共轭先验中的参数即Beta分布中的两个参数a,b 是非常有意义的hyperparameter的解释，前面已经提到了，a,b
可以理解为在观测样本 (m1, m2)的基础上的先验知识，或者可以理解为伪计数，即在我们的先验知识中， x=1和x=0分别应该出现多少次，并且，这个先验知识的取值，对于后验概率的计算有比较大的影响。

二项分布的参数p的后验概率分布仍然符合Beta分布可以通过下面的公式推到进行证明：