1. Gamma 分布
参考https://zhuanlan.zhihu.com/p/37976562
模型假设事件单位时间内发生 次,则发生 次所经过的时间
1.1 Gamma 函数
Gamma函数定义为
形象理解为用一个伽马刀,对
动了一刀,于是指数为
,动完刀需要扶着梯子
才能走下来。
Gamma 函数是阶乘的推广具有性质
1.2 Gamma 函数的可视化
import numpy as np
from scipy.special import gamma
import matplotlib.pyplot as plt
import pylab
fig = plt.figure(figsize=(12,8))
# The Gamma function
x = np.linspace(-5, 5, 1000)
plt.plot(x, gamma(x), ls='-', c='k', label='$\Gamma(x)$')
# (x-1)! for x = 1, 2, ..., 6
x2 = np.linspace(1,6,6)
y = np.array([1, 1, 2, 6, 24, 120])
pylab.plot(x2, y, marker='*', markersize=12, markeredgecolor='r',
markerfacecolor='r', ls='',c='r', label='$(x-1)!$')
plt.title('Gamma Function')
plt.ylim(-50,50)
plt.xlim(-5, 5)
plt.xlabel('$x$')
plt.legend()
plt.show()
从图可以看出Gamma函数在 是凸函数,不仅如此, 也是凸函数。
fig = plt.figure(figsize=(12,8))
# The Gamma function
x = np.linspace(0, 15, 1000)
plt.plot(x, np.log(gamma(x)), ls='-', c='k', label='$log\Gamma(x)$')
plt.title('Log$\Gamma(x)$ Function')
plt.ylim(-1,50)
plt.xlim(-1, 15)
plt.xlabel('$x$')
plt.legend()
plt.show()
1.3. Gamma 分布可视化
import numpy as np
from scipy.stats import gamma
from matplotlib import pyplot as plt
alpha_values = [1, 2, 3, 3, 3]
beta_values = [0.5, 0.5, 0.5, 1, 2]
color = ['b','r','g','y','m']
x = np.linspace(1E-6, 10, 1000)
fig, ax = plt.subplots(figsize=(12, 8))
for k, t, c in zip(alpha_values, beta_values, color):
dist = gamma(k, 0, t)
plt.plot(x, dist.pdf(x), c=c, label=r'$\alpha=%.1f,\ \beta=%.1f$' % (k, t))
plt.xlim(0, 10)
plt.ylim(0, 2)
plt.xlabel('$x$')
plt.ylabel(r'$p(x|\alpha,\beta)$')
plt.title('Gamma Distribution')
plt.legend(loc=0)
plt.show()
2. 分布
模型:
-
假设 服从 。
-
把 个随机变量排序后得到顺序统计量 。
-
第 小的 所服从的分布
2.1 模型的推导
落入两个以上的数
2.2. 分布的可视化
import numpy as np
from scipy.stats import beta
from matplotlib import pyplot as plt
alpha_values = [1/3,2/3,1,1,2,2,4,10,20]
beta_values = [1,2/3,3,1,1,6,4,30,20]
colors = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple',
'tab:brown', 'tab:pink', 'tab:gray', 'tab:olive']
x = np.linspace(0, 1, 1002)[1:-1]
fig, ax = plt.subplots(figsize=(14,9))
for a, b, c in zip(alpha_values, beta_values, colors):
dist = beta(a, b)
plt.plot(x, dist.pdf(x), c=c,label=r'$\alpha=%.1f,\ \beta=%.1f$' % (a, b))
plt.xlim(0, 1)
plt.ylim(0, 6)
plt.xlabel('$x$')
plt.ylabel(r'$p(x|\alpha,\beta)$')
plt.title('Beta Distribution')
ax.annotate('Beta(1/3,1)', xy=(0.014, 5), xytext=(0.04, 5.2),
arrowprops=dict(facecolor='black', arrowstyle='-'))
ax.annotate('Beta(10,30)', xy=(0.276, 5), xytext=(0.3, 5.4),
arrowprops=dict(facecolor='black', arrowstyle='-'))
ax.annotate('Beta(20,20)', xy=(0.5, 5), xytext=(0.52, 5.4),
arrowprops=dict(facecolor='black', arrowstyle='-'))
ax.annotate('Beta(1,3)', xy=(0.06, 2.6), xytext=(0.07, 3.1),
arrowprops=dict(facecolor='black', arrowstyle='-'))
ax.annotate('Beta(2,6)', xy=(0.256, 2.41), xytext=(0.2, 3.1),
arrowprops=dict(facecolor='black', arrowstyle='-'))
ax.annotate('Beta(4,4)', xy=(0.53, 2.15), xytext=(0.45, 2.6),
arrowprops=dict(facecolor='black', arrowstyle='-'))
ax.annotate('Beta(1,1)', xy=(0.8, 1), xytext=(0.7, 2),
arrowprops=dict(facecolor='black', arrowstyle='-'))
ax.annotate('Beta(2,1)', xy=(0.9, 1.8), xytext=(0.75, 2.6),
arrowprops=dict(facecolor='black', arrowstyle='-'))
ax.annotate('Beta(2/3,2/3)', xy=(0.99, 2.4), xytext=(0.86, 2.8),
arrowprops=dict(facecolor='black', arrowstyle='-'))
#plt.legend(loc=0)
plt.show()
2.3. 分布的应用
棒球击球率
那么我们简单说个Beta-Binomial共轭的应用。用一句话来说,beta分布可以看作一个概率的概率分布,当你不知道一个东西的具体概率是多少时,它可以给出了所有概率出现的可能性大小。
举一个简单的例子,熟悉棒球运动的都知道有一个指标就是棒球击球率(batting average),就是用一个运动员击中的球数除以击球的总数,我们一般认为0.266是正常水平的击球率,而如果击球率高达0.3就被认为是非常优秀的。现在有一个棒球运动员,我们希望能够预测他在这一赛季中的棒球击球率是多少。传统的频率学派会直接计算棒球击球率,用击中的数除以击球数,但是如果这个棒球运动员只打了一次,而且还命中了,那么他就击球率就是100%了,这显然是不合理的,因为根据棒球的历史信息,我们知道这个击球率应该是0.215到0.36之间才对。对于这个问题,我们可以用一个二项分布表示(一系列成功或失败),一个最好的方法来表示这些经验(在统计中称为先验信息)就是用beta分布,这表示在我们没有看到这个运动员打球之前,我们就有了一个大概的范围。beta分布取
import numpy as np
from scipy.stats import beta
from matplotlib import pyplot as plt
x = np.linspace(0, 1, 1002)[1:-1]
fig, ax = plt.subplots(figsize=(10,6))
dist = beta(81, 219)
plt.plot(x, dist.pdf(x), c='b',label=r'$\alpha=%.1f,\ \beta=%.1f$' % (81, 219))
plt.xlim(0, .6)
plt.ylim(0, 16)
plt.xlabel('$x$')
plt.ylabel(r'$p(x|\alpha,\beta)$')
plt.title('Beta Distribution')
plt.legend(loc=0)
plt.show()
可以看到这个分布其实没多大变化,这是因为只打了1次球并不能说明什么问题。但是如果我们得到了更多的数据,假设一共打了300次,其中击中了100次,200次没击中,那么这一新分布就是:
import numpy as np
from scipy.stats import beta
from matplotlib import pyplot as plt
x = np.linspace(0, 1, 1002)[1:-1]
fig, ax = plt.subplots(figsize=(10,6))
dist = beta(81, 219)
plt.plot(x, dist.pdf(x), c='b',label=r'$\alpha=%.1f,\ \beta=%.1f$' % (81, 219))
dist = beta(181, 419)
plt.plot(x, dist.pdf(x), c='r',label=r'$\alpha=%.1f,\ \beta=%.1f$' % (181, 419))
plt.xlim(0, .6)
plt.ylim(0, 22)
plt.xlabel('$x$')
plt.ylabel(r'$p(x|\alpha,\beta)$')
plt.title('Beta Distribution')
plt.legend(loc=0)
plt.show()
3. Dirichlet分布
3.1 Dirichlet分布的模型
模型:
- 假设 服从 。
- 把 个随机变量排序后得到顺序统计量 。
- 所服从的分布就是Dirichlet分布
3.2 可视化Dirichlet分布
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as tri
#生成等边三角形
corners = np.array([[0, 0], [1, 0], [0.5, 0.75**0.5]])
triangle = tri.Triangulation(corners[:, 0], corners[:, 1])
#每条边中点位置
midpoints = [(corners[(i + 1) % 3] + corners[(i + 2) % 3]) / 2.0 for i in range(3)]
def xy2bc(xy, tol=1.e-3):
#将三角形顶点的笛卡尔坐标映射到重心坐标系
s = [(corners[i] - midpoints[i]).dot(xy - midpoints[i]) / 0.75 for i in range(3)]
return np.clip(s, tol, 1.0 - tol)
#有了重心坐标,可以计算Dirichlet概率密度函数的值
class Dirichlet(object):
def __init__(self, alpha):
from math import gamma
from operator import mul
from functools import reduce
self._alpha = np.array(alpha)
self._coef = gamma(np.sum(self._alpha)) / reduce(mul, [gamma(a) for a in self._alpha]) #reduce:sequence连续使用function
def pdf(self, x):
#返回概率密度函数值
from operator import mul
from functools import reduce
return self._coef * reduce(mul, [xx ** (aa - 1) for (xx, aa)in zip(x, self._alpha)])
def draw_pdf_contours(dist, nlevels=200, subdiv=8, **kwargs):
import math
#细分等边三角形网格
refiner = tri.UniformTriRefiner(triangle)
trimesh = refiner.refine_triangulation(subdiv=subdiv)
pvals = [dist.pdf(xy2bc(xy)) for xy in zip(trimesh.x, trimesh.y)]
plt.tricontourf(trimesh, pvals, nlevels, **kwargs)
plt.axis('equal')
plt.xlim(0, 1)
plt.ylim(0, 0.75**0.5)
plt.axis('off')
让我们看看几个对称的Dirichlet分布,第一个 产生一个均匀分布,所有点在simplex上概率相同。
draw_pdf_contours(Dirichlet([1, 1, 1]))
对于 ,分布集中在角落和和simplex边界。(黄色概率最大,蓝色概率最小)
draw_pdf_contours(Dirichlet([0.999, 0.999, 0.999]))
对于 ,分布集中在simplex中心。
draw_pdf_contours(Dirichlet([5, 5, 5]))