世界上有哪些代码量很少，但很牛逼很经典的算法或项目案例？

前面的回答有提到FFT，这个算法也很牛逼，另外我觉得的PID算法也很强，特别是在工业控制中，使用非常广泛，下面是我之前学习PID写的笔记，可以参考一下；

1 前言

控制系统通常根据有没有反馈会分为开环系统和闭环系统，在闭环系统的控制中，PID算法非常强大，其三个部分分别为；

• P：比例环节；
• I：积分环节；
• D：微分环节；

PID算法可以自动对控制系统进行准确且迅速的校正，因此被广泛地应用于工业控制系统。

2 开环控制

首先来看开环控制系统，如下图所示，隆哥蒙着眼，需要走到虚线旗帜所表示的目标位置，由于缺少反馈（眼睛可以感知当前距离和位置，由于眼睛被蒙上没有反馈，所以这也是一个开环系统），最终隆哥会较大概率偏离预期的目标，可能会运行到途中实线旗帜所表示的位置。

开环系统的整体结构如下所示；

这里做一个不是很恰当的比喻；

• Input：告诉隆哥目标距离的直线位置（10米）；
• Controller：隆哥大脑中计算出到达目标所需要走多少步；
• Process：双腿作为执行机构，输出了相应的步数，但是最终仍然偏离了目标；

看来没有反馈的存在，很难准确到达目标位置。

3 闭环控制

所以为了准确到达目标位置，这里就需要引入反馈，具体如下图所示；

在这里继续举个不怎么恰当的比喻；隆哥重获光明之后，基本可以看到目标位置了；

• 第一步Input：告诉隆哥目标距离的直线位置（10米）；
• 第二步Controller：隆哥大脑中计算出到达目标所需要走多少步；
• 第三步Process：双腿作为执行机构，输出了相应的步数，但是最终仍然偏离了目标；
• 第四步Feedback：通过视觉获取到目前已经前进的距离，（比如前进了2米，那么还有8米的偏差）；
• 第五步err：根据偏差重新计算所需要的步数，然后重复上述四个步骤，最终隆哥达到最终的目标位置。

4 PID

4.1 系统架构

虽然在反馈系统下，隆哥最终到达目标位置，但是现在又来了新的任务，就是又快又准地到达目标位置。所以这里隆哥开始采用PID Controller，只要适当调整P，I和D的参数，就可以到达目标位置，具体如下图所示；

隆哥为了最短时间内到达目标位置，进行了不断的尝试，分别出现了以下几种情况；

• 跑得太快，最终导致冲过了目标位置还得往回跑；
• 跑得太慢，最终导致到达目标位置所用时间太长；

经过不断的尝试，终于找到了最佳的方式，其过程大概如下图所示；

这里依然举一个不是很恰当的比喻；

• 第一步：得到与目标位置的距离偏差（比如最开始是10米，后面会逐渐变小）；
• 第二步：根据误差，预估需要多少速度，如何估算呢，看下面几步；

P比例则是给定一个速度的大致范围，满足下面这个公式；

因此比例作用相当于某一时刻的偏差（err）与比例系数的乘积，具体如下所示；

绿色线为上述例子中从初始位置到目标位置的距离变化； 红色线为上述例子中从初始位置到目标位置的偏差变化，两者为互补的关系；

I积分则是误差在一定时间内的和，满足以下公式；

如下图所示；

红色曲线阴影部分面积即为积分作用的结果，其不断累积的误差，最终乘以积分系数就得到了积分部分的输出；

D微分则是误差变化曲线某处的导数，或者说是某一点的斜率，因此这里需要引入微分；

从图中可知，当偏差变化过快，微分环节会输出较大的负数，作为抑制输出继续上升，从而抑制过冲。

综上，，分别增加其中一项参数会对系统造成的影响总结如下表所示；

4.2 理论基础

上面扯了这么多，无非是为了初步理解PID在负反馈系统中的调节作用，下面开始推导一下算法实现的具体过程；PID控制器的系统框图如下所示；

因此不难得出输入和输出的关系；

是比例增益； 是积分增益； 是微分增益；

4.3 离散化

在数字系统中进行PID算法控制，需要对上述算法进行离散化；假设系统采样时间为 则将输入序列化得到；

将输出序列化得到；

• 比例项：
• 积分项：
• 微分项：

所以最终可以得到式①，也就是网上所说的位置式PID：

将式①再做一下简化；

最终得到增量式PID的离散公式如下：

4.4 伪算法

这里简单总结一下位置式PID实现的伪算法；

       previous_error := 0  //上一次偏差 integral := 0   //积分和  //循环  //采样周期为dt loop:  //setpoint 设定值  //measured_value 反馈值     error := setpoint − measured_value //计算得到偏差     integral := integral + error × dt //计算得到积分累加和     derivative := (error − previous_error) / dt //计算得到微分     output := Kp × error + Ki × integral + Kd × derivative //计算得到PID输出     previous_error := error //保存当前偏差为下一次采样时所需要的历史偏差     wait(dt) //等待下一次采用     goto loop        

5 C++实现

这里是位置式PID算法的C语言实现；

pid.cpp

       #ifndef _PID_SOURCE_ #define _PID_SOURCE_  #include <iostream> #include <cmath> #include "pid.h"  using namespace std;  class PIDImpl {     public:         PIDImpl( double dt, double max, double min, double Kp, double Kd, double Ki );         ~PIDImpl();         double calculate( double setpoint, double pv );      private:         double _dt;         double _max;         double _min;         double _Kp;         double _Kd;         double _Ki;         double _pre_error;         double _integral; };   PID::PID( double dt, double max, double min, double Kp, double Kd, double Ki ) {     pimpl = new PIDImpl(dt,max,min,Kp,Kd,Ki); } double PID::calculate( double setpoint, double pv ) {     return pimpl->calculate(setpoint,pv); } PID::~PID()  {     delete pimpl; }   /**  * Implementation  */ PIDImpl::PIDImpl( double dt, double max, double min, double Kp, double Kd, double Ki ) :     _dt(dt),     _max(max),     _min(min),     _Kp(Kp),     _Kd(Kd),     _Ki(Ki),     _pre_error(0),     _integral(0) { }  double PIDImpl::calculate( double setpoint, double pv ) {          // Calculate error     double error = setpoint - pv;      // Proportional term     double Pout = _Kp * error;      // Integral term     _integral += error * _dt;     double Iout = _Ki * _integral;      // Derivative term     double derivative = (error - _pre_error) / _dt;     double Dout = _Kd * derivative;      // Calculate total output     double output = Pout + Iout + Dout;      // Restrict to max/min     if( output > _max )         output = _max;     else if( output < _min )         output = _min;      // Save error to previous error     _pre_error = error;      return output; }  PIDImpl::~PIDImpl() { }  #endif       

pid.h

       #ifndef _PID_H_ #define _PID_H_  class PIDImpl; class PID {     public:         // Kp -  proportional gain         // Ki -  Integral gain         // Kd -  derivative gain         // dt -  loop interval time         // max - maximum value of manipulated variable         // min - minimum value of manipulated variable         PID( double dt, double max, double min, double Kp, double Kd, double Ki );          // Returns the manipulated variable given a setpoint and current process value         double calculate( double setpoint, double pv );         ~PID();      private:         PIDImpl *pimpl; };  #endif       

pid_example.cpp

       #include "pid.h" #include <stdio.h>  int main() {      PID pid = PID(0.1, 100, -100, 0.1, 0.01, 0.5);      double val = 20;     for (int i = 0; i < 100; i++) {         double inc = pid.calculate(0, val);         printf("val:% 7.3f inc:% 7.3f
", val, inc);         val += inc;     }      return 0; }       

编译并测试；

       g++ -c pid.cpp -o pid.o # To compile example code: g++ pid_example.cpp pid.o -o pid_example       

6 总结

本文总结了PID控制器算法在闭环系统中根据偏差变化的具体调节作用，每个环节可能对系统输出造成什么样的变化，给出了位置式和离散PID算法的推导过程，并给出了位置式算法的C++程序实现。文章中只是最简易的PID算法，并没做优化，在实际的应用中，还需要做积分抗饱和，微分环节的优化，具体需要参考实际系统的情况来定。

由于作者能力和水平有限，文中难免存在错误和纰漏，请不吝赐教。

手动码字，感觉走心的话，点赞，关注支持一下吧； @极客小麦

快速傅里叶变换（FFT），核心算法用递归法几十行搞定。

计算机应用最为重要的算法之一，广泛应用于各种各样的工业和科学实践中。第一次看到递归算出快速傅里叶的时候，惊了很久。后来经过数学家的很多努力，有了新算法将复杂度从N^2，降到N*log(N)，觉得可以称之为人类智慧的精华了。

有空可以贴下源代码。

2020.1.15 科研狗这几天忙着科研，大家别急，今天一定抽空把算法和源代码贴出来。

第一次更新（1.15）：

众所周知，数学中有一个傅里叶变换，是将函数在实空间的定义域变换到倒空间的定义域上，这个变换非常重要，因为一些实空间不易发现的信息，在倒空间中会表现的非常明显，所以傅里叶变换是很多信号、物理的分析基础。它的形式如下：

以及其逆变换：

但是如果在计算机想要编程实现傅里叶变换就要不可避免地面临两个近似：第一，现实中我们无法将积分的上下界延伸到无穷，所以这里一定是一个有界积分；第二，无论什么连续的公式，只要是想要被计算机实现，就必须进行离散化，所以这里的连续积分（ ）就近似为离散求和（ ）。而基于这两种朴素近似下的，就是所谓的“离散傅里叶变换”（Discrete Fourier Transform，简称DFT），它的形式也很直观：

看到这里，其实你如果利用DFT简单地编程，其实已经可以实现傅里叶变换了。但是我们为什么不用这个算法呢？是因为复杂度。

复杂度是算法中的一个概念，简单来说就是复杂度越高，那么这个算法计算大体系就越困难，那有多困难呢？我举个例子，现在算法基础库中对于矩阵的对角化应该是最基础和重要的了，矩阵的对角化的计算复杂度是 ，就是说对角化的时间就随着矩阵维数的三次方增加，这是十分恐怖的增长。目前市面的通用的对角化库，比如 lapcak，对角化的极限维数一般也就在 10000, 这个量级。那DFT的复杂度是多少呢？你可以从公式上看到，如果你现在要对一个包含N个数的数组进行DFT，那你每算其中一个变换后的数，都要进行N次的运算，那这个算法的复杂度也就是 ，也就是说，用DFT最多也就只能处理大概一个拥有大概 个数的数组。但是在生产科研实践中，这实在不是一个大的上限，所以这个复杂度的硬伤一直制约着DFT的发展，而这就是属于当时制约整个人类文明进步的“卡脖子”的难题，如果FFT没有出现，那我们的世界现在绝对不是现在这个样子。

时间来到1965年，James Cooley 和 John Tukey 发现了一种算法，这种的算法的复杂度降低为了 ，而这种算法就称之为“快速傅里叶变换”（Fast Fourier transform， 简称FFT）。但其实，就这个问题，Danielson 和 Lanczos 在1942年就已经开始了相关的讨论。他们发现任何一个DFT其实可以重新改写成两个DFT的和：

数学家进一步发现这一分为二的DFT组，一个是原来是偶数项（even，第0,2,4,6...项）组成的，而另一个是由原来的奇数项（odd）组成的。如果你只是看到这一步，你已经发现了一个大事情，就是原来用这个公式来算的话，求一个傅里叶变换不再需要进行N次相乘，而只需要odd组的N/2次运算，even组的直接不用做任何运算直接拉下来就好了。那这样计算量不就减少一半了吗？但是其实事情并没有那么简单，我们刚才证明了任何一个DFT组都能被一分为二，而且计算量减半。但是我们现在一分为二得到的even组和odd组，他们不又是DFT组吗？我们可以继续之前的操作，将他们分别再次分为两个数组的求和！

一分为二，二分为四，四分为八，一直分到每个数组只剩一个数字为止！而且每进行一次操作，就能使得其中一半的计算量被减少。这样我们就把计算量从 降到了 。你可能现在对这个算法的突破性还没有概念，我现在来给你构建下概念，

随着N的增加，FFT算法对于DFT的优势会不断显现，一开始1024个数，FFT比DFT快100倍；32768个数，就快了2000多倍，到了我们之前提到的极限1百万个数，快了近5万倍！FFT将我们的计算分析能力成万倍的提升。提一句，现在FFT最稳定最快的库是FFTW，这个W是“in West”，就是说“西方的FFT”。我当时刚看懂FFT时，不服气，一心想写个东方版，FFTE（FFT in East），后来和我自己的程序和FFTW性能一比，我人都傻了，果然好的FFT程序还是超级难写的。

这个算法关键的关键就是将DFT组不断奇偶细分，细分到最后时，如何确定每个数组前面的相位系数。我希望今天可以在第二次更新中把如何确定细分后的序列，以及蝶形算法更完。

我先把我当时研一第一次接触到FFT编写的Fortran源码贴出来吧。你品，细品。短短几十行，如果你能花几天的时间理解这个递归程序，那么这将是2020年获得的第一个受益终身的知识。

       ! Cooley-Tukey FFT recursive subroutine fft(x,sgn)      implicit none     integer, intent(in) :: sgn     complex(8), dimension(:), intent(inout) :: x     complex(8) :: t     integer :: N     integer :: i     complex(8), dimension(:), allocatable :: even, odd      N=size(x)      if(N .le. 1) return      allocate(odd((N+1)/2))     allocate(even(N/2))      ! divide     odd =x(1:N:2)     even=x(2:N:2)      ! conquer     call fft(odd, sgn)     call fft(even, sgn)      ! combine     do i=1,N/2         t=exp(cmplx(0.0d0,sgn*2.0d0*pi*(i-1)/N))*even(i)         x(i)     = odd(i) + t         x(i+N/2) = odd(i) - t     end do      deallocate(odd)     deallocate(even)  end subroutine fft     

有的，比如 洗牌算法：这个算法真的很牛逼很经典，而且代码很少。

评论区好多大佬，膜拜ing。。

先来思考一个问题：有一个大小为 100 的数组，里面的元素是从 1 到 100 按顺序排列，怎样随机的从里面选择 1 个数？

最简单的方法是利用系统的方法 Math.random() * 100 ，这样就可以拿到一个 0 到 99 的随机数，然后去数组找对应的位置就即可。

接下来在思考一个问题： 有一个大小为100的数组，里面的元素是从 1 到 100 按顺序排列，怎样随机的从里面选择 50 个数？

注意数字不能重复！

注意数字不能重复！

注意数字不能重复！

如果根据上面的思路，你第一想法是：随机 50 次不就行了？

但是，这样做有个很明显的 bug ：数字是会重复的。

修改一下？

弄一个数组，把每一次随机的数都放到数组里，下一次随机就看这个数组里面有没有这数，有的话就继续随机，直到这个数组里面有 50 个数字就停止。

这样是可以的！

但，还是有个小问题，考虑一下极端情况：有一个大小为100的数组，里面的元素是从 1 到 100 按顺序排列，怎样随机的从里面选择 99 个数。

如果按照上面的方法操作，越往后选择的数字跟前面已经挑选的数字重复的概率越高，这就会造成如果数组很大，选择的数字数目也很大的话，重复次数在量级上会很大。

这个时候就需要换一个思路，如果先将数组里面的元素打乱，那么按顺序选择前 50 个不就可以了？

是的！

但我们得注意什么叫乱？

一副扑克有 54 张牌，有 54! 种排列方式。所谓的打乱指的是，你所执行的操作，应该能够 等概率地生成 这 54! 种结果中的一种。

洗牌算法就能做到这一点。

洗牌算法

Fisher–Yates shuffle 算法由 Ronald Fisher 和 Frank Yates 于 1938 年提出，在 1964 年由 Richard Durstenfeld 改编为适用于电脑编程的版本。

这个算法很牛逼却很好理解，通俗的解释就是：将最后一个数和前面任意 n-1 个数中的一个数进行交换，然后倒数第二个数和前面任意 n-2 个数中的一个数进行交换。。。

小程序实现代码

       for (var i = this.rowCount * this.colCount - 1; i >= 0 ; i--){   var iX = parseInt(i / this.colCount);   var iY = i % this.colCount;    var randNumber = this.rangeRandom(0, i + 1);    var randX = parseInt(randNumber / this.colCount);   var randY = randNumber % this.colCount;   //交换两个位置   var temp = tmpMineMap[iX][iY];   tmpMineMap[iX][iY] = tmpMineMap[randX][randY];   tmpMineMap[randX][randY] = temp; }     

在整个过程中，这个算法保证了每一个元素出现在每一个位置的概率是相等的。

这个算法真的很牛逼很经典，而且代码很少。

我做了一个视频进行辅助理解。

为了避免知乎大佬觉得我吹逼，先贴一下自己的 GitHub 地址，目前 70，000 star，全球排名 51 名。

https:// github.com/MisterBooo

算法是一种技能，是可以通过科学合理的方式训练出来的能力。

在想刷题之前，得从心里认识到接受刷题很重要，才能坚持去刷题。

江湖有个传言：国内刷 LeetCode，最多够你吃 1 年老本；湾区刷 LeetCode ，够你吃 10 年老本了。

为什么湾区的刷题性价比这么高呢？

你想想，电面考 4 道题，一道题值 5 万！单位是 Dollar ！

刷到就是赚到！！

想想是不是很刺激，有没有动力开始刷题了！可以提速刷题了！

就目前互联网的情况来说，无论是面国外大厂还是面国内大厂，如果想换工作都要去刷题，一面二面不丢你几道 Hard 题，都对不住你偷偷摸摸找个会议室假装开会实则面试的鸡贼。

同时，还得认识到一点，面试能力和你平时的工作能力其实差别挺大的。

有些人技术挺厉害的，但没有刷题，一面二面都过不了，而某些小镇刷题家，还真就靠刷题拿下了 Google、微软、脸书等大厂offer。

国内大厂也有这种趋势，比如字节，一大半都是面试题。

要不是他提前先看视频刷题，妥妥得凉凉。

所以，刷题很重要。

（PS：感谢大家耐心的阅读，算法是程序员的重中之重，必须攻克，大厂面试必考，顺便送一份阿里大佬刷Leetcode总结的算法笔记，如果你能吃透，那我相信80%的技术面试都会不在话下：

BAT大佬写的Leetcode刷题笔记，看完秒杀80%的算法题！

这本书的目录，非常经典：

刷题大概可以分为 4 个阶段。

1、纯小白，不知道怎么刷题，对很多概念都很陌生，各种数据结构和知识点几乎完全不懂，打开 LeetCode 第一题，满头问号。

有人相爱、有人夜里开车看海、有人 LeetCode 第一题都做不出来。

2、算法上基本已经入门，Easy 可以做出来，Medium 纠结半天也能有头绪，但基础不牢，比如字符转字符串还得 Google 一下。

3、刷了几百道题后，总结了自己的解题模板，参加周赛有时候甚至可以全部完成。

4、开始以 beat 100% 作为 AC 的目标了。

就目前的算法面试大环境来说，能达到第二阶段，中小公司可以应付过去了，到达第三阶段，字节、腾讯算法面试环节妥妥没问题了。

怎么样到达第三阶段？

给一下我的一些小建议吧。

1、如果目标是国内大厂，那么一定要刷足够的题，不需要把 LeetCode 上 2500 道算法题都刷完，但至少刷 200 道算法高频题，这些高频题我都写了题解同时也录制了视频，

在这个链接总结了：https://www.algomooc.com/1659.html

2、面试前一周以看题为主，因为刷题也刷不了几题，多看看自己总结或者别人总结的模板，比如回溯算法模板，掌握后，几十道回溯题都不在话下。

一些模板：

3、刷题过程需要注意难度要循序渐进，算法训练是一个系统工程，需要循序渐进，太过于急功近利，反而容易因做不出难题而产生挫败感，带来反效果。

如果你本身有基础，熟练度高，那你刷简单的 LeetCode 应该是几分钟一题，几分钟一题的，花不了你多少时间。

如果你刷简单都花费很长时间，说明熟练度不够，就更应该从简单开始，然后过度到中等，再过度到困难。

并且，目前国内大厂的算法考察，基本不会超过 LeetCode 中等难度，上限难度基本都是 LeetCode 中等题里面的中等难度，所以不要太去纠结难题怪题偏题。

把高频题掌握就行了：https://www.algomooc.com/1659.html

再退一步，如果你觉得 LeetCode 的题目太难，可以先从《剑指 Offer》上的算法题开始学起。

为了帮助大家更好的入门学习算法，经过半年的积累，我给大家卷了《剑指 Offer》系列的三十道题目，结合动画的形式录制了视频，相信能帮助你更好的刷题。

领取地址：

4、按算法分类来选题，比如一个时间段，只刷链表题，刷得差不多的时候，接下来再刷二叉树的题。

这样做有几个很明显的好处。

一、持续地刷同个类型的题目，可以不断地巩固和加深理解，可以总结出自己的思考路径或者解题模板。

比如链表题目，就会去思考虚拟头节点、双指针、快慢指针。

二、可以更全面地接触这个数据结构，算法的各个变种，这会促使你对这个数据结构，算法的理解更加全面和深刻，学习的效率会更高。

我一直认为读书是世界上性价比最高的成长方式，书很便宜但分量很重，是让我们摆脱平庸走向卓越的方式之一。

对于计算机专业的学生而言，读计算机经典书籍不光能让你快速提升知识和能力，更会让你在校招之际如虎添翼。

书籍下载：计算机必看经典书籍（含下载方式）

最后，再给大家送上点干货！

下面这是一个高赞回答合集，建议大家点赞&收藏，Mark住别丢了，大学期间绝对用得上。

1、怎么学好数据结构，看下面这个回答，已经获得了 21000+ 的赞和 50000+的收藏。

2、如何系统地学习算法，看下面这个回答，已经获得了 11000+ 的赞和 26000+的收藏。

3、新手该如何使用 GitHub，看下面这个回答，如果在大学期间就知道使用 GitHub ，那么能力远超同龄人。

4、想成为一名优秀的程序员，那么这些程序员平时都喜欢逛的论坛怎么说你也得收藏一些吧。

5、无论别人怎么说，我都是坚定不移的选择计算机专业。

6、如何系统地学习 C++ ，这个回答能帮你找到路线。

7、想要准备 Java 面试，那么这些面试题必须掌握。

赶紧点赞和收藏吧~

没有人提王垠的40行代码吗？

简单地说，这段代码做了两件事，一件事是CPS，也就是自动尾递归，第二件事则是用Scheme语言写了一个Scheme的解释器。通过他给出的cps函数，我可以用Scheme这个语言的符号系统重新定义所有Scheme的关键字，并执行正确的程序语义。换言之，它可以让这个语言自己解释自己。本质上，他的代码是在模仿当初 John McCarthy 发明 Lisp 语言时给出的代码，但用了Scheme风格重写了一遍。

-----------------------------------

关于这段代码的内容大家可以到这个问题下去看

王垠的「40 行代码」真如他说的那么厉害吗？ - 知乎 https://www.zhihu.com/question/20822815

也欢迎大家到我的公众号不恰柠檬来玩呀~

机器学习十大算法之一：EM算法。

能评得上十大之一，让人听起来觉得挺NB的。什么是NB啊，我们一般说某个人很NB，是因为他能解决一些别人解决不了的问题。那么EM算法能解决什么问题呢？或者说EM算法是因为什么而来到这个世界上，还吸引了那么多世人的目光。

进入正题之前先给大家推荐一本书--《Python机器学习及实践－从零开始通往KAGGLE竞赛之路》，对于想系统学习机器学习的同学，推荐去看一下，电子版下载链接如下：

链接：https://pan.baidu.com/s/1rzaFErwbtQ6CQdTSRnaURg

提取码：5gou

进入正题：

一、最大似然

• 假设我们遇到的是下面这样的问题：
  
• 假设我们需要调查我们学校的男生和女生的身高分布。你怎么做啊？你说那么多人不可能一个一个去问吧，肯定是抽样了。假设你在校园里随便地活捉了100个男生和100个女生。他们共200个人（也就是200个身高的样本数据，为了方便表示，下面，我说“人”的意思就是对应的身高）都在教室里面了。那下一步怎么办啊？你开始喊：“男的左边，女的右边，其他的站中间！”。然后你就先统计抽样得到的100个男生的身高。假设他们的身高是服从高斯分布的。但是这个分布的均值u和方差∂2我们不知道，这两个参数就是我们要估计的。记作θ=[u, ∂]T。
• 用数学的语言来说就是：在学校那么多男生（身高）中，我们独立地按照概率密度p(x|θ)抽取100了个（身高），组成样本集X，我们想通过样本集X来估计出未知参数θ。这里概率密度p(x|θ)我们知道了是高斯分布N(u,∂)的形式，其中的未知参数是θ=[u, ∂]T。抽到的样本集是X={x1,x2,…,xN}，其中xi表示抽到的第i个人的身高，这里N就是100，表示抽到的样本个数。
• 由于每个样本都是独立地从p(x|θ)中抽取的，换句话说这100个男生中的任何一个，都是我随便捉的，从我的角度来看这些男生之间是没有关系的。那么，我从学校那么多男生中为什么就恰好抽到了这100个人呢？抽到这100个人的概率是多少呢？因为这些男生（的身高）是服从同一个高斯分布p(x|θ)的。那么我抽到男生A（的身高）的概率是p(xA|θ)，抽到男生B的概率是p(xB|θ)，那因为他们是独立的，所以很明显，我同时抽到男生A和男生B的概率是p(xA|θ)* p(xB|θ)，同理，我同时抽到这100个男生的概率就是他们各自概率的乘积了。用数学家的口吻说就是从分布是p(x|θ)的总体样本中抽取到这100个样本的概率，也就是样本集X中各个样本的联合概率，用下式表示：

• 这个概率反映了，在概率密度函数的参数是θ时，得到X这组样本的概率。因为这里X是已知的，也就是说我抽取到的这100个人的身高可以测出来，也就是已知的了。而θ是未知了，则上面这个公式只有θ是未知数，所以它是θ的函数。这个函数放映的是在不同的参数θ取值下，取得当前这个样本集的可能性，因此称为参数θ相对于样本集X的似然函数（likehood function）。记为L(θ)。
• 这里出现了一个概念，似然函数。还记得我们的目标吗？我们需要在已经抽到这一组样本X的条件下，估计参数θ的值。怎么估计呢？似然函数有啥用呢？那咱们先来了解下似然的概念。

直接举个例子：

• 某位同学与一位猎人一起外出打猎，一只野兔从前方窜过。只听一声枪响，野兔应声到下，如果要你推测，这一发命中的子弹是谁打的？你就会想，只发一枪便打中，由于猎人命中的概率一般大于这位同学命中的概率，看来这一枪是猎人射中的。
• 这个例子所作的推断就体现了极大似然法的基本思想。
• 再例如：下课了，一群男女同学分别去厕所了。然后，你闲着无聊，想知道课间是男生上厕所的人多还是女生上厕所的人比较多，然后你就跑去蹲在男厕和女厕的门口。蹲了五分钟，突然一个美女走出来，你狂喜，跑过来告诉我，课间女生上厕所的人比较多，你要不相信你可以进去数数。呵呵，我才没那么蠢跑进去数呢，到时还不得上头条。我问你是怎么知道的。你说：“5分钟了，出来的是女生，女生啊，那么女生出来的概率肯定是最大的了，或者说比男生要大，那么女厕所的人肯定比男厕所的人多”。看到了没，你已经运用最大似然估计了。你通过观察到女生先出来，那么什么情况下，女生会先出来呢？肯定是女生出来的概率最大的时候了，那什么时候女生出来的概率最大啊，那肯定是女厕所比男厕所多人的时候了，这个就是你估计到的参数了。
• 从上面这两个例子，你得到了什么结论？
• 回到男生身高那个例子。在学校那么男生中，我一抽就抽到这100个男生（表示身高），而不是其他人，那是不是表示在整个学校中，这100个人（的身高）出现的概率最大啊。那么这个概率怎么表示？哦，就是上面那个似然函数L(θ)。所以，我们就只需要找到一个参数θ，其对应的似然函数L(θ)最大，也就是说抽到这100个男生（的身高）概率最大。这个叫做θ的最大似然估计量，记为：

• 有时，可以看到L(θ)是连乘的，所以为了便于分析，还可以定义对数似然函数，将其变成连加的：

• 好了，现在我们知道了，要求θ，只需要使θ的似然函数L(θ)极大化，然后极大值对应的θ就是我们的估计。这里就回到了求最值的问题了。怎么求一个函数的最值？当然是求导，然后让导数为0，那么解这个方程得到的θ就是了（当然，前提是函数L(θ)连续可微）。那如果θ是包含多个参数的向量那怎么处理啊？当然是求L(θ)对所有参数的偏导数，也就是梯度了，那么n个未知的参数，就有n个方程，方程组的解就是似然函数的极值点了，当然就得到这n个参数了。
• 最大似然估计你可以把它看作是一个反推。多数情况下我们是根据已知条件来推算结果，而最大似然估计是已经知道了结果，然后寻求使该结果出现的可能性最大的条件，以此作为估计值。比如，如果其他条件一定的话，抽烟者发生肺癌的危险时不抽烟者的5倍，那么如果现在我已经知道有个人是肺癌，我想问你这个人抽烟还是不抽烟。你怎么判断？你可能对这个人一无所知，你所知道的只有一件事，那就是抽烟更容易发生肺癌，那么你会猜测这个人不抽烟吗？我相信你更有可能会说，这个人抽烟。为什么？这就是“最大可能”，我只能说他“最有可能”是抽烟的，“他是抽烟的”这一估计值才是“最有可能”得到“肺癌”这样的结果。这就是最大似然估计。
• 好了，极大似然估计就讲到这，总结一下：
  
• 极大似然估计，只是一种概率论在统计学的应用，它是参数估计的方法之一。说的是已知某个随机样本满足某种概率分布，但是其中具体的参数不清楚，参数估计就是通过若干次试验，观察其结果，利用结果推出参数的大概值。最大似然估计是建立在这样的思想上：已知某个参数能使这个样本出现的概率最大，我们当然不会再去选择其他小概率的样本，所以干脆就把这个参数作为估计的真实值。

• 求最大似然函数估计值的一般步骤：

（1）写出似然函数；

（2）对似然函数取对数，并整理；

（3）求导数，令导数为0，得到似然方程；

（4）解似然方程，得到的参数即为所求；

二、EM算法

• 好了，重新回到上面那个身高分布估计的问题。现在，通过抽取得到的那100个男生的身高和已知的其身高服从高斯分布，我们通过最大化其似然函数，就可以得到了对应高斯分布的参数θ=[u, ∂]T了。那么，对于我们学校的女生的身高分布也可以用同样的方法得到了。
• 再回到例子本身，如果没有“男的左边，女的右边，其他的站中间！”这个步骤，或者说我抽到这200个人中，某些男生和某些女生一见钟情，已经好上了，纠缠起来了。咱们也不想那么残忍，硬把他们拉扯开。那现在这200个人已经混到一起了，这时候，你从这200个人（的身高）里面随便给我指一个人（的身高），我都无法确定这个人（的身高）是男生（的身高）还是女生（的身高）。也就是说你不知道抽取的那200个人里面的每一个人到底是从男生的那个身高分布里面抽取的，还是女生的那个身高分布抽取的。用数学的语言就是，抽取得到的每个样本都不知道是从哪个分布抽取的。
• 这个时候，对于每一个样本或者你抽取到的人，就有两个东西需要猜测或者估计的了，一是这个人是男的还是女的？二是男生和女生对应的身高的高斯分布的参数是多少？
• 只有当我们知道了哪些人属于同一个高斯分布的时候，我们才能够对这个分布的参数作出靠谱的预测，例如刚开始的最大似然所说的，但现在两种高斯分布的人混在一块了，我们又不知道哪些人属于第一个高斯分布，哪些属于第二个，所以就没法估计这两个分布的参数。反过来，只有当我们对这两个分布的参数作出了准确的估计的时候，才能知道到底哪些人属于第一个分布，那些人属于第二个分布。
• 这就成了一个先有鸡还是先有蛋的问题了。鸡说，没有我，谁把你生出来的啊。蛋不服，说，没有我，你从哪蹦出来啊。（呵呵，这是一个哲学问题。当然了，后来科学家说先有蛋，因为鸡蛋是鸟蛋进化的）。为了解决这个你依赖我，我依赖你的循环依赖问题，总得有一方要先打破僵局，说，不管了，我先随便整一个值出来，看你怎么变，然后我再根据你的变化调整我的变化，然后如此迭代着不断互相推导，最终就会收敛到一个解。这就是EM算法的基本思想了。
• 不知道大家能否理解其中的思想，我再来啰嗦一下。其实这个思想无处在不啊。
  
• 例如，小时候，老妈给一大袋糖果给你，叫你和你姐姐等分，然后你懒得去点糖果的个数，所以你也就不知道每个人到底该分多少个。咱们一般怎么做呢？先把一袋糖果目测的分为两袋，然后把两袋糖果拿在左右手，看哪个重，如果右手重，那很明显右手这代糖果多了，然后你再在右手这袋糖果中抓一把放到左手这袋，然后再感受下哪个重，然后再从重的那袋抓一小把放进轻的那一袋，继续下去，直到你感觉两袋糖果差不多相等了为止。呵呵，然后为了体现公平，你还让你姐姐先选了。
• EM算法就是这样，假设我们想估计知道A和B两个参数，在开始状态下二者都是未知的，但如果知道了A的信息就可以得到B的信息，反过来知道了B也就得到了A。可以考虑首先赋予A某种初值，以此得到B的估计值，然后从B的当前值出发，重新估计A的取值，这个过程一直持续到收敛为止。
• EM的意思是“Expectation Maximization”，在我们上面这个问题里面，我们是先随便猜一下男生（身高）的正态分布的参数：如均值和方差是多少。例如男生的均值是1米7，方差是0.1米（当然了，刚开始肯定没那么准），然后计算出每个人更可能属于第一个还是第二个正态分布中的（例如，这个人的身高是1米8，那很明显，他最大可能属于男生的那个分布），这个是属于Expectation一步。有了每个人的归属，或者说我们已经大概地按上面的方法将这200个人分为男生和女生两部分，我们就可以根据之前说的最大似然那样，通过这些被大概分为男生的n个人来重新估计第一个分布的参数，女生的那个分布同样方法重新估计。这个是Maximization。然后，当我们更新了这两个分布的时候，每一个属于这两个分布的概率又变了，那么我们就再需要调整E步……如此往复，直到参数基本不再发生变化为止。
• 这里把每个人（样本）的完整描述看做是三元组yi={xi,zi1,zi2}，其中，xi是第i个样本的观测值，也就是对应的这个人的身高，是可以观测到的值。zi1和zi2表示男生和女生这两个高斯分布中哪个被用来产生值xi，就是说这两个值标记这个人到底是男生还是女生（的身高分布产生的）。这两个值我们是不知道的，是隐含变量。确切的说，zij在xi由第j个高斯分布产生时值为1，否则为0。例如一个样本的观测值为1.8，然后他来自男生的那个高斯分布，那么我们可以将这个样本表示为{1.8, 1, 0}。如果zi1和zi2的值已知，也就是说每个人我已经标记为男生或者女生了，那么我们就可以利用上面说的最大似然算法来估计他们各自高斯分布的参数。但是它们未知，因此我们只能用EM算法。
• 咱们现在不是因为那个恶心的隐含变量（抽取得到的每个样本都不知道是从哪个分布抽取的）使得本来简单的可以求解的问题变复杂了，求解不了吗。那怎么办呢？人类解决问题的思路都是想能否把复杂的问题简单化。好，那么现在把这个复杂的问题逆回来，我假设已经知道这个隐含变量了，哎，那么求解那个分布的参数是不是很容易了，直接按上面说的最大似然估计就好了。那你就问我了，这个隐含变量是未知的，你怎么就来一个假设说已知呢？你这种假设是没有根据的。呵呵，我知道，所以我们可以先给这个给分布弄一个初始值，然后求这个隐含变量的期望，当成是这个隐含变量的已知值，那么现在就可以用最大似然求解那个分布的参数了吧，那假设这个参数比之前的那个随机的参数要好，它更能表达真实的分布，那么我们再通过这个参数确定的分布去求这个隐含变量的期望，然后再最大化，得到另一个更优的参数，……迭代，就能得到一个皆大欢喜的结果了。
• 这时候你就不服了，说你老迭代迭代的，你咋知道新的参数的估计就比原来的好啊？为什么这种方法行得通呢？有没有失效的时候呢？什么时候失效呢？用到这个方法需要注意什么问题呢？呵呵，一下子抛出那么多问题，搞得我适应不过来了，不过这证明了你有很好的搞研究的潜质啊。呵呵，其实这些问题就是数学家需要解决的问题。在数学上是可以稳当的证明的或者得出结论的。那咱们用数学来把上面的问题重新描述下。（在这里可以知道，不管多么复杂或者简单的物理世界的思想，都需要通过数学工具进行建模抽象才得以使用并发挥其强大的作用，而且，这里面蕴含的数学往往能带给你更多想象不到的东西，这就是数学的精妙所在啊）

三、EM算法推导

• 假设我们有一个样本集{x(1),…,x(m)}，包含m个独立的样本。但每个样本i对应的类别z(i)是未知的（相当于聚类），也即隐含变量。故我们需要估计概率模型p(x,z)的参数θ，但是由于里面包含隐含变量z，所以很难用最大似然求解，但如果z知道了，那我们就很容易求解了。
• 对于参数估计，我们本质上还是想获得一个使似然函数最大化的那个参数θ，现在与最大似然不同的只是似然函数式中多了一个未知的变量z，见下式（1）。也就是说我们的目标是找到适合的θ和z让L(θ)最大。那我们也许会想，你就是多了一个未知的变量而已啊，我也可以分别对未知的θ和z分别求偏导，再令其等于0，求解出来不也一样吗？

• 本质上我们是需要最大化（1）式（对（1）式，我们回忆下联合概率密度下某个变量的边缘概率密度函数的求解，注意这里z也是随机变量。对每一个样本i的所有可能类别z求等式右边的联合概率密度函数和，也就得到等式左边为随机变量x的边缘概率密度），也就是似然函数，但是可以看到里面有“和的对数”，求导后形式会非常复杂（自己可以想象下log(f1(x)+ f2(x)+ f3(x)+…)复合函数的求导），所以很难求解得到未知参数z和θ。那OK，我们可否对（1）式做一些改变呢？我们看（2）式，（2）式只是分子分母同乘以一个相等的函数，还是有“和的对数”啊，还是求解不了，那为什么要这么做呢？咱们先不管，看（3）式，发现（3）式变成了“对数的和”，那这样求导就容易了。我们注意点，还发现等号变成了不等号，为什么能这么变呢？这就是Jensen不等式的大显神威的地方。

Jensen不等式：

• 设f是定义域为实数的函数，如果对于所有的实数x。如果对于所有的实数x，f(x)的二次导数大于等于0，那么f是凸函数。当x是向量时，如果其hessian矩阵H是半正定的，那么f是凸函数。如果只大于0，不等于0，那么称f是严格凸函数。
• Jensen不等式表述如下：
• 如果f是凸函数，X是随机变量，那么：E[f(X)]>=f(E[X])
• 特别地，如果f是严格凸函数，当且仅当X是常量时，上式取等号。
• 如果用图表示会很清晰：
  

• 图中，实线f是凸函数，X是随机变量，有0.5的概率是a，有0.5的概率是b。（就像掷硬币一样）。X的期望值就是a和b的中值了，图中可以看到E[f(X)]>=f(E[X])成立。
• 当f是（严格）凹函数当且仅当-f是（严格）凸函数。
  
• Jensen不等式应用于凹函数时，不等号方向反向。

• 回到公式（2），因为f(x)=log x为凹函数（其二次导数为-1/x2<0）。
• （2）式中

的期望，（考虑到E(X)=∑x*p(x)，f(X)是X的函数，则E(f(X))=∑f(x)*p(x)），又，所以就可以得到公式（3）的不等式了（若不明白，请拿起笔，呵呵）：

• OK，到这里，现在式（3）就容易地求导了，但是式（2）和式（3）是不等号啊，式（2）的最大值不是式（3）的最大值啊，而我们想得到式（2）的最大值，那怎么办呢？
• 现在我们就需要一点想象力了，上面的式（2）和式（3）不等式可以写成：似然函数L(θ)>=J(z,Q)，那么我们可以通过不断的最大化这个下界J，来使得L(θ)不断提高，最终达到它的最大值。

• 见上图，我们固定θ，调整Q(z)使下界J(z,Q)上升至与L(θ)在此点θ处相等（绿色曲线到蓝色曲线），然后固定Q(z)，调整θ使下界J(z,Q)达到最大值（θt到θt+1），然后再固定θ，调整Q(z)……直到收敛到似然函数L(θ)的最大值处的θ*。这里有两个问题：什么时候下界J(z,Q)与L(θ)在此点θ处相等？为什么一定会收敛？
• 首先第一个问题，在Jensen不等式中说到，当自变量X是常数的时候，等式成立。而在这里，即：

• 再推导下，由于

（因为Q是随机变量z(i)的概率密度函数），则可以得到：分子的和等于c（分子分母都对所有z(i)求和：多个等式分子分母相加不变，这个认为每个样例的两个概率比值都是c），则：

• 至此，我们推出了在固定参数θ后，使下界拉升的Q(z)的计算公式就是后验概率，解决了Q(z)如何选择的问题。这一步就是E步，建立L(θ)的下界。接下来的M步，就是在给定Q(z)后，调整θ，去极大化L(θ)的下界J（在固定Q(z)后，下界还可以调整的更大）。那么一般的EM算法的步骤如下：

EM算法（Expectation-maximization）：

• 期望最大算法是一种从不完全数据或有数据丢失的数据集（存在隐含变量）中求解概率模型参数的最大似然估计方法。

EM的算法流程：

• 初始化分布参数θ；

重复以下步骤直到收敛：

• E步骤：根据参数初始值或上一次迭代的模型参数来计算出隐性变量的后验概率，其实就是隐性变量的期望。作为隐藏变量的现估计值：

• M步骤：将似然函数最大化以获得新的参数值：

• 这个不断的迭代，就可以得到使似然函数L(θ)最大化的参数θ了。那就得回答刚才的第二个问题了，它会收敛吗？
• 感性的说，因为下界不断提高，所以极大似然估计单调增加，那么最终我们会到达最大似然估计的最大值。理性分析的话，就会得到下面的东西：

• 具体如何证明的，看推导过程参考：Andrew Ng《The EM algorithm》:
  http://www.cnblogs.com/jerrylead/archive/2011/04/06/2006936.html

四、EM算法另一种理解

• 坐标上升法（Coordinate ascent）：

• 图中的直线式迭代优化的路径，可以看到每一步都会向最优值前进一步，而且前进路线是平行于坐标轴的，因为每一步只优化一个变量。
• 这犹如在x-y坐标系中找一个曲线的极值，然而曲线函数不能直接求导，因此什么梯度下降方法就不适用了。但固定一个变量后，另外一个可以通过求导得到，因此可以使用坐标上升法，一次固定一个变量，对另外的求极值，最后逐步逼近极值。对应到EM上，E步：固定θ，优化Q；M步：固定Q，优化θ；交替将极值推向最大。

五、EM的应用

• EM算法有很多的应用，最广泛的就是GMM混合高斯模型、聚类、HMM等等。具体可以参考JerryLead的cnblog中的Machine Learning专栏：
• （EM算法）The EM Algorithm
• 混合高斯模型（Mixtures of Gaussians）和EM算法
• K-means聚类算法

更多计算机相关学习资料可以去我的个人网站：https://tanqingbo.cn/CSBook001/

推荐阅读：