从数据集中进行分布参数估计: 以伯努利分布为例

大数据机器学习课程第一次作业

第一题

题目如下

说明伯努利模型的极大似然估计以及贝叶斯估计中的统计学习方法三要素。伯努利模型是定义在取值为0与1的随机变量上的概率分布。假设观测到伯努利模型n次独立的数据生成结果，其中k次的结果为1，这时可以用极大似然估计或贝叶斯估计来估计结果为1的概率。

审题

• 题目出处？（题目中的概念是哪一本书的定义？）
  • 李航《统计学习方法》
• 李航对于“模型”的定义是什么？
  • 见下一个疑问
• 李航的书里是怎么定义“统计学习方法三要素”的？
  • 方法=模型+策略+算法
  • 模型 model
    • 书上定义：所要学习的条件概率分布或者决策函数；需要给出一个假设空间来包括所有可能的这个分布或者函数。
    • 相当于定义了输入输出，定义搜索的空间
    • 这一道题中“伯努利模型”意思是这个模型要学习的就是一个不知道参数的伯努利分布
  • 策略 strategy
    • 书上定义：按照什么样的准则从假设空间中学习或者选择最优的模型。损失函数度量一次预测的好坏，风险函数度量平均意义下的好坏。
    • 相当于定义了“How to measure?”
    • 风险函数是联合分布下期望的损失函数值
  • 算法 algorithm
    • 书上定义：根据策略，通过计算方法求解最优模型，通过不同的最优化算法，比如解析解或者数值计算。
    • 相当于定义了How to search/Optimize“
• 伯努利分布 Bernoulli distribution 复习
  • =两点分布=0-1分布
  • 伯努利试验：单次实验的结果或1或0；与泊松实验不同

解题

1. 伯努利模型的极大似然估计

1.1 三要素的识别

• 模型

  • 所要学习的分布：\(X\sim Bern(p)\)
    • P(X=1) = p, P(X=0) = 1-p
    • 统一一下形式: \(P(x) = p^x(1-p)^{(1-x)}\)， \(x\in\set{0, 1}\)
      • 注：如果p=1或者0，这个式子定义不良，我们需要补充规定\(0^0=1\)。
  • 参数空间: \(p\in[0, 1]\)
  • 假设空间：\(\mathcal{F}=\{P_p|P_p(x)=p^x(1-p)^{(1-x)}, p\in [0, 1]\}\)
    • 注：李航书16页中的集合表示是错误的
      • 如果按照李航书中的表示，这里我们会写成\(\mathcal{F}=\{P|P_p(x)=p^x(1-p)^{(1-x)}, p\in [0, 1]\}\)
      • 然而，这个集合只有一个元素“P”，P是固定的单一函数（是关于p和x的二元函数），形式已经完全确定了（就是上面的伯努利分布的形式）。
      • 只有一个元素自然不能表示“所有可能的条件概率分布或决策函数”，我们必须用\(P_p\)写在左边，右边指出了p具有取值范围，这才正确地用集合表示法获得了所有分布的一个集合。
  • 注：这里条件概率分布没有条件（或者说条件为空），因为我们是在做density estimation

• 策略

  • 策略是说，给定一个猜测的参数p，还有观测数据\(x_{1:n}\)，可以评价p猜测的是否合理。
  • 最大似然估计用似然概率的大小来评价
  • 似然概率是指，给定p之后，观测数据出现的概率
  • \(Like(p) = P(x_{1:n} | p)\)， 由于数据i.i.d.,
  • \(Like(p) =\prod_{i=1}^n P(x_i|p) = \prod_{i=1}^n p^{x_i}(1-p)^{(1-x_i)}\)
  • 题目中，没有告诉我们每个\(x_i\)是多少，只是告诉我们\(\sum_{i=1}^{n}x_i = k\)
  • 不过这个也足够我们化简上式为\(Like(p) = p^{k}(1-p)^{(n-k)}\)
  • 注：\(\binom{n}{k}p^k(1-p)^{(n-k)}\)是典型的错误答案，虽然不影响求解的结果
    • 。错误答案中的概率表示做n次伯努利实验得到k次正例的所有情况
    • 但是本题中我们只进行了一次“n次伯努利实验”，有具体的实验结果\(x_{1:n}\)， 不是进行了多次”n次伯努利实验”然后统计\(k_{1:m}\)，是其中具体的一次情况。
    • 那样的话我们的模型就变成了二项分布，而不是本题的伯努利分布模型。

• 算法

  • 极大似然估计的算法就是找到\(p\)让似然概率最大化。
  • \(\displaystyle \mathop{\arg\max}_{p} Like(p)= p^{k}(1-p)^{(n-k)}\)
  • 对于本模型可以获得解析解，具体1.2所示。

1.2 具体的求解

求函数\(Like(p)\)的最小值，首先显然\(Like(p)\)是连续可导的，所以最小值一定在critical point或者边缘0, 1上。

\(Like^\prime(p)= kp^{k-1}(1-p)^{(n-k)}-(n-k)p^k(1-p)^{(n-k-1)}\)

我们还可以学习一下python的sympy库辅助我们求导，防止人工求导抄错了

from sympy import *

p, n, k = symbols('p n k')
like = p**k * (1-p)**(n-k)
like_prime = diff(like, p)
like_prime

\(\displaystyle \frac{k p^{k} \left(1 - p\right)^{- k + n}}{p} + \frac{p^{k} \left(1 - p\right)^{- k + n} \left(k - n\right)}{1 - p}\)

令\(Like^{\prime}(p) = 0\) 。这些p称为critical point，这些点可能是极值点可能不是。

\(kp^{k-1}(1-p)^{(n-k)}=(n-k)p^k(1-p)^{(n-k-1)}\)

\(k(1-p) = (n-k)p\)

\(k-kp = (n-k)p\)

\(k = np\)

\(p=\frac{n}{k}\) 为唯一解

同样的，我们可以用Python做验证

eq = Eq(like_prime, 0)
solutions:list = solve(eq, p) 
assert len(solutions)!=0
po = solutions[0]
po

\(\displaystyle \frac{k}{n}\)

代入并且比较Like(0), Like(1)和Like(p)， 发现\(Like(p)\ge Like(1)\and Like(p)\ge Like(0)\)，所以\(Like(p)\)就是最小值（之一，如果n=k或者k=0）。

此时，\(Like(p) = \left(\frac{k}{n}\right)^{k} \left(\frac{- k + n}{n}\right)^{- k + n}\)

Python可以代入p为求解出来的p，得到上面的表达式。

like_po = like.subs(p, po)
like_po

\(\displaystyle \left(\frac{k}{n}\right)^{k} \left(- \frac{k}{n} + 1\right)^{- k + n}\)

2. 伯努利模型的贝叶斯估计

注：在PRML书中，贝叶斯估计也叫做Maximum a Posterior(MAP, 极大后验估计)

2.1 三要素的识别

• 模型
  • 与上文1.1一致，都是伯努利模型，不再赘述。
• 策略
  • 策略是说，给定一个猜测的参数p，还有观测数据\(x_{1:n}\)，可以评价p猜测的是否合理。
  • 最大似然估计用p的似然概率（实际上是观测数据的出现概率）来估计，而贝叶斯估计用p的后验概率来估计。
  • 后验概率是\(P(p|x_{1:n})\), 根据贝叶斯公式
  • \(P(p|x_{1:n}) =\frac{P(x_{1:n}|p)*P(p)}{P(x_{1:n})} = \frac{likelihood*prior}{evidence}\)
• 算法
  • 极大后验估计的算法就是找到\(p\)让后验概率最大化。
  • \(\displaystyle \mathop{\arg\max}_{p} P(p|x_{1:n})= \frac{P(x_{1:n}|p)*P(p)}{P(x_{1:n})}\)
  • 对于本模型，对于不同的先验概率P(p)会有不同的情况，有一些可以获得解析解，具体2.2所示。

2.2 具体的求解

由于题目中没有给出p的先验概率分布，

2.2.1 不妨认为参数p服从均匀分布，\(p\sim U[0, 1]\), 也就是说p的先验概率为均匀分布，那么$P(p) = =1,p$

• 注，这P是连续随机变量的概率密度函数

\(P(p|x_{1:n})\)的分母是一个正的相对于p的常数，对于求最大值而言可以去除。

那么此时原优化目标变为 \(\displaystyle \mathop{\arg\max}_{p} P(x_{1:n}|p)*P(p) = P(x_{1:n}|p)\)， 贝叶斯估计退化为极大似然估计，我们得到的结论和过程与上文1.2所述完全一致。

2.2.2 如果认为参数p服从\(\beta\)分布，我们能获得更加有意思的结果，这个在PRML书中有讲到。

\(\beta\)分布的密度函数是这样的\(P(p) = \beta(p;a,b) = \frac{p^{a-1}(1-p)^{b-1}}{C}\), a和b是这个分布的参数，C是使得概率密度函数积分为1的一个常数，在这里不重要。

\(\displaystyle \mathop{\arg\max}_{p} P(x_{1:n}|p)*P(p) = P(x_{1:n}|p)*\beta(p;a,b)\)

求解极值问题等价于求解log之后的极值问题。

a, b = symbols('a b')
beta = p**(a-1)*(1-p)**(b-1)
posterior = like*beta
posterior = log(posterior)
posterior

\(\displaystyle \log{\left(p^{k} p^{a - 1} \left(1 - p\right)^{b - 1} \left(1 - p\right)^{- k + n} \right)}\)

类似于1.2，我们来求导数。

posterior_prime = diff(posterior, p)
eq = Eq(posterior_prime, 0)
solutions:list = solve(eq, p) 
solutions

[0**(1/(a + k - 1)),
 (a + k - 1)/(a + b + n - 2),
 1 - 0**(1/(k - n)),
 1 - 0**(1/(b - k + n - 1))]

可以得到critical point为0, 1和 \(\frac{a + k - 1}{a + b + n - 2}\)，同样的可以发现第三个解比较通用，涵盖了前两个解的情况，所以贝叶斯估计的结果是\(p = \frac{a + k - 1}{a + b + n - 2}\)。如果a=b=1,那么和最大似然估计的结果一样。

第二题

通过经验风险最小化推导极大似然估计。证明模型是条件概率分布，当损

失函数是对数损失函数时，经验风险最小化等价于极大似然估计。

审题

• 题目出处？（题目中的概念是哪一本书的定义？）
  • 李航《统计学习方法》
• 李航书中“经验风险”和“对数损失函数”是怎么定义的？
  • 经验风险\(R_{emp}(f)=\frac{1}{N}\sum_{i=1}^{N}L(y_i, f(x_i))\)
  • 对数损失函数\(Loss(Y, P(Y|X)) = -logP(Y|X)\)

解题

根据题意，当损失函数是对数损失函数时，经验风险为

\(R_{emp}(f)=\frac{1}{N}\sum_{i=1}^{N}L(y_i, f(x_i))=\frac{1}{N}\sum_{i=1}^{N}-logP(y_i|x_i)\)

经验风险最小化是指

\(f = \mathop{\arg\min}_{f} R_{emp}(f) = \mathop{\arg\min}_{f} -\sum_{i=1}^{N}logP(y_i|x_i)=\mathop{\arg\max}_{f}\sum_{i=1}^{N}logP(y_i|x_i)\)

由于log函数的性质，

\(f = \mathop{\arg\max}_{f}\prod_{i=1}^{N}P(y_i|x_i) = \mathop{\arg\max}_{f}\prod_{i=1}^{N}\frac{P(y_i, x_i)}{P(x_i)}\)，由于evidence概率\(P(x_i)\)与模型的参数无关，可以认为是一个常数，所以在优化中可以去除，注意这是重要的一个推导。

\(f = \mathop{\arg\max}_{f}\prod_{i=1}^{N}P(y_i, x_i) = \mathop{\arg\max}_{f}\prod_{i=1}^{N}P(y_i, x_i| f)\)， 注意，在上文中李航书中的条件概率省略了f（而PRML书是没有省略的），比如上文的\(P(Y|X)\)实际上是\(P(Y|X, f)\)，这里我们为了和似然做比较，显式地加回来。

由于数据i.i.d.

\(f = \mathop{\arg\max}_{f}P(Y,X|f)\), 而这就是极大似然估计，给定一个f，求出现这个数据集中的(X, Y)的出现概率，这个叫做似然。因此得证。

附加题

考虑一个回归模型 $ f \(,它的目标变量为\)t=f(x,w,^2)+$，其中 \(\varepsilon\) 是一个随机噪声，\(\varepsilon\) 的概率密度为： \[p(x)=\frac{q}{2(2\sigma ^2)^{\frac{1}{q}}\Gamma (\frac{1}{q})}\exp(-\frac{|x|^q}{2\sigma ^2})\] 给定观测数据集 \(Data=\{(X,t)\}= \{(x_1,t_1)...(x_N,t_N)\}\),求f 关于参数 w 和 \(\sigma^2\) 的对数似然函数。

审题

• 题目出处？
  • 推测是PRML上的
• \(\Gamma\)是什么函数？
  • Γ函数 - 维基百科，自由的百科全书 (wikipedia.org)
  • 是阶乘函数的一个扩展

解题

本题我们不知道f的形式，但是题目告诉了我们ground truth，也就是t的真实情况是一个已知的f函数加上一个噪声，f函数的未知参数需要我们根据极大似然估计来求出。

似然函数本来是问我们，假如已知了\(w,\sigma^2\)，对应数据集Data出现的概率是多少，也就是\(P(Data|w,\sigma^2)\)。

如果有了\(w,\sigma^2\), f也知道，我们算出来的预测\(\hat{t}=f(x,w,\sigma^2)\)与t之间有差值，如果我们\(w,\sigma^2\)就是真相，那么对应的差值就是真正的\(\varepsilon\)，那么它就应该服从题目描述的那个分布，这个情况下，我们就能算出来\(P(\varepsilon)\)，这个\(P(\varepsilon)\)就是我们给定\(w,\sigma^2\)情况下数据出现的概率。

\(LogLike(w,\sigma^2) = logP(Data|w,\sigma^2) = \sum_{i=1}^{N}logP(\varepsilon_i|w,\sigma^2)) = \sum_{i=1}^{N}log p(t_i-f(x_i, w,\sigma^2); w,\sigma^2)\)

带入题目给的概率密度p，这个式子变为

$LogLike(w,^2) = _{i=1}^{N} ( (-) ) $

利用log的性质

\(LogLike(w,\sigma^2) = \sum_{i=1}^{N} \left( \log(q) - \log(2(2\sigma^2)^{\frac{1}{q}}\Gamma(\frac{1}{q})) - \frac{|t_i - f(x_i, w, \sigma^2)|^q}{2\sigma^2} \right)\)

最后一项和求和有关，其他都没有关系

\(LogLike(w,\sigma^2) = N \log(q) - N \log(2(2\sigma^2)^{\frac{1}{q}}\Gamma(\frac{1}{q})) - \frac{1}{2\sigma^2} \sum_{i=1}^{N} |t_i - f(x_i, w, \sigma^2)|^q\)

由于没法继续进行化简，这个式子就是我们给出的答案。

当我们做极大似然估计的时候，q和N是常数，所以损失函数为

\(Loss(w,\sigma^2) = \frac{1}{2\sigma^2} \sum_{i=1}^{N} |t_i - f(x_i, w, \sigma^2)|^q + N \log(2(2\sigma^2)^{\frac{1}{q}}\Gamma(\frac{1}{q}))\)

我们可以看到，这是用了误差绝对值的q次方来惩罚，然后后面增加了一个正则项防止\(\sigma^2\)的值太大。