机器学习小组知识点6:最大似然估计法(MLE)

为了照顾大家不同的基础,我只能写的不是很数学化,因此就用大家都能理解的东西。

最大似然估计中采样需满足一个很重要的假设,就是所有的采样都是独立同分布的。下面我们具体描述一下最大似然估计:
最大似然估计提供了一种“给定观察数据来评估模型参数”的方法,即:“模型已定,参数未知”。简单而言,假设我们要统计全国人口的身高,首先假设这个身高服从服从正态分布,但是该分布的均值与方差未知。我们没有人力与物力去统计全国每个人的身高,但是可以通过采样,获取部分人的身高,然后通过最大似然估计来获取上述假设中的正态分布的均值与方差。
备注:
1. 什么是独立?独立所谓的就是不受其他的影响,例如:上例中的人的身高不能相互影响,当然你不能来硬的说,姚明的女儿个子肯定很高,怎么可能独立?这里我们要这样想,样本太大了,那么他们爷俩的影响微不足道,毕竟姚明跟你跟我关系不大,身高上几乎不影响。
2. 什么是同分布?这个现实生活中就是身高的样本来自于中国这个区域,来自同一个地域,所以分布应该可以认为是来自一个分布。说到底,用我统计的同学话说,他们就是用同一个样本:中国人民的身高。

最大似然估计
首先,假设x_1,x_2,\cdots,x_n为独立同分布的采样,\theta为模型参数,f为我们使用的模型,遵循上面的独立分布假设。参数为\theta的模型f产生的上述采样为

    \[f(x_1,x_2,\cdots,x_n)=f(x_1|\theta)*f(x_2|\theta)*\cdots *f(x_n|\theta)\]

解释下为什么能这样乘?就是他们独立同分布! 不影响,可以拆开。
这样模型就固定了,就是上式。回到上面所说的“模型固定,参数未知”说法,此时我们已知的为x_1,x_2,\cdots,x_n(为什么?我们采样了,挨个挨个测量了他们的身高嘛),未知的为\theta(为什么?我不知道国人的平均身高和方差,所以要求。这里只是假设了个参数,有点儿像待定系数法),所以”极大似然估计”定义为:

    \[L(\theta|x_1,x_2,\cdots,x_n)=f(x_1,x_2,\cdots,x_n)=\Pi_{i=1}^nf(x_i|\theta)\]

在实际应用中对两边取对数,(为什么取对数?凸的,单调递增的函数我喜欢!),得到

    \[ln L(\theta|x_1,x_2,\cdots,x_n)=\sum_{i=1}^nln f(x_i|\theta)\]

于是有

    \[\hat{l}=\frac{1}{n}ln L\]

其中ln L(\theta|x_1,x_2,\cdots,x_n)为对数似然,而\hat{l}为平均对数似然。为我们平时所称的最大似然为最大的对数平均似然,即:

    \[\hat{\theta}=arg\ max\ \hat{l}(\theta|x_1,x_2,\cdots,x_n)\]

举个别人博客中的例子,假如有一个罐子,里面有黑白两种颜色的球,数目多少不知,两种颜色的比例也不知。我 们想知道罐中白球和黑球的比例,但我们不能把罐中的球全部拿出来数。现在我们可以每次任意从已经摇匀的罐中拿一个球出来,记录球的颜色,然后把拿出来的球 再放回罐中。这个过程可以重复,我们可以用记录的球的颜色来估计罐中黑白球的比例。假如在前面的一百次重复记录中,有七十次是白球,请问罐中白球所占的比例最有可能是多少?很多人马上就有答案了:70%。而其后的理论支撑是什么呢?

我们假设罐中白球的比例是p,那么黑球的比例就是1-p。因为每抽一个球出来,在记录颜色之后,我们把抽出的球放回了罐中并摇匀,所以每次抽出来的球的颜 色服从同一独立分布。这里我们把一次抽出来球的颜色称为一次抽样。题目中在一百次抽样中,七十次是白球的概率是P(Data | M),这里Data是所有的数据,M是所给出的模型,表示每次抽出来的球是白色的概率为p。如果第一抽样的结果记为x1,第二抽样的结果记为x_2\cdots,那么Data = (x_1,x_2,…,x_{100})。这样,
P(Data | M)= P(x1,x2,…,x_{100}|M)= P(x_1|M)P(x_2|M)\cdots P(x_{100}|M)= p^{70}(1-p)^{30}.
那么p在取什么值的时候,P(Data |M)的值最大呢?将p^{70}(1-p)^{30}p求导,并其等于零。
解得p=0.7
在边界点p=0,1P(Data|M)=0。所以当p=0.7时,P(Data|M)的值最大。这和我们常识中按抽样中的比例来计算的结果是一样的。

由上可知最大似然估计的一般求解过程:

(1) 写出似然函数;

(2) 对似然函数取对数,并整理;

(3) 求导数 ;

(4) 解似然方程
注意:
最大似然估计函数不一定是唯一的,也可能不存在。

总结:
极大似然估计(MLE)是建立在**极大似然原理上**的参数估计,在假定**总体服从某种分布**的情况下,我们可以通过一些样本,来估计总体的参数。

这就给我们抛出来几个问题,什么是极大似然原理,什么是参数估计,又是如何进行极大似然估计的。

1.极大似然原理
极大似然原理,简单地描述就是“存在即合理,所见即真实”,是人们对于事物认知的一种假设。**对于一个有多种可能的事件,我们总是认为,我们看到的那些结果(样本),是该事件概率最大的可能**。
比如,对于外星生物,人类最会理所应当地假设,我们见到的第一个外星生物的长相是普遍外形生物的模样。
更多的原理解释可以看[维基百科](https://zh.wikipedia.org/wiki/%E6%9C%80%E5%A4%A7%E4%BC%BC%E7%84%B6%E4%BC%B0%E8%AE%A1)。

2.参数估计
在统计里面,我们习惯于用一些参数来描述事物整体的分布情况。比如“男性平均寿命74岁”,这里面的74便是**描述整体模样的一个参数**。
这就是我们上面说的\theta\sigma
而通过我们取样的一千个男性的平均寿命来估计全体男性的平均寿命,这样的过程就称为参数估计。

机器学习小组知识点4&5:批量梯度下降法(BGD)和随机梯度下降法(SGD)的代码实现Matlab版

这里趁着脑子还清醒就把代码敲出来了:
亲测是没有bug的,欢迎大家来找虫。


clear all
close all
clc

%% initialization
% input x & y
x1 = [1 3 2104];        y1 = 400;
x2 = [1 3 1600];        y2 = 330;
x3 = [1 3 2400];        y3 = 369;
x4 = [1 2 1416];        y4 = 232;
x5 = [1 4 3000];        y5 = 540;

X = [x1;x2;x3;x4;x5];
Y = [y1;y2;y3;y4;y5];

N = 5;      % the number of sample

% parameter
theta_new = zeros(length(x1),1);

% iteration number
itn = 2e2;

% tolerance
tol = 1e-20;

% learning rate(stepsize)
alpha_BGD = 0.00000001;
alpha_SGD = 0.0000001;
%% main function

i = 1;       %itn index
theta_old = theta_new; 
error = 1.0;
k = 1; % SGD sample index

enenrgy = zeros(itn,1);

option = input('Enter the case(number) you want: \n 1 -> BGD. \n 2 -> SGD.\n ' )

while((i <= itn) && (error >= tol))

    switch option        

        case 1            
            % BGD
            sum_temp = 0;
                for jj = 1:N
                    sum_temp = sum_temp + ( Y(jj) - X(jj,:)*theta_old )*X(jj,:)';
                end
            theta_new = theta_old + alpha_BGD*sum_temp;

        case 2
            % SGD            
            theta_new = theta_old + alpha_SGD*( Y(k) - X(k,:)*theta_old )*X(k,:)';  

            if k == 5
                k = 1;
            else
                k = k+1;
            end

       end


    error = norm(X*theta_new-Y)^2;
    theta_old = theta_new;

    energy(i) = error;


    i = i + 1;
end

figure,plot(1:length(energy),energy),title('error(energy)');
itn = i-1
theta = theta_new

如果输出图的话,大家会看到SGD误差振荡下降。说明了SGD的噪音太多然而BGD却一路向下。

机器学习小组知识点3:最小二乘法(LSM)

上篇博客介绍了最小均方算法(LMS),其实里面的东西包含的很多,其中有最小二乘法,梯度下降以及随机梯度下降法。这篇博客着重介绍最小二乘法的推导,来源以及做一点儿推广。下面进入正题:

最小二乘法的闭形式推导
在上篇博客我们引入了J(\theta)成本函数的具体形式,这里我们要推导出关于\theta的“闭形式”,数学上也称为解析解的形式。下面我们要重新将J写成矩阵乘向量的形式。

给定一个训练集,定义“设计矩阵”X是一个m*n的矩阵(实际上就是样本输入变量写成的矩阵,下面就看到了),这里要注意实际上是m*(n+1)维矩阵,至于为什么n+1维矩阵呢?这里要回头看看上篇博客里面的截矩项,他多占了其中的一个维度。不过为了叙述方便,我们还是统一为n维。
我们先给出X的矩阵形式,

    \[\left[ \begin{matrix} x^{(1)^T} \\ x^{(2)^T}\\ \vdots \\ x^{(m)^T} \end{matrix} \right]\]

这里我们还是保留了老传统,样本数量是m个,特征的维度总数为n维。另外,里面的x都是列向量。
接下来,相对应的就是数据项,也就是我们样本已经有的观测数据项,再清楚点儿就是我们监督学习里面起到“监督”二字的关键,那么他本身应该是对应x的维度的,那么这样的话,我们就能够得到m维的向量y了,见下面形式Y为:

    \[\left[ \begin{matrix} y^{(1)} \\ y^{(2)}\\ \vdots \\ y^{(m)} \end{matrix} \right]\]

现在因为h_\theta(x^{(i)})=(x^{(i)})^T\theta,那么我们简单的写出数据拟合的形式即

    \[X\theta-Y= {\left[ \begin{matrix} x^{(1)^T}\theta \\ x^{(2)^T}\theta\\ \vdots \\ x^{(m)^T}\theta \end{matrix} \right]}- \left[ \begin{matrix} y^{(1)} \\ y^{(2)}\\ \vdots \\ y^{(m)} \end{matrix} \right]\]

    \[=\left[   \begin{matrix}    h_\theta(x^{(1)})-y^{(1)}\\    h_\theta(x^{(2)})-y^{(2)}\\    \vdots \\    h_\theta(x^{(m)})-y^{(m)}   \end{matrix}    \right]\]

因此,利用一个事实,“对于一个向量z,我们有z^Tz=\sum_iz^2_i”,那么我们就利用这个事实可以得到:

    \[\frac{1}{2}(X\theta-Y)^T(X\theta-Y)=\frac{1}{2}\sum_{i=1}^m (h_\theta(x^{(i)})-y^{(i)})^2\]

    \[=J(\theta)\]

因此我们对\theta进行求导可以得到**正规方程**,与上篇博客相互照应:

    \[X^TX\theta=X^TY\]

所以我们就能得到关于\theta的闭形式(解析形式)为

    \[\theta=(X^TX)^{-1}X^TY\]

这里注意,需要矩阵求导的相关知识,如果大家不熟悉,强烈要求讲的话,我再敲一个专题,不过同学们看文献[1]就应该知道了。

最小二乘法概率解释(来源)
这一部分解释的就是为啥是最小二乘法,而不是最小一乘法呢?要知道最小二乘法这些东西都是当年统计的人们天天搞的呀,那跟统计有何关系呢?这里默认为高斯老人家也可以称之为统计学家了。
首先让我们假设目标变量y和输入变量x是满足线性的关系(照应我们上篇博客),满足以下的方程:

    \[y^{(i)}=\theta^{T}x^{(i)}+\epsilon^{(i)}\]

其中这个\epsilon^{(i)}来表示误差项,或者是随机的噪音。图像处理的人看这个“一般”就是默认高斯噪音了。那么为什么是高斯噪音呢?其实高斯噪音适用的范围是有很多独立的因素对模型造成的误差导致的,生活中很多场所一般用的是高斯噪音。与之区分的是椒盐噪音,拉普拉斯噪音以及莱斯噪音(肖凯博士告诉我的)。有点儿扯远了,下面我们继续假设\epsilon^{(i)}来源于高斯分布(正态分布),也就是均值为0,方差为\sigma^2,且噪音对应没个样本是相互独立的。
这里多说两点:

1. 用高斯噪音,大家只是默认了很多情况。方差其实不一定都是同一个,现在很多研究人员就在搞混合噪音,方差也是混合的阶段,但是就目前来看,也只是停留在社会初级阶段,还不是一个热门。
2. 为啥均值是0啊,吕大忽悠?我没骗你,因为可以把均值放到前面数据拟合项的截矩呢,你说是不是?
回归正题,这里我们就假设\epsilon^{(i)}服从N(0,\sigma^2)的正态分布。那么对应的概率密度函数是:

    \[p(\epsilon{(i)})=\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(\epsilon^{(i)^2})}{2\sigma^2})\]

那也就说明

    \[p(y^{(i)}|x^{(i)};\theta)=\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(y^{(i)}-\theta^Tx^{(i)})^2}{2\sigma^2})\]

这里p(y^{(i)}|x^{(i)};\theta)表示y^{(i)}的分布是由x^{(i)}和参数\theta共同决定的。
上面是一个样本y^{(i)}的分布,那么问题来了,那一堆样本,m个到底怎么办啊?好在他们都是相互独立的,又是同分布的(IID),就是来源于同一个样本空间(娘胎),这里我们就要用到极大似然估计方法(这里抛个砖,下篇文章介绍),

(1)   \[L(\theta)=\Pi^{m}_{i=1}p(y^{(i)}|x^{(i)};\theta)=\Pi^{m}_{i=1}\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(y^{(i)}-\theta^Tx^{(i)})^2}{2\sigma^2})\]

现在我们给了这个方程,接下来就是最大化(1)式了,这里还是极大似然估计的用法。
由于(1)式是联乘,那么自然反应就是用log函数了,为啥?log是凸的又是单调递增的不影响原来函数的凸性质和递增的性质呀!
接下来就是推导过程:

    \[l(\theta)=logL(\theta)\]

    \[=log\Pi_{i=1}^m\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(y^{(i)}-\theta^Tx^{(i)})^2}{2\sigma^2})\]

    \[=\sum_{i=1}^mlog\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(y^{(i)}-\theta^Tx^{(i)})^2}{2\sigma^2})\]

    \[=mlog\frac{1}{\sqrt(2\pi)\sigma}-\frac{1}{2\sigma^2}\sum_{i=1}^m(y^{(i)}-\theta^Tx^{(i)})^2\]

因此,我们最大化l(\theta)得到同样的结果就是最小化下面的式子

    \[\frac{1}{2}\sum_{i=1}^m(y^{(i)}-\theta^Tx^{(i)})^2\]

这是什么啊?对!没错,就是J(\theta)呀!这里我们终于解释清楚了,为啥要用最小二乘了吧?
总结:
1. 最小二乘不能瞎用,他针对的是高斯噪音影响的数据,大多数出现的情况还是高斯噪音,就是很多因素同时是独立的影响的时候,我们就用它。不能细说,细说就得去看统计上的高斯分布了。
2. 最小二乘法是来源于极大似然估计方法,其实两者是可以互推的。所以,是等价的。
3. \theta你这里是固定的值啊,为哈不能看成依赖于\sigma变量?吕大忽悠解释下,这个还真不能细细地去挖,因为现在他们统计部门还在争执呢!有两个学派,一个是看成常数,另外一派把他们看成以来\sigma随机变量。但是备注下,Andrew Ng提到对于我们这个问题,这两派得到的结果是一样的。我曾经深挖过,现在还是忘了!感兴趣的可以看看知乎的大神们的详细解释。

———-
参考文献:
[1] http://cs229.stanford.edu/notes/cs229-notes1.pdf

机器学习小组知识点2:最小均方算法(LMS)

感谢:向原文作者Nanshu Wang致以崇高的敬意!(http://nanshu.wang/post/2015-02-10/)
声明: 博主只是推了一遍然后细节更多的阐释了下,再次向原文作者致谢!

  1. 有监督学习(Supervised Learning)
  2. 线性回归(Linear Regression)
  3.  LMS算法
  4. 正规方程(The Normal Equations)
  5. 正规方程与梯度下降的比较

 

  1. 有监督学习
    先理清几个概念:
    1. x^i: x^i表示”输入”变量(“input” variables),也称为特征(features)。
    2. y^i: y^i表示”输出”变量(“output” variables),也称为目标值(target)。
    3. 一对(x^i,y^i)称为一个训练样本(training example),用作训练的数据集就是就是一组m个训练样本(x^i,y^i);i=1,…,m;被称为训练集(training set)。
    4. X表示输入变量的取值空间,Y表示输出变量的取值空间。那么h:X→Y是训练得到的映射函数,对于每个取值空间X的取值,都能给出取值空间Y上的一个预测值。函数hh的含义为假设(hypothesis)。
    5. 当预测值y为连续值时,则有监督学习问题是回归(regression)问题;预测值y为离散值时,则为分类(classification)问题。 图形化表示整个过程:20161015104654923
  2. 线性回归
    先简单将y表示为x的线性函数:

        \[h(x)=\sum^n_{i=0}\theta_ix_i=\theta^Tx\]

    1. \theta称为参数(parameters),也叫做权重(weights),参数决定了XY的射映空间。
    2. 用x_0=1来表示截距项(intercept term)。 有了训练集,如果通过学习得到参数\theta? 一种方法是,让预测值h(x)h(x)尽量接近真实值y,定义成本函数(cost function):

        \[J(\theta)=\frac{1}{2}\sum^m_{i=1}(h_\theta(x_i)-y_i)^2\]

  3. LMS算法 为了找到使成本函数J(\theta)最小的参数\theta,采用搜索算法:给定一个\theta 的初值,然后不断改进,每次改进都使J(\theta)更小,直到最小化J(\theta)\theta的值收敛。 考虑梯度下降(gradient descent)算法:从初始\theta开始,不断更新:

        \[\theta_j:=\theta_j-\alpha\frac{\partial}{\partial\theta_j}J(\theta)\]

    注意,更新是同时对所有j=0,…,n\theta_j值进行。\alpha被称作学习率(learning rate),也是梯度下降的长度,若\alpha取值较小,则收敛的时间较长;相反,若\alpha取值较大,则可能错过最优值。 算法具体的细节请见下面图片:(因为我实在不想敲公式了。。。。)

20161015111329598

BTW,是的,我手动加滤镜了。向我女票致敬!

图片内容补充:
1. 第一种“整体上”的梯度下降算法就是“批量梯度下降算法”(batch gradient descent)。
2. 第二种“局部上”的梯度下降算法就是“随机梯度下降算法”(stochastic gradient descent),这是因为选择下降的样本只是一部分,**更多的参考资料显示仅仅一个,并且就是对应的样本**。
3. 随机梯度下降和批量梯度下降不同点在于,批量梯度下降每一步更新\theta值,都需要遍历全部的训练样本,而随机梯度下降在遇到每个训练样本时,更新\theta之后继续处理下一个样本,每个样本只遍历一次,算法的学习时间比批量梯度下降快很多。但是,随机梯度下降可能永远不会收敛到全局最优值,而是在成本(惩罚)函数J(\theta)最优值周围附近摇摆。但是在实际问题中,接近最优值的参数值可能已经是足够好的结果了,特别是对于数据量非常大的训练集来说,随机梯度下降是比批量梯度下降更好的选择。
4. 在实际使用梯度下降算法时,将输入变量归一到同一取值范围,能够减少算法的迭代次数。这是因为\theta在小的取值范围内会下降很快,但在大的取值范围内会下降较慢。并且当输入变量取值范围不够均衡时,\theta更容易在最优值周围波动。采用特征缩放(feature scaling)和均值归一化(mean normalization)可以避免这些问题。选择学习率\alpha的值,要观察J(\theta)值在每次迭代后的变化。已经证明了如果\alpha的取值足够小,则J(\theta)每次迭代后的值都会减少。如果J(\theta)在某次迭代后反而增加了,说明学习率\alpha的值应该减小,因为错过了最优值。Andrew Ng推荐的一个经验是每次将\alpha减少3倍。
5. 这里想看一个样本的推导,可见(http://nanshu.wang/post/2015-02-10/)。

D.正规方程
梯度下降是最小化J(\theta)的一种方式,正规方程是另一种求解参数\theta的方法,这种方法可以直接求出最优值参数结果,不需要迭代更新,也不需要事先对数据进行归一化预处理。这种方法实际上是直接求出J(\theta)的导数,并令其为0。 **注明:明白最小二乘的同学,一看就知道什么东西了。** ! 20161015113629410

正规方程求解\theta的时间复杂度为O(n^3)n是特征数量。当特征数量很大时,正规方程求解会很慢。Andrew Ng给出的一个经验参考是:当n>10,000时,采用梯度下降比正规方程更好。

5.正规方程与梯度下降的比较

正规方程 梯度下降法
不需要调整参数 不需要调整参数α
不需要迭代 需要迭代更新权重θ
n n

注意,正规方程还会遇到X^TX不可逆的情况,通常这是因为输入变量中存在线性相关的变量或者是因为特征太多(nm)。解决方法是去掉线性相关的冗余变量,或者删掉一些特征。这里就是数据分析里面的回归系数的假设性检验了。感兴趣的同学,可以参考《线性统计模型》这本书。
更多的细节,会在小组讨论中提出来。