EdmondFrank's 时光足迹

この先は暗い夜道だけかもしれない それでも信じて進むんだ。星がその道を少しでも照らしてくれるのを。
或许前路永夜,即便如此我也要前进,因为星光即使微弱也会我为照亮前途。
——《四月は君の嘘》

文本建模算法入门-1



文本建模算法入门(一)

0x1 文本建模

文本并不像其他数值型的数据一样可以比较轻易的通过运算、函数、方程、矩阵等来表达他们之间的相互关系。

在处理一篇文本的时候,假设每一个文本储存为一篇文档,那么从人的视角来看,这可以说是一段有序的词序列
统计学家将这些序列的生成,生动地描绘成了一个“上帝的游戏”,即人类产生的所有语料的文本都可以看作是:一个伟大的上帝在天堂中掷骰子形成的。我们所看到的文本其实就是这个游戏掷了若干次后产生的序列。

那么,在这个游戏中,我们最需要关注的两个核心问题就出现了:

  1. 上帝有怎样的骰子?
  2. 上帝是怎么掷的骰子?

对应这两个问题,各大学家持着不同的观点,于是便有了以下三种模型:

  1. Unigram Model
  2. Topic Model(PLSA)
  3. LDA

0x11 Unigram Model

Unigram Model是非常简单直接的,它假设

  1. 上帝只有一个V面的骰子,每一个面对应一个词,同时每个面的概率不一样。(这里可以通过“老千骰子”来理解下,有些面因为被做过了“手脚”,所以抛到的几率就大了,那如果没个面都这样,那么每个面抛到的几率也就不一样了)
  2. 每抛一次骰子,抛出的面就对应有一个词,那么抛n次骰子后,按抛掷顺序产生的序列就生成了一篇n个词的文档

enter image description here

那现在我们把上帝这个唯一的骰子各个面的概率记为,然后我们把掷这个V面骰子的实验记作

那么对应一篇文档d,该文档生成的概率就是

而文档和文档之间我们认为是独立的,所以如果语料中有多篇文档,则该语料生成的概率就是

在Unigram Model中,我们又假设了文档之间是独立可交换的。即,词与词之间的顺序对文档表示并不造成影响,一篇文档相当于一个袋子,里面装着一些词。这样的模型也称为词袋模型(Bag-of-words)。

那么,如果语料中的总词频是N,在N个词中,如果我们关注每一个词的发生次数,则正好是一个多项分布

此时,语料的概率是

0x110 贝叶斯Unigram Model假设

在贝叶斯统计学派看来,上帝拥有唯一一个固定的骰子是不合理的。于是他们觉得以上模型的骰子不应该唯一固定,而应该也是一个随机变量。

那这样我们就可以将整个掷骰子的游戏过程更新成以下形式:

  1. 上帝有一个装着无穷多骰子的罐子,里面有各种骰子,每个骰子有V个面。每一个面对应一个词
  2. 上帝从罐子抽出一个骰子,然后用这个骰子不断的抛,每抛一次骰子,抛出的面就对应有一个词,那么抛n次骰子后,按抛掷顺序产生的序列就生成了一篇n个词的文档
    enter image description here

在以上的假设之下,由于我们事先并不知道上帝用了哪个骰子,所以每个骰子都是有可能被使用的,只是使用的概率由先验分布决定,对每一个具体的骰子,由该骰子产生数据的概率是,所以最终数据产生的概率就是对每个骰子上产生的数据概率进行积分累加求和

在贝叶斯分析的框架之下,此处的先验分布概率其实就是一个多项分布的概率,其中一个比较好的选择即多项分布对应的共轭分布:Dirichlet分布

0x12 PLSA Topic Model

再来回顾Unigram Model我们发现:这个模型的假设过于简单,和人类真实的书写有着较大的差距

从人类视角看,我们在日常构思文章中,我们往往要先确定自己文章的主旨,包含了哪些主题,然后再围绕着这些主题展开阐述。

例如,一篇关于现代教育的文章,它可能就包含了这些主题:教育方法、多媒体技术、互联网等。篇幅上可能以教育方法为主,而其他为辅。然后,在不同的主题里面就包含了许多主题领域内常见的关键词。例如,谈到互联网时,我们会提及Web、Tcp等。

这样一种直观的想法就在PLSA模型中进行了明确的体现,如果我们再利用这个想法更新“掷骰子”的游戏就有以下情形:

  1. 上帝有两种类型骰子,一类是doc-topic骰子,骰子共k个面代表了k个主题;一类是topic-word骰子,骰子V个面,每面对应主题内的一个词
  2. 生成文章的时候,首先要先创造一个特定的doc-topic骰子,使得骰子内的主题围绕文章要阐述的主题
  3. 投掷doc-topic骰子,得到一个主题z
  4. 根据得到的主题z,找到对应的topic-word骰子,投掷它得到一个词
  5. 不断重复3,4步,直至文章生成完成

enter image description here

在上面的游戏规则中,文档与文档之间顺序无关,同一个文档内的词的顺序也是无关的。所以还是一个bag-of-words模型。

那么,在第m篇文档Dm中每个词的生成概率为:

:对应游戏中K个topic-word骰子中第z个骰子对应的词列表

:文档对应的第z个主题,即对应的doc-topic

所以整篇文档的生成概率就为:

0x13 LDA(Latent Dirichlet Allocation)

就像Unigram Model 加入贝叶斯框架那样,doc-topic和topic-word都是模型的参数,即随机变量。于是类似对Unigram Model的改造对以上两个骰子模型加入先验分布。然后由于都对应着多项分布,因此先验分布依旧选择Drichlet分布。于是得到的这个新模型就是LDA(Latent Dirichlet Allocation)模型
enter image description here

在LDA模型中,上帝的游戏规则就又被更新成如下情形:

  1. 上帝有两个罐子,第一个装着都哦doc-topic骰子,第二个装着topic-word骰子
  2. 上帝随机的从第二个罐子中独立抽出K个topic-doc骰子,编号1-K
  3. 每次生成新文档时,从第一个罐子随机抽一个doc-topic骰子
  4. 投掷得到的doc-topic骰子。得到一个主题编号z
  5. 在K个主题骰子里面选择编号为z的骰子,投掷骰子,得到一个词
  6. 重复4,5步,直至文档生成完成

0x2 后记

至此,入门篇(一)的内容就写到这了。由于是第一篇,文章中避过了许多较为复杂的数学证明和计算,尤其是关于LDA模型的。这样是为了,先建立对文本建模思路和过程的直观认识,而不是一上来就深究细节。

加上笔者目前也是在学习阶段,后面的细节再慢慢地一一补充,大家共同进步 !!
\(^_^)/

SLAM的复兴-无人驾驶的未来?

SLAM的前世今生

SLAM,全称(simultaneous localization and mapping),中文译做,即时定位与地图构建。是指运动物体根据传感器的信息,一边计算自身位置,一边构建环境地图的过程,解决机器人等在未知环境下运动时的定位与地图构建问题。目前,SLAM的主要应用于机器人、无人机、无人驾驶、AR、VR等领域。其用途包括传感器自身的定位,以及后续的路径规划、运动性能、场景理解。

作为一种基础技术,最早可以追溯到当初的军事潜艇定位技术。而今年来随着扫地机器人的盛行,再次令他名声大噪。加上最近的无人驾驶与机器人的寻迹的崛起以及基于三维视觉的VSLAM的出现又让它越来越显主流。

目前用在SLAM上的Sensor主要分两大类,分别是激光雷达和摄像头。

SLAM到底做了些什么?

既然SLAM技术代表了即时定位与地图构建,那么可以说这就是一个机器人或设备尝试去根据他周围的环境创建地图的过程,并且还可以在该地图中实时定位自己。

这不是一件容易的事,这项技术现在还处在技术研究和设计的前沿。而为了成功实施SLAM有一大障碍是不可i避免的,那就是地图定位问题,两个问题同时引入就变成了我们常说的典型的鸡与蛋问题。为了能够成功地顺利的根据环境创建地图,我们的设备就必须先知道它的方向和位置 ;然而,这些定位信息设备也只能从预先存在的环境地图中获得。

面对这一障碍,SLAM技术通常通过使用GPS数据构建预先存在的环境地图来克服这一复杂的鸡与蛋问题。随着机器人或设备在环境中移动,生成的地图会被迭代地改进。而这项技术的真正挑战是准确性。随着机器人或设备在空间中移动与寻迹,测量评估必须不断进行,并且该技术必须考虑设备移动和测量方法不准确而引起的“噪音”。

enter image description here

这个gif是宾大的教授kumar做的特别有名的一个demo,是在无人机上利用二维激光雷达做的SLAM。

VSLAM又是什么?

VSLAM(基于视觉的定位与建图):随着计算机视觉的迅速发展,视觉SLAM因为信息量大,适用范围广等优点受到广泛关注。

(1)基于深度摄像机的Vslam,跟激光SLAM类似,通过收集到的点云数据,能直接计算障碍物距离;

(2)基于单目、鱼眼相机的VSLAM方案,利用多帧图像来估计自身的位姿变化,再通过累计位姿变化来计算距离物体的距离,并进行定位与地图构建;

enter image description here

以上为百度 VSLAM demo展示。

VSLAM能够为自动驾驶带来些什么?

结合自动驾驶的场景,可以推出VSLAM的应用点主要是:

  • gps缺失场景下的长时间定位,如室内,楼房中。
  • 补偿行驶过程中gps信号不稳定造成的定位跳跃,如山洞,高楼群,野外山区等。

VSLAM的精度及鲁棒性越高,适用的场景越宽广。如果VSLAM能在任何场景无限长定位保持高精度,那其他定位技术就可以下岗了,虽然按目前看来这个可能性很小,因此需要尝试结合IMU,编码器等设备融合。至于最终的技术形态,目前还没有定论。

在无人驾驶汽车上,目前比较显著的瓶颈还在计算方面,数据收发,障碍物感知,融合定位,路径规划,每个模块都需要占据相当的计算资源。而无论视觉还是激光SLAM本身就是对计算力消耗极大的一个模块,如果加上高频率的IMU进行融合优化,则计算力更加捉襟见肘。因此是否值得为视觉SLAM分配原本就那么有限的计算资源,如何识别场景对SLAM模块进行激活都是需要仔细衡量的。

我的算法天梯之路之-穷举搜索

穷举搜索的例子-Google方程式

 问题描述

有一个由字符组成的等式,WWWDOT-GOOGLE = DOTCOM,每个字符代表一个0-9之间的数字,WWWDOT、GOOGLE和DOTCOM都是合法的数字,不能以零开头。请找出一组字符与数字的对应关系,使得他们可以互相转换,并且替换之后能够满足等式。

问题分析

从排列组合的角度来看,这道题是一道典型的排列组合问题,题目中一共出现了9个字母。

如果不考虑0开头的情况下,这样的组合应该有10x9x8x7x6x5x4x3x2=3628800种组合。

在这样的情况之下,计算机的穷举处理应该是毫无压力的。

问题求解

首先为了能够表示这样一种可变的字符元素列表,我们需要自定义一种数据结构。首先我们知道这个自定义结构中应该包含有三个属性,分别是 字符本身,字符代表的数字以及是否为数字的最高位,因为最高位是不能为0的,所以在这里我们要区别对待。

1
2
3
4
5
6
typedef struct tagCharItem
{
  char c;
  int value;
  bool leading;
}CHAR_ITEM;

接着我们初始化这个列表。=

1
2
3
4
5
CHAR_ITEM charItem[]={
{'W',-1,true},{'D',-1,true},{'O',-1,true},
{'T',-1,false},{'G',-1,true},{'L',-1,false},
{'E',-1,false},{'C',-1,false},{'M',-1,false}
};

因为这是一个组合问题,那么两个字母就不能被指定为相同的数字,这样我们需要定义额外的占用标识。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
typedef struct tagCharValue
{
  bool isused;
  int value;
}CHAR_VALUE;

//穷举算法实现如下
void SearchingResult(CHAR_ITEM ci[],CHAR_VALUE cv[],
int index,CharListReadyFuncPtr callback)
{
if(index == max_char_count)
{
  callback(ci);
  return;
}
for(int i = 0;i<max_number_count;++i)
{
  if(IsValueValid(ci[index],cv[i]))
  {
      cv[i].isused= true;
      ci[index].value = cv[i].value;
      SearchingResult(ci,cv,index+1,callback);
      cv[i].isused = false;
  }
}
}

根据题目要求,W,G,D三者都不能为0,为了加快穷举速度可以对他们为0的情况进行剪枝。
IsValueValid是评估函数,在剪枝操作之后,callback的被调用次数可以减少约30%。

1
2
3
4
5
6
7
8
9
10
11
12
13
//callback代码实现
void OnCharListReady(CHAR_ITEM ci[])
{
  char * minuend = "WWWDOT";
  char * subtrahend = "GOOGLE";
  char * diff = "DOTCOM";
  int m = MakeIntegerValue(ci,minuend);
  int s = MakeIntegerValue(ci,subtrahend);
  int d = MakeIntegerValue(ci,diff);
  if((m-s)==d){
      std::cout << m << "-" <<s<< "=" << d << std::endl;
  }
}

问题答案
777589 - 188103 = 589486
777589 - 188106 = 589483

隐马尔可夫模型(HMM)的原理



隐马尔可夫模型(HMM)的原理

简介

隐马尔可夫模型是一个基于马尔可夫链的统计模型。

马尔可夫链安德烈·马尔可夫(Andrey Markov,1856-1922)得名,是数学中具有马尔可夫性质的离散时间随机过程。该过程中,在给定当前知识或信息的情况下,过去(即当期以前的历史状态)对于预测将来(即当期以后的未来状态)是无关的。

虽然从维基百科上摘取下来的概念和定义看着十分的晦涩难懂,但是模型背后的思想是非常简单的,即:首先假设你的系统可以建模为马尔可夫链,然后,系统所发出的信号(输出的可见结果)仅取决于系统当前的状态。那么,隐马尔可夫模型代表的一种情景就是:系统的状态对你而言是不可见的,你仅仅只能观测到系统所发出的信号或者说是系统所输出的结果。

举个通俗易懂的栗子。

假设你有一个住在国外的朋友,他通常会根据天气来安排他的日常活动。你不知道他的国家那边的天气如何(系统的状态),但你确可以在跟他的聊天中知道他今天进行了什么活动(系统的输出),然后这就是一个简单的隐马尔可夫模型。

我们将整个模型
markov.jpg

现在有三个比较关键的问题有待我们解决:
1. 知道整个模型后,你朋友告诉你他这三天的活动是:散步(Walk),购物(Shop),清洁屋子(Clean),那么根据模型,计算产生这些行为序列的概率是多少?
2. 知道整个模型后,朋友让你根据他的活动猜一猜他那边这三天天气怎么样
3. 朋友告诉你三天里他做了些什么,然后让你找出他活动规律的模型。

为了解决以上的问题我们首先要了解一些有关HMM的基本元素先:

初始概率分布:初始概率分布即事件初始时发生的概率,我们这里隐藏的状态是天气,然后我们在图中可以看出的初始概率有:天气为晴天的概率为0.8 ; 天气为雨天的概率为0.2

转移概率矩阵P:转移功率矩阵就可以通过一副图来描述了。

晴天 雨天
晴天 0.7
雨天 0.6

其中,列参数表示第一天的状态,行参数表示第二天状态。表格中的第一行的含义就是已知第一天为晴天,那么第二天为晴天的概率是0.7,而为雨天的概率只是0.3。

观测量的概率分布B:在这个问题中观测量就是朋友的活动,其概率分布就分别表示为朋友在晴天和雨天的情况下进行散步,购物,打扫屋子各项活动的可能性(概率)。

现在,我们再次回到上面提到的三个问题。其实,每个问题的解决解决在历史上早有前人给出了算法。

问题1 -> Forward Algorithm,向前算法 或 Backward Algorithm,向后算法。

问题2 -> Viterbi Algorithm,维特比算法。

问题3 -> Baum-Welch Algorithm,鲍姆-维尔奇算法。

算法思路

假设前提:
已知,雨天,朋友选择去散步,购物,收拾的概率分别是0.1,0.4,0.5, 而如果是晴天,选择去散步,购物,收拾的概率分别是0.6,0.3,0.1。
三天活动序列:散步(Walk),购物(Shop),清洁屋子(Clean)

Forward Algorithm

然后我们先计算 t = 1 时,发生 “散步” 行为的概率,如果下雨,则;如果为晴天,则

t = 2 时,发生 “购物” 的概率,如果 t = 2 下雨,则

如果为晴天,则

t = 3 时的算法也可以依此类推,

所以,最终:

从上面的例子可以看出,向前算法计算了每个时间点时,每个状态的发生观测序列的概率,看似复杂,但在T变大时,复杂度也会随之降低。

Backward Algorihm

既然,向前算法是在时间 t = 1 的时候,一步一步向前计算。那么反过来,向后算法就是从最后一个状态往前推。
假设最初时

那么:


其中第一项则是转移概率,第二天下雨转到第三天下雨的概率为0.4;第二项则是观测概率,第三天下雨的状况下,在家收拾的概率为0.5;第三项就是我们定义的向后变量(backward variable)。

同理也可推得其他数据,并且最终答案与向前算法的求解相同。

Viterbi Algorithm(维比特算法)

利用动态规划求解概率最大的路径(最优路径)。利用动态规划,可以解决任何一个图中的最短路径问题。而维特比算法是针对一个特殊的图——篱笆网络的有向图(Lattice)的最短路径提出的。

我们假设用符号来表示系统的第 i 种状态的第j个可能的值。如果把每个状态按照不同的值展开,就可以得到以下的篱笆网络(Lattice):

lattice.png

那么从第一个状态到最后一个状态的任何一条路径(path)都可能产生我们观察到的输出序列Y。当然,这些路径的可能性不一样,而我们要做的就是找到最可能的这条路径。对于每一条给定的路径我们都可以用公式:

计算出它的可能性,但是随着组合增多,它使得序列状态数的增长呈指数爆炸式。

为了解决这个问题,需要一个最好能和状态数成正比的算法。也就是我们要讲的维特比算法。

维特比算法的基础可以概括成三点:

  1. 如果概率最大的路径P(或者说最优路径)经过某个点,如果上图中的,那么这条路径上从起始点S 到 的这一段子路径Q,一定是S到的最短路径。否则用S到的最短路径R来代替Q,便构成了一条比P更短的路径。

  2. 从S到E的路径必定经过第i个时刻的某个状态,假定第i个时刻有k个状态,那么如果记录了从S到第i个状态的所以k个节点的最短路径,最终最短路径必经过其中的一条。这样,在任何时刻,只要考虑非常有限条最短路径即可。

  3. 结合上述两点,假定当我们从状态 i 进入状态 i + 1 时,从 S 到状态 i 上各个节点的最短路径已经找到,并且记录在这些节点上,那么在计算从起点 S 到第 i + 1 状态的某个节点的最短路径时,只要考虑从S到前一个状态 i 所有的 k 个节点的最短路径,以及从这 k 个节点到,j 的距离即可。

实现过程

  1. 从点S出发,对于第一个状态的各个节点,不妨假定有个。计算出 S 到它们的距离,其中代表任意状态 1 的节点。因为只有一步,所以这些距离都是S到它们各自的最短距离。
  2. 对于第二个状态的所有节点,要计算出从S到它们的最短距离。我们知道,对于特定的节点,从S到它的路径可以经过状态1的中任何一个节点,当然,对应的路径长度就是。由于 j 有种可能性,我们要一一计算,然后找到最小值。即:
  3. 这样对于第二状态的每个节点,需要次乘法计算。假定这个状态有个节点,把S这些节点的距离都算一遍,就有次计算。

接下来,类似地按照上述方法从第二个状态走到第三个状态,一直走到最后一个状态,就得到了整个网格从头到尾的最短路径。每一步计算的复杂度都和相邻两个状态各自的节点数目的乘积成正比,即。如果假定在这个隐含马尔可夫链中的节点最多的状态有D个节点,也就是说整个网络的宽度为D,那么任何一布的复杂度不超过,由于网络长度是N,所以整个维特比算法的复杂度是 ——本段推理节选自<<数学之美>>第二版

贝叶斯景象图



贝叶斯景象图

理论说明

均匀分布

对于一个含有N个未知元素的贝叶斯推断问题,我们隐式地为其先验分布创建了一个N维空间。先验分布上某一点的概率,都投射到某个高维的面或曲线上,其形状由先验分布决定。比如,假定有两个未知元素 ,其先验分布都是(0,5)上的均匀分布,那么先验分布存在于一个边长为5的正方形空间,而其概率面就是正方形上方的一个平面(由于假定了均匀分布,因此每一点概率相同)。

代码绘图演示

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
%matplotlib inline
import scipy.stats as stats
from IPython.core.pylabtools import figsize
import numpy as np
figsize(12.5, 4)

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

jet = plt.cm.jet
fig = plt.figure()
x = y = np.linspace(0, 5, 100)
X, Y = np.meshgrid(x, y)

plt.subplot(121)
uni_x = stats.uniform.pdf(x, loc=0, scale=5)
uni_y = stats.uniform.pdf(y, loc=0, scale=5)
M = np.dot(uni_x[:, None], uni_y[None, :])
im = plt.imshow(M, interpolation='none', origin='lower',
                cmap=jet, vmax=1, vmin=-.15, extent=(0, 5, 0, 5))

plt.xlim(0, 5)
plt.ylim(0, 5)
plt.title("Landscape formed by Uniform priors.")

ax = fig.add_subplot(122, projection='3d')
ax.plot_surface(X, Y, M, cmap=plt.cm.jet, vmax=1, vmin=-.15)
ax.view_init(azim=390)
plt.title("Uniform prior landscape; alternate view");

指数分布

再者,如果 的先验分布为Exp(3)和Exp(10),那么对应的空间便是二维平面上,各维都取正值确定的范围,而对应的概率面的形状就是一个从(0,0)点向正值方向流淌的瀑布。

以下的示例图就描绘了这样的情形,其中颜色越是趋向于暗红的位置,其先验概率就越高。反过来,颜色越是趋向于深蓝的位置,其先验概率就越低。

代码绘图演示

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
figsize(12.5, 5)
fig = plt.figure()
plt.subplot(121)

exp_x = stats.expon.pdf(x, scale=3)
exp_y = stats.expon.pdf(x, scale=10)
M = np.dot(exp_x[:, None], exp_y[None, :])
CS = plt.contour(X, Y, M)
im = plt.imshow(M, interpolation='none', origin='lower',
                cmap=jet, extent=(0, 5, 0, 5))
#plt.xlabel("prior on $p_1$")
#plt.ylabel("prior on $p_2$")
plt.title("$Exp(3), Exp(10)$ prior landscape")

ax = fig.add_subplot(122, projection='3d')
ax.plot_surface(X, Y, M, cmap=jet)
ax.view_init(azim=390)
plt.title("$Exp(3), Exp(10)$ prior landscape; \nalternate view");

这些二维空间的例子很简单,我们的大脑可以轻易想象得到。但实际中,先验分布所在的空间和其概率面往往具有更高的维度。

观测值对先验分布的影响 

在实际中,观测样本对空间不会有影响,但它会改变概率面的形状,将其在某些局部区域拉伸或挤压,以表明参数的真实性在哪里。更多的数据意味着对概率面更多的拉伸和挤压,使得最初的概率面形状变得十分奇怪。反之数据越少,那么最初的形状就保留得越好。不管如何,最后得到的概率面描述了后验分布的形状。

假如我们现在想对两个参数为的泊松分布进行估计。那么我们将要分别比较用均匀分布和指数分布来对的先验分布进行假设的不同效果。

代码绘图演示

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
# create the observed data

# sample size of data we observe, trying varying this (keep it less than 100 ;)
N = 1

# the true parameters, but of course we do not see these values...
lambda_1_true = 1
lambda_2_true = 3

#...we see the data generated, dependent on the above two values.
data = np.concatenate([
    stats.poisson.rvs(lambda_1_true, size=(N, 1)),
    stats.poisson.rvs(lambda_2_true, size=(N, 1))
], axis=1)
print("observed (2-dimensional,sample size = %d):" % N, data)

# plotting details.
x = y = np.linspace(.01, 5, 100)
likelihood_x = np.array([stats.poisson.pmf(data[:, 0], _x)
                        for _x in x]).prod(axis=1)
likelihood_y = np.array([stats.poisson.pmf(data[:, 1], _y)
                        for _y in y]).prod(axis=1)
L = np.dot(likelihood_x[:, None], likelihood_y[None, :])
observed (2-dimensional,sample size = 1): [[0 2]]
In [4]:
figsize(12.5, 12)
# matplotlib heavy lifting below, beware!
plt.subplot(221)
uni_x = stats.uniform.pdf(x, loc=0, scale=5)
uni_y = stats.uniform.pdf(x, loc=0, scale=5)
M = np.dot(uni_x[:, None], uni_y[None, :])
im = plt.imshow(M, interpolation='none', origin='lower',
                cmap=jet, vmax=1, vmin=-.15, extent=(0, 5, 0, 5))
plt.scatter(lambda_2_true, lambda_1_true, c="k", s=50, edgecolor="none")
plt.xlim(0, 5)
plt.ylim(0, 5)
plt.title("Landscape formed by Uniform priors on $p_1, p_2$.")

plt.subplot(223)
plt.contour(x, y, M * L)
im = plt.imshow(M * L, interpolation='none', origin='lower',
                cmap=jet, extent=(0, 5, 0, 5))
plt.title("Landscape warped by %d data observation;\n Uniform priors on $p_1, p_2$." % N)
plt.scatter(lambda_2_true, lambda_1_true, c="k", s=50, edgecolor="none")
plt.xlim(0, 5)
plt.ylim(0, 5)

plt.subplot(222)
exp_x = stats.expon.pdf(x, loc=0, scale=3)
exp_y = stats.expon.pdf(x, loc=0, scale=10)
M = np.dot(exp_x[:, None], exp_y[None, :])

plt.contour(x, y, M)
im = plt.imshow(M, interpolation='none', origin='lower',
                cmap=jet, extent=(0, 5, 0, 5))
plt.scatter(lambda_2_true, lambda_1_true, c="k", s=50, edgecolor="none")
plt.xlim(0, 5)
plt.ylim(0, 5)
plt.title("Landscape formed by Exponential priors on $p_1, p_2$.")

plt.subplot(224)
# This is the likelihood times prior, that results in the posterior.
plt.contour(x, y, M * L)
im = plt.imshow(M * L, interpolation='none', origin='lower',
                cmap=jet, extent=(0, 5, 0, 5))

plt.scatter(lambda_2_true, lambda_1_true, c="k", s=50, edgecolor="none")
plt.title("Landscape warped by %d data observation;\n Exponential priors on \
$p_1, p_2$." % N)
plt.xlim(0, 5)
plt.ylim(0, 5);

四张图里的黑点代表了参数的真实取值,左下图为均匀先验得到的后验分布图。虽然观测值相同,但是两种假设下的后验分布形状是不一样的。其主要原因是因为观测点的位置在两者的假设的前提先验概率是不一样的。这样,我们可以知道,即便只有一个观测值,形成的山峰也试图要包括参数值的真实位置。当然,在真正的推断中,仅用一个观测值显然也是十分不科学的,这里仅仅为了方便阐述而已。

本文参考自《Probabilistic-Programming-and-Bayesian-Methods-for-Hackers》

Python Keras + LSTM 进行单变量时间序列预测

Python Keras + LSTM 进行单变量时间序列预测

首先,时间序列预测问题是一个复杂的预测模型问题,它不像一般的回归预测模型。时间序列预测的输入变量是一组按时间顺序的数字序列。它既具有延续性又具有随机性,所以在建模难度上相对回归预测更大。

但同时,正好有一种强大的神经网络适合处理这种存在依赖关系的序列问题:RNN(Recurrent neural networks)。在过去几年中,应用 RNN 在语音识别,语言建模,翻译,图片描述等问题上已经取得一定成功,并且应用领域还在扩展。

LSTM网络

Long Short-Term Memory 网络亦称LSTM 网络,是一种在深度学习中应用的循环神经网络。可以学习长期依赖信息。LSTM 由Hochreiter & Schmidhuber (1997)提出,并在近期被Alex Graves进行了改良和推广。在很多问题,LSTM 都取得相当巨大的成功,并得到了广泛的使用。LSTM 通过刻意的设计来避免长期依赖问题。记住长期的信息在实践中是 LSTM 的默认行为,而非需要付出很大代价才能获得的能力。

具体应用

下面以一个洗发水销售的例子,来实现LSTM。 首先,你可以在这里下载到本文需要用的数据集。这是一个描述了3年内洗发水的月度销售数量的数据集。

数据读取

1
2
3
4
5
6
7
8
9
10
11
from pandas import read_csv
from pandas import datetime
from matplotlib import pyplot as plt
def parser(x):
    return datetime.strptime('190'+x,"%Y-%m")
series = read_csv('sales-of-shampoo-over-a-three-ye.csv',
                 header=0,parse_dates=[0],index_col=0,
                  squeeze=True, date_parser=parser)
print(series.head())
series.plot()
plt.show()

数据划分

首先我们把数据集划分成两个部分即:训练集和测试集。 那么我们该如何划分呢?因为我们今天研究的是时间序列分析,所以在数据集的划分上我们也应该按照时间来划分。我们可以将前两年的数据作为我们的训练集而将最后一年的数据作为测试集。

1
2
3
# split data into train and test
X = series.values
train, test = X[0:-12], X[-12:]

这里我们假设一个滚动预测的情景,又称前向模型验证(walk-forward model validation)。其原理很简单,举例来说就像当公司的预测期长达一年时,预测会将已过去的月份排除,而将预测期末的月份补上。好比一月份过去后,我们将其从预测中移除,同时次年的一月份就会作为收尾被添加到预测中以便预测总能保持12个月的完整性。

这样通过使用每月新的洗发水销售量来进行下个月的预测,我们就像模拟了一个更接近于真实世界的场景。

最后,我们将所有在测试集上的预测结果收集起来并计算出他们与真实值的均方根误差(RMSE)以此来作为评估我们模型的基准。

持续模型预测(Persistence Model Forecast)

持续性预测的基本思路就是从先前的(t-1)时间序列的结果用于预测当前时间(t)的取值。 那么根据以上的思路,我们可以通过滚动预测的原理从训练集的历史数据中获取最后一次观察值并使用它来预测当前时间的可能取值。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
from pandas import read_csv
from pandas import datetime
from sklearn.metrics import mean_squared_error
from math import sqrt
from matplotlib import pyplot
# load dataset
def parser(x):
  return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('sales-of-shampoo-over-a-three-ye.csv',
                  header=0, parse_dates=[0], index_col=0,
                  squeeze=True, date_parser=parser)
# split data into train and test
X = series.values
train, test = X[0:-12], X[-12:]
# walk-forward validation
history = [x for x in train]
predictions = list()
for i in range(len(test)):
  # make prediction
  predictions.append(history[-1])
  # observation
  history.append(test[i])
# report performance
rmse = sqrt(mean_squared_error(test, predictions))
print('RMSE: %.3f' % rmse)
# line plot of observed vs predicted
pyplot.plot(test)
pyplot.plot(predictions)
pyplot.show()

persistence_rmse.png

通过持续模型的预测,我们得到了一个最基础的预测模型以及RMSE(baseline)为了提升我们预测模型的效果,下面让我们进入正题来构建LSTM模型来对数据集进行时间序列预测。

数据处理

为了能够构建一个LSTM模型对训练集进行训练,我们首先要对数据进行一下处理:

  1. 将时间序列问题转化成监督学习问题
  2. 平稳时间序列
  3. 数据标准化

将时间序列转换成监督学习

对于一个时间序列问题,我们可以通过使用从最后一个(t-1)时刻的观测值作为输入的特征X和当前时刻(t)的观测值作为输出Y来实现转换。

因为,需要转换的是一组时间序列数据,所以无法组合成像真正的监督学习那样有明确一对一映射的输入输出关系。尤其是在数据集的最开始或最后时,两个位置总有一个位置无法在训练集中找到对应关系。为了解决这样的问题,我们通常的做法是,在最开始时将输入特征置为0,而它对应的输出就是时间序列的第一个元素。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
from pandas import read_csv
from pandas import datetime
from pandas import DataFrame
from pandas import concat

# frame a sequence as a supervised learning problem
def timeseries_to_supervised(data, lag=1):
  df = DataFrame(data)
  columns = [df.shift(i) for i in range(1, lag+1)]
  columns.append(df)
  df = concat(columns, axis=1)
  df.fillna(0, inplace=True)
  return df

# load dataset
def parser(x):
  return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('sales-of-shampoo-over-a-three-ye.csv', header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser)
# transform to supervised learning
X = series.values
supervised = timeseries_to_supervised(X, 1)
print(supervised)

输出结果: 0 0 0 0.0 266.0 1 266.0 145.9 2 145.9 183.1 3 183.1 119.3 4 119.3 180.3 5 180.3 168.5 6 168.5 231.8 7 231.8 224.5 8 224.5 192.8 9 192.8 122.9 10 122.9 336.5 11 336.5 185.9 12 185.9 194.3 13 194.3 149.5 14 149.5 210.1 15 210.1 273.3 16 273.3 191.4 17 191.4 287.0 18 287.0 226.0 19 226.0 303.6 20 303.6 289.9 21 289.9 421.6 22 421.6 264.5 23 264.5 342.3 24 342.3 339.7 25 339.7 440.4 26 440.4 315.9 27 315.9 439.3 28 439.3 401.3 29 401.3 437.4 30 437.4 575.5 31 575.5 407.6 32 407.6 682.0 33 682.0 475.3 34 475.3 581.3 35 581.3 646.9

平稳时间序列

虽然不明显,但我们仍可以看出这个洗发水销售数据集在时间上呈上升趋势。因此我们说这个时间序列数据是非平稳的。那么,不平稳怎么办?

答案就是:差分。(有关差分的介绍点击此处

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
from pandas import read_csv
from pandas import datetime
from pandas import Series

# create a differenced series
def difference(dataset, interval=1):
  diff = list()
  for i in range(interval, len(dataset)):
      value = dataset[i] - dataset[i - interval]
      diff.append(value)
  return Series(diff)

# invert differenced value
def inverse_difference(history, yhat, interval=1):
  return yhat + history[-interval]

# load dataset
def parser(x):
  return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('sales-of-shampoo-over-a-three-ye.csv', header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser)
print(series.head())
# transform to be stationary
differenced = difference(series, 1)
print(differenced.head())
# invert transform
inverted = list()
for i in range(len(differenced)):
  value = inverse_difference(series, differenced[i], len(series)-i)
  inverted.append(value)
inverted = Series(inverted)
print(inverted.head())
differenced.plot()
plt.show()

diff.png

经过一阶差分处理后,从图上看还是挺平稳的。

标准化数据

在数据输入前进行标准化可以非常有效的提升收敛速度和效果。尤其如果我们的激活函数是sigmoid或者tanh,其梯度最大的区间是0附近,当输入值很大或者很小的时候,sigmoid或者tanh的变化就基本平坦了(sigmoid的导数sig(1-sig)会趋于0),也就是进行梯度下降进行优化的时候,梯度会趋于0,而倒是优化速度很慢。

如果输入不进行归一化,由于我们初始化的时候一般都是0均值的的正太分布或者小范围的均匀分布(Xavier),如果输入中存在着尺度相差很大的特征,例如(10000,0.001)这样的,很容易导致激活函数的输入w1x1+w2x2+b变的很大或者很小,从而引起梯度趋于0。

而LSTM的默认激活函数就是tanh函数,它的输出范围在-1 到 1 之间,同时这是时间序列数据的首选范围。因此我们可以使用MinMaxScaler类将数据集转换到范围[-1,1]。像其他scikit用于转换数据的方法类一样,它需要以行和列的矩阵格式提供的数据。因此,在转换之前,我们必须重塑NumPy数组。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
from pandas import read_csv
from pandas import datetime
from pandas import Series
from sklearn.preprocessing import MinMaxScaler
# load dataset
def parser(x):
  return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('sales-of-shampoo-over-a-three-ye.csv', header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser)
print(series.head())
# transform scale
X = series.values
X = X.reshape(len(X), 1)
scaler = MinMaxScaler(feature_range=(-1, 1))
scaler = scaler.fit(X)
scaled_X = scaler.transform(X)
scaled_series = Series(scaled_X[:, 0])
print(scaled_series.head())
# invert transform
inverted_X = scaler.inverse_transform(scaled_X)
inverted_series = Series(inverted_X[:, 0])
print(inverted_series.head())

输出结果: Month 1901-01-01 266.0 1901-02-01 145.9 1901-03-01 183.1 1901-04-01 119.3 1901-05-01 180.3 Name: Sales of shampoo over a three year period, dtype: float64 0 -0.478585 1 -0.905456 2 -0.773236 3 -1.000000 4 -0.783188 dtype: float64 0 266.0 1 145.9 2 183.1 3 119.3 4 180.3 dtype: float64

构建LSTM模型

长短期记忆网络(LSTM)是一种递归神经网络(RNN)。 这类网络的的优点是它能学习并记住较长序列,并不依赖预先指定的窗口滞后观察值作为输入。 在Keras中,这被称为stateful,在定义LSTM网络层时将“stateful”语句设定为“True”。

LSTM层要求输入矩阵格式为:[样本,时间步长,特征]

鉴于训练数据集的形式定义为X输入和y输出,必须先将其转化为样本/时间步长/特征的形式。

完整代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
from pandas import DataFrame
from pandas import Series
from pandas import concat
from pandas import read_csv
from pandas import datetime
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import MinMaxScaler
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM
from math import sqrt
from matplotlib import pyplot
import numpy

# date-time parsing function for loading the dataset
def parser(x):
  return datetime.strptime('190'+x, '%Y-%m')

# frame a sequence as a supervised learning problem
def timeseries_to_supervised(data, lag=1):
  df = DataFrame(data)
  columns = [df.shift(i) for i in range(1, lag+1)]
  columns.append(df)
  df = concat(columns, axis=1)
  df.fillna(0, inplace=True)
  return df

# create a differenced series
def difference(dataset, interval=1):
  diff = list()
  for i in range(interval, len(dataset)):
      value = dataset[i] - dataset[i - interval]
      diff.append(value)
  return Series(diff)

# invert differenced value
def inverse_difference(history, yhat, interval=1):
  return yhat + history[-interval]

# scale train and test data to [-1, 1]
def scale(train, test):
  # fit scaler
  scaler = MinMaxScaler(feature_range=(-1, 1))
  scaler = scaler.fit(train)
  # transform train
  train = train.reshape(train.shape[0], train.shape[1])
  train_scaled = scaler.transform(train)
  # transform test
  test = test.reshape(test.shape[0], test.shape[1])
  test_scaled = scaler.transform(test)
  return scaler, train_scaled, test_scaled

# inverse scaling for a forecasted value
def invert_scale(scaler, X, value):
  new_row = [x for x in X] + [value]
  array = numpy.array(new_row)
  array = array.reshape(1, len(array))
  inverted = scaler.inverse_transform(array)
  return inverted[0, -1]

# fit an LSTM network to training data
def fit_lstm(train, batch_size, nb_epoch, neurons):
  X, y = train[:, 0:-1], train[:, -1]
  X = X.reshape(X.shape[0], 1, X.shape[1])
  model = Sequential()
  model.add(LSTM(neurons, batch_input_shape=(batch_size, X.shape[1], X.shape[2]), stateful=True))
  model.add(Dense(1))
  model.compile(loss='mean_squared_error', optimizer='adam')
  for i in range(nb_epoch):
      model.fit(X, y, epochs=1, batch_size=batch_size, verbose=0, shuffle=False)
      model.reset_states()
  return model

# make a one-step forecast
def forecast_lstm(model, batch_size, X):
  X = X.reshape(1, 1, len(X))
  yhat = model.predict(X, batch_size=batch_size)
  return yhat[0,0]

# load dataset
series = read_csv('sales-of-shampoo-over-a-three-ye.csv', header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser)

# transform data to be stationary
raw_values = series.values
diff_values = difference(raw_values, 1)

# transform data to be supervised learning
supervised = timeseries_to_supervised(diff_values, 1)
supervised_values = supervised.values

# split data into train and test-sets
train, test = supervised_values[0:-12], supervised_values[-12:]

# transform the scale of the data
scaler, train_scaled, test_scaled = scale(train, test)

# fit the model
lstm_model = fit_lstm(train_scaled, 1, 3000, 4)
# forecast the entire training dataset to build up state for forecasting
train_reshaped = train_scaled[:, 0].reshape(len(train_scaled), 1, 1)
lstm_model.predict(train_reshaped, batch_size=1)

# walk-forward validation on the test data
predictions = list()
for i in range(len(test_scaled)):
  # make one-step forecast
  X, y = test_scaled[i, 0:-1], test_scaled[i, -1]
  yhat = forecast_lstm(lstm_model, 1, X)
  # invert scaling
  yhat = invert_scale(scaler, X, yhat)
  # invert differencing
  yhat = inverse_difference(raw_values, yhat, len(test_scaled)+1-i)
  # store forecast
  predictions.append(yhat)
  expected = raw_values[len(train) + i + 1]
  print('Month=%d, Predicted=%f, Expected=%f' % (i+1, yhat, expected))

#report performance
rmse = sqrt(mean_squared_error(raw_values[-12:], predictions))
print('Test RMSE: %.3f' % rmse)
# line plot of observed vs predicted
pyplot.plot(raw_values[-12:])
pyplot.plot(predictions)
pyplot.show()

lstm_pred.png

最后运行结果打印出测试数据集12个月份中每个月份的预期和预测销量。示例还打印了所有预测值得均方根误差。该模型显示洗发水月度销量的均方根误差为111.925,好于持续性模型得出的对应结果136.761。

另外,神经网络的一个难题是初始条件不同,它们给出结果就不同。一种解决办法是修改Keras使用的随机数种子值以确保结果可复制。另一种办法是使用不同的实验设置控制随机初始条件。

如何免费使用谷歌GPU训练神经网络



完全云端运行:免费使用谷歌GPU训练神经网络

背景

对,你没有听错,高大上的GPU,现在不花钱也能用上了。这是Google的一项免费云端机器学习服务,全名Colaboratory。

Colaboratory 是一个 Google 研究项目,旨在帮助传播机器学习培训和研究成果。它是一个 Jupyter 笔记本环境,不需要进行任何设置就可以使用,并且完全在云端运行。Colaboratory 笔记本存储在 Google 云端硬盘中,并且可以共享,就如同您使用 Google 文档或表格一样。Colaboratory 可免费使用,而且最重要的还提供免费的英伟达Tesla K80 GPU。还有这等好事?事不宜迟,本文马上介绍如何使用 Google CoLaboratory 训练神经网络。

准备工作

在Google Drive上创建文件夹

Colab用的数据都存储在Google Drive云端硬盘上,所以,我们需要先指定要在Google Drive上用的文件夹。

比如说,可以在Google Drive上创建一个“app”文件夹,或者其他什么名字,也可以选择Colab笔记本默认的文件夹。

新建Colab笔记本

在刚刚创建的app文件夹里点击右键,选择“更多”,然后从菜单里选择“Colaboratory”,这样就新建出了一个Colab笔记本。

若是更多选项中没有“Colaboratory”选项,可以点击“关联更多应用”选项,然后在打开的页面中,搜索“Colaboratory”,然后再点关联应用,再次点击右键就可以在“更多”选项中看到“Colaboratory”选项了。

设置免费GPU

新建Colaboratory成功后,在笔记本里点Edit>Notebook settings(编辑>笔记本设置),或者Runtime>Change runtime type(运行时>改变运行时类型),然后在Hardware accelerator(硬件加速器)一栏选择GPU。

然后,Google Colab就可以用了。

关联Google Drive

为了能让Colaboratory使用到你的Google Drive的文件,我们需要先运行下面这些代码,来安装必要的库、执行授权。

!apt-get install -y -qq software-properties-common python-software-properties module-init-tools
!add-apt-repository -y ppa:alessandro-strada/ppa 2>&1 > /dev/null
!apt-get update -qq 2>&1 > /dev/null
!apt-get -y install -qq google-drive-ocamlfuse fuse
from google.colab import auth
auth.authenticate_user()
from oauth2client.client import GoogleCredentials
creds = GoogleCredentials.get_application_default()
import getpass
!google-drive-ocamlfuse -headless -id={creds.client_id} -secret={creds.client_secret} < /dev/null 2>&1 | grep URL
vcode = getpass.getpass()
!echo {vcode} | google-drive-ocamlfuse -headless -id={creds.client_id} -secret={creds.client_secret}

运行的时候应该会看到下图所示的结果:

看见那个链接之后,点击它,复制验证码并粘贴到文本框里。(这里其实是调用了Google Drive的SDK来访问你的Google Drive,而这个验证码就相当于access_key了)

授权完成后,就可以挂载Google Drive了:

!mkdir -p drive
!google-drive-ocamlfuse drive

测试GPU

这时,我们在本地电脑上创建一个.py文件来测试一下,挂载是否成功以及GPU是否在工作吧。

echo "import tensorflow as tf\nprint(tf.test.gpu_device_name())" > test.py

然后将test.py上传到我们开始时创建的app的文件夹里。

然后在Colaboratory笔记本中运行一下代码:

!python3 drive/app/test.py

不出意外的话,就会输出类似以下的结果:

/usr/local/lib/python3.6/dist-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  from ._conv import register_converters as _register_converters
2018-02-18 12:37:05.172726: I tensorflow/stream_executor/cuda/cuda_gpu_executor.cc:898] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero
2018-02-18 12:37:05.172988: I tensorflow/core/common_runtime/gpu/gpu_device.cc:1208] Found device 0 with properties: 
name: Tesla K80 major: 3 minor: 7 memoryClockRate(GHz): 0.8235
pciBusID: 0000:00:04.0
totalMemory: 11.17GiB freeMemory: 503.62MiB
2018-02-18 12:37:05.173016: I tensorflow/core/common_runtime/gpu/gpu_device.cc:1308] Adding visible gpu devices: 0
2018-02-18 12:37:05.457665: I tensorflow/core/common_runtime/gpu/gpu_device.cc:989] Creating TensorFlow device (/device:GPU:0 with 243 MB memory) -> physical GPU (device: 0, name: Tesla K80, pci bus id: 0000:00:04.0, compute capability: 3.7)
/device:GPU:0

到这里的话,那么恭喜你,你的GPU环境基本可以用了,只要把你的项目文件夹上传到你的app文件夹下,搭建好深度学习的库环境,就可以通过类似上面的操作进行神经网络训练了。

Tips

如何安装库?

安装Keras:

!pip install -q keras
import keras

安装PyTorch:

!pip install -q http://download.pytorch.org/whl/cu75/torch-0.2.0.post3-cp27-cp27mu-manylinux1_x86_64.whl torchvision
import torch

安装OpenCV:

!apt-get -qq install -y libsm6 libxext6 && pip install -q -U opencv-python
import cv2

安装XGBoost:

!pip install -q xgboost==0.4a30
import xgboost

安装GraphViz:

!apt-get -qq install -y graphviz && pip install -q pydot
import pydot

安装7zip Reader:

!apt-get -qq install -y libarchive-dev && pip install -q -U libarchive
import libarchive

安装其他库:

用!pip install或者!apt-get install命令。

元胞自动机与生命游戏(Game of Life)



元胞自动机与《生命游戏》(Game of Life)

背景

笔者最近读到一篇交通流仿真的论文里,提到了一个挺有意思的模型 -- 元胞自动机

在好奇心的驱使之下查询了不少资料,所以今天就来跟大家来分享一下“元胞自动机”这个模型以及它和“康威《生命游戏》”的关系。

为了让话题更加有趣,我们先从《生命游戏》开始谈起。

什么是《生命游戏》?

生命游戏由英国数学家约翰·何顿·康威提出,它其实是一个零玩家游戏,它包括一个二维矩形世界,这个世界中的每个方格居住着一个活着的或死了的细胞。

而整个《生命游戏》是贯彻着一条生命游戏定律的,即:如果一个生命,其周围的同类生命太少,会因为得不到帮助而死亡;如果太多,则会因为得不到足够的生命资源而死亡。 ——英国数学家约翰·康威

一个细胞在下一个时刻生死取决于相邻八个方格中活着的或死了的细胞的数量。如果相邻方格活着的细胞数量过多,这个细胞会因为资源匮乏而在下一个时刻死去;相反,如果周围活细胞过少,这个细胞会因太孤单而死去。实际中,你可以设定周围活细胞的数目怎样时才适宜该细胞的生存。如果这个数目设定过低,世界中的大部分细胞会因为找不到太多的活的邻居而死去,直到整个世界都没有生命;如果这个数目设定过高,世界中又会被生命充满而没有什么变化。

实际中,这个数目一般选取2或者3;这样整个生命世界才不至于太过荒凉或拥挤,而是一种动态的平衡。

这样的话,游戏的规则就是:当一个方格周围有2或3个活细胞时,方格中的活细胞在下一个时刻继续存活;即使这个时刻方格中没有活细胞,在下一个时刻也会“诞生”活细胞。在这个游戏中,还可以设定一些更加复杂的规则,例如当前方格的状况不仅由父一代决定,而且还考虑祖父一代的情况。你还可以作为这个世界的上帝,随意设定某个方格细胞的死活,以观察对世界的影响。

在游戏的进行中,杂乱无序的细胞会逐渐演化出各种精致、有形的结构;这些结构往往有很好的对称性,而且每一代都在变化形状。一些形状已经锁定,不会逐代变化。

有时,一些已经成形的结构会因为一些无序细胞的“入侵”而被破坏。但是形状和秩序经常能从杂乱中产生出来。

《生命游戏》和元胞自动机的关系以及什么是元胞自动机?

生命游戏的原理就是基于元胞自动机的,或者说《生命游戏》就是元胞自动机的一个展示。

元胞自动机(Cellular Automata,简称CA,也有人译为细胞自动机、点格自动机、分子自动机或单元自动机)。是一时间和空间都离散的动力系统。

散布在规则格网 (Lattice Grid)中的每一元胞(Cell)取有限的离散状态,遵循同样的作用规则,依据确定的局部规则作同步更新。大量元胞通过简单的相互作用而构成动态系统的演化。不同于一般的动力学模型,元胞自动机不是由严格定义的物理方程或函数确定,而是用一系列模型构造的规则构成。凡是满足这些规则的模型都可以算作是元胞自动机模型。

因此,元胞自动机是一类模型的总称,或者说是一个方法框架。其特点是时间、空间、状态都离散,每个变量只取有限多个状态,且其状态改变的规则在时间和空间上都是局部的。

初等元胞自动机( Elementary Cellular Automata, ECA)的基本要素如下

  • 空间:一维直线上等间距的点。可为某区间上的整数点的集合。
  • 状态集:S={s1,s2} 即只有两种不同的状态。这两种不同的状态可将其分别编码为0 与 1;若用图形表示,则可对应“黑”与“白” 或者其他两种不同的颜色。
  • 邻居:取邻居半径r=1,即每个元胞最多只有“左邻右舍”两个邻居。
  • 演化规则:任意设定, 最多2^8=256

元胞以相邻的8个元胞为邻居。即Moore邻居;一个元胞的生死由其在该时刻本身的生死状态和周围八个邻居的状态。

为了解释它,我们可以把计算机中的宇宙想象成是一堆方格子构成的封闭空间,尺寸为N的空间就有N*N个格子。而每一个格子都可以看成是一个生命体,每个生命都有生和死两种状态,如果该格子生就显示蓝色,死则显示白色。每一个格子旁边都有邻居格子存在,如果我们把3x3的9个格子构成的正方形看成一个基本单位的话,那么这个正方形中心的格子的邻居就是它旁边的8个格子。

每个格子的生死遵循下面的原则:

  1. 如果一个细胞周围有3个细胞为生(一个细胞周围共有8个细胞),则该细胞为生(即该细胞若原先为死,则转为生,若原先为生,则保持不变)。
  2. 如果一个细胞周围有2个细胞为生,则该细胞的生死状态保持不变;
  3. 在其它情况下,该细胞为死(即该细胞若原先为生,则转为死,若原先为死,则保持不变)

设定图像中每个像素的初始状态后依据上述的游戏规则演绎生命的变化,由于初始状态和迭代次数不同,将会得到令人叹服的优美图案。

这样就把这些若干个格子(生命体)构成了一个复杂的动态世界。运用简单的3条作用规则构成的群体会涌现出很多意想不到的复杂行为,这就是复杂性科学的研究焦点。

元胞自动机的应用

在实际应用过程中,有的元胞自动机模型对其中的某些特征进行了扩展,有的在规则设计中引入随机因素,如:森林火灾模型。 又如,在交通、通讯发达的今天, 研究流行病或计算机病毒的传播问题时, 我们还可以将空间背景换成复杂网络的结点,用网络邻接点作为邻居。

这样的调整显然比仍旧使用二维欧氏空间、采用欧氏距离的模型更加符合实际情况。 在大型场所人群紧急疏散问题模拟研究中,可以考虑年龄、性别等因素,即元胞不是同质的,更加有利于使模拟系统接近真实系统。

《生命游戏实现》Python版

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
import numpy as np
import matplotlib.pyplot as plt


class GameOfLife(object):

  def __init__(self, cells_shape):
    """
    Parameters
    ----------
    cells_shape : 一个元组,表示画布的大小。

    Examples
    --------
    建立一个高20,宽30的画布
    game = GameOfLife((20, 30))
    
    """

    # 矩阵的四周不参与运算
    self.cells = np.zeros(cells_shape)

    real_width = cells_shape[0] - 2
    real_height = cells_shape[1] - 2

    self.cells[1:-1, 1:-1] = np.random.randint(2, size=(real_width, real_height))
    self.timer = 0
    self.mask = np.ones(9)
    self.mask[4] = 0

  def update_state(self):
    """更新一次状态"""
    buf = np.zeros(self.cells.shape)
    cells = self.cells
    for i in range(1, cells.shape[0] - 1):
      for j in range(1, cells.shape[0] - 1):
        # 计算该细胞周围的存活细胞数
        neighbor = cells[i-1:i+2, j-1:j+2].reshape((-1, ))
        neighbor_num = np.convolve(self.mask, neighbor, 'valid')[0]
        if neighbor_num == 3:
          buf[i, j] = 1
        elif neighbor_num == 2:
          buf[i, j] = cells[i, j]
        else:
          buf[i, j] = 0
    self.cells = buf
    self.timer += 1

  def plot_state(self):
    """画出当前的状态"""
    plt.title('Iter :{}'.format(self.timer))
    plt.imshow(self.cells)
    plt.show()

  def update_and_plot(self, n_iter):
    """更新状态并画图
    Parameters
    ----------
    n_iter : 更新的轮数
    """
    plt.ion()
    for _ in range(n_iter):
      plt.title('Iter :{}'.format(self.timer))
      plt.imshow(self.cells)
      self.update_state()
      plt.pause(0.2)
    plt.ioff()


if __name__ == '__main__':
  game = GameOfLife(cells_shape=(60, 60))
  game.update_and_plot(200)

效果图 



使用树莓派实现24小时不间断直播



使用树莓派进行24小时不间断直播

开始

多余的话就不多说了,今天本文为大家介绍两种使用树莓派来做直播服务器的方法。

方案一 ffmpeg + ffserver搭建流媒体服务器

首先
我们用到的工具有:

硬件方面:

  • 树莓派主板一块
  • 兼容树莓派的USB摄像头一个

软件方面:

  • ffmpeg,负责媒体文件的转码工作,把你服务器上的源媒体文件转成要发出去的流媒体文件。
  • ffserver,负责响应客户端的流媒体请求,把流媒体数据发送给客户端,相当与一个小型的服务端。

具体的工作方式就如下图所示:

多个输入源被“喂”到广播服务器,这些多媒体内容就会分发到多个客户端。上图的目的是显示地表明你的流系统能够被分成多个块部署到网络上,允许你广播不同的在线内容,而不需要改变流媒体系统的结构。

配置

无论是树莓派官方摄像头模块还是其他兼容的USB摄像头,连接好摄像头之后,运行命令去启用摄像头:

sudo raspi-config

ffserver.conf,ffserver启动时的配置文件,在这个文件中主要是对网络协议,缓存文件feed1.ffm(见下述)和要发送的流媒体文件的格式参数做具体的设定。

feed1.ffm,可以看成是一个流媒体数据的缓存文件,ffmpeg把转码好的数据发送给ffserver,如果没有客户端连接请求,ffserver把数据缓存到该文件中。

下面就是一个ffserver.conf的一个例子

Port 8090                      # Port to bind the server to
BindAddress 0.0.0.0
MaxHTTPConnections 2000
MaxClients 1000
MaxBandwidth 10000             # Maximum bandwidth per client
                               # set this high enough to exceed stream bitrate
CustomLog -
NoDaemon                       # Remove this if you want FFserver to daemonize after start

<Feed feed1.ffm>               # This is the input feed where FFmpeg will send
   File ./feed1.ffm            # video stream.
   FileMaxSize 64M              # Maximum file size for buffering video
   ACL allow 127.0.0.1         # Allowed IPs
</Feed>

<Stream test.webm>              # Output stream URL definition
   Feed feed1.ffm              # Feed from which to receive video
   Format webm

   # Audio settings
   AudioCodec vorbis
   AudioBitRate 64             # Audio bitrate

   # Video settings
   VideoCodec libvpx
   VideoSize 720x576           # Video resolution
   VideoFrameRate 25           # Video FPS
   AVOptionVideo flags +global_header  # Parameters passed to encoder
                                       # (same as ffmpeg command-line parameters)
   AVOptionVideo cpu-used 0
   AVOptionVideo qmin 10
   AVOptionVideo qmax 42
   AVOptionVideo quality good
   AVOptionAudio flags +global_header
   PreRoll 15
   StartSendOnKey
   VideoBitRate 400            # Video bitrate
</Stream>

<Stream status.html>            # Server status URL
   Format status
   # Only allow local people to get the status
   ACL allow localhost
   ACL allow 192.168.0.0 192.168.255.255
</Stream>

<Redirect index.html>    # Just an URL redirect for index
   # Redirect index.html to the appropriate site
   URL http://www.ffmpeg.org/
</Redirect>

ffserver启动时默认查看 /etc/ffserver.conf 配置文件,你可以通过-f选项控制查阅的配置文件。

ffserver -f ffserver.conf

运行结果如下所示的话,那么ffserver就算是启动成功了。

打开http://localhost:8090/status.html可以看到当前server中各个流的状态。

接入视频流

ffserver启动之后,就可以向
http://localhost:8090/feed1.ffm接入视频流。

注意,这里不需要指定编码格式,FFserver会重新编码。

视频流的来源可以是文件、摄像头或者录制屏幕。

接入视频文件

ffmpeg -i testvideo.mp4 http://localhost:8090/feed1.ffm

接入录制屏幕

ffmpeg -f x11grab -r 25 -s 640x512 -i :0.0 -f alsa -i pulse http://localhost:8090/feed1.ffm

这里有两个-f,第一个指的是视频流,第二个指的是音频流。视频流是抓取屏幕形成视频,-r设置帧率为25帧/s,-s设置抓取图像大小为640x512,-i设置录制视频的初始坐标。音频流设置为alsa(Advanced Linux Sound Architecture),从Linux系统中获取音频。这其中这样ffmpeg可以录制屏幕feed到feed1.ffm中。

接入摄像头直播

ffmpeg -f video4linux2 -s 640x480 -r 25 -i /dev/video0 -f alsa -i pulse http://localhost:8090/feed1.ffm

方案二 avconv 和 GStreamer 用于采集摄像头捕获的视频流并推送到 RTMP 服务

首先
我们用到的工具有:

硬件方面:

  • 树莓派主板一块
  • 兼容树莓派的USB摄像头一个

软件方面:

  • avconv 和 GStreamer 用于采集摄像头捕获的视频流并推送到 RTMP 服务
  • NGINX 和 RTMP 模块,用于接收视频流,同时提供视频发布功能

安装&配置

因为这里我们要用到nginx的rtmp模块作为服务端,而系统自带的apt安装的nginx是没有这个模块的,所以我们需要先安装后移除nginx然后再手动编译(安装是为了下载好相关依赖)。

sudo apt-get update
#安装 nginx
sudo apt-get -y install nginx
#移除 nginx
sudo apt-get -y remove nginx
sudo apt-get clean
#清空 nginx 的配置文件
sudo rm -rf /etc/nginx/*
#安装编译用的模块
sudo apt-get install -y curl build-essential libpcre3 libpcre3-dev libpcre++-dev zlib1g-dev libcurl4-openssl-dev libssl-dev
#创建存放网页的目录给 nginx 使用
sudo mkdir -p /var/www
#创建编译用的目录
mkdir -p ~/nginx_src
cd ~/nginx_src
#下载 nginx 源码包
wget http://nginx.org/download/nginx-1.11.8.tar.gz
#下载 nginx-rtmp-module 源码包
wget https://github.com/arut/nginx-rtmp-module/archive/master.zip
tar -zxvf nginx-1.11.8.tar.gz
unzip master.zip
cd nginx-1.11.8
#设定编译参数
./configure --prefix=/var/www --sbin-path=/usr/sbin/nginx --conf-path=/etc/nginx/nginx.conf --pid-path=/var/run/nginx.pid --error-log-path=/var/log/nginx/error.log --http-log-path=/var/log/nginx/access.log --with-http_ssl_module --without-http_proxy_module --add-module=/home/pi/nginx_src/nginx-rtmp-module-master
#开始编译安装
make
sudo make install

配置 nginx

sudo gedit /etc/nginx/nginx.conf

在末尾添加

rtmp {
    server {
        listen 1935;
        chunk_size 4096;
        application live {
            live on;
            record off;
        }
    }
}

重启 nginx 服务。

sudo service nginx start

安装 avconv 和 GStreamer

sudo apt-get update
sudo apt-get install libav-tools
#安装 GStreamer
sudo apt-get install gstreamer1.0-tools
#安装 GStreamer 扩展组件
sudo apt-get  install libgstreamer1.0-0 libgstreamer1.0-0-dbg libgstreamer1.0-dev liborc-0.4-0 liborc-0.4-0-dbg liborc-0.4-dev liborc-0.4-doc gir1.2-gst-plugins-base-1.0 gir1.2-gstreamer-1.0 gstreamer1.0-alsa gstreamer1.0-doc gstreamer1.0-omx gstreamer1.0-plugins-bad gstreamer1.0-plugins-bad-dbg gstreamer1.0-plugins-bad-doc gstreamer1.0-plugins-base gstreamer1.0-plugins-base-apps gstreamer1.0-plugins-base-dbg gstreamer1.0-plugins-base-doc gstreamer1.0-plugins-good gstreamer1.0-plugins-good-dbg gstreamer1.0-plugins-good-doc gstreamer1.0-plugins-ugly gstreamer1.0-plugins-ugly-dbg gstreamer1.0-plugins-ugly-doc gstreamer1.0-pulseaudio gstreamer1.0-tools gstreamer1.0-x libgstreamer-plugins-bad1.0-0 libgstreamer-plugins-bad1.0-dev libgstreamer-plugins-base1.0-0 libgstreamer-plugins-base1.0-dev

采集与呈现视频流

gst-launch-1.0 -v v4l2src device=/dev/video0 ! 'video/x-raw, width=1024, height=768, framerate=30/1' ! queue ! videoconvert ! omxh264enc ! h264parse ! flvmux ! rtmpsink location='rtmp://树莓派的IP地址/live live=1' &

采用以上命令就可以在后台采集USB摄像头拍摄的直播内容并推送到rtmp服务端上了。

呈现直播视频画面

1、使用 RTMP 播放器播放视频流
例如 VLC 等播放器(桌面版和手机版均有)支持 RTMP 视频流播放,填入 rtmp://树莓派的IP地址/live 即可播放。不过这个软件有数十秒的缓冲延迟,需要设定缓冲时间来缩短延迟。

2、推送至斗鱼直播平台观看
你可能注意到了 GStreamer 这个命令中有 location 这个参数。这个参数是设定采集到的视频流推向哪里,通过设定这个参数可以将视频流推向任何支持 RTMP 协议的服务器。

斗鱼平台同样采用了 RTMP 协议传输直播视频,首先获取斗鱼的 RTMP 推流地址。开启了直播室之后可以获得推流码。注意,斗鱼的推流码是有时限的,取到推流码需要尽快使用以免过期。

谈谈Google自动编程框架AutoML



谈谈Google自动编程框架AutoML

背景

首先,近年在Google AI First 的战略领导之下,Google 陆陆续续发布了不少AI相关产品。那么,我们今天就来看看最近许多公众号和报道中提到的AutoML系统。

据最新的报道:Google AutoML 系统自主编写机器学习代码,其效率在某种程度上竟然超过了专业的研发工程师。

但AutoML的目标并不是要将人类从开发过程中剥离出去,也不是要开发全新的人工智能,而是让人工智能继续维持某种速度来改变世界。

看完这篇报道后许多程序员开始担心未来程序员的工作将很快会被替代。但在笔者看来,AutoML的出现更像是机器学习领域中早就该出现的一个“编译器”,而对于被替代的担心,纯属无稽之谈。就像人写的汇编代码几乎都不可能好过由编程自动生成的优化后的汇编代码一样。

下面就让我们一起来看看AutoML到底是何方神圣吧。

什么是AutoML?

根据AutoML官网上的介绍:
机器学习(Machine Learning, ML)近年来取得了相当大的成功,越来越多的学科需要依赖它。然而,这个成功的关键是需要人类机器学习工程师完成以下的工作:

  • 预处理数据
  • 选择适当的功能
  • 选择一个适当的模型选择系列
  • 优化模型超参
  • 后处理机器学习模型
  • 严格分析所得的结果

由于这些任务的复杂性通常超过了非机器学习专家的能力,机器学习应用的快速增长产生了对于现成的机器学习方法的需求,而且这些现成的机器学习方法简单易使用且不需要专业的知识。我们称以机器学习的渐进自动化为目标的研究领域为AutoML(Automatic Machine Learning, AutoML)

​虽然它的最终用户面向那些没有专业机器学习知识的人,但AutoML依然向机器学习专业人士提供一些新的工具,如:

  • 执行深层表示的架构搜索
  • 分析超参数的重要性

遵循优化编程的范例,AutoML主张开发可以用数据驱动的方式自动实例化的灵活软件包。

AutoML的架构

AutoML网络的设计从卷积架构的初始版本进行多年的仔细实验和细化完成的。

在团队一个名为「AutoML」的项目中(如图所示),左边有一个名为「控制器」(the controller)的 RNN,它设计出一个「child」的模型架构,而后者能够通过某些特定任务进行训练与评估。随后,反馈的结果(feedback)得以返回到控制器中,并在下一次循环中提升它的训练设定。这一过程重复上千次——生成新的架构、测试、再把反馈输送给控制器再次学习。最终,控制器会倾向于设计那些在数据集中能获得更高准确性的架构,而反之亦然。

谷歌团队将这一方法应用于深度学习的两大数据集中,专注图像识别的 CIFAR-10 与语言建模的 Penn Treebank。在两个数据集上,系统自行设计的模型性能表现与目前机器学习专家所设计的领先模型不相上下。

让机器自行选择架构(machine-chosen architecture),与人类在设计神经网络的时候有一些共通之处,比如都采用了合并输入,并借鉴了此前的隐藏层。但其中也有一些亮点,比如机器选择的架构包含乘法组合 ( multiplicative combination),如右图最左边(机器设计)的蓝色标签为「elem_mult」。对于循环神经网络而言,出现组合的情况并不多见,可能因为人类研究者并没有发现明显的优势。有意思的地方在于,此前人类设计者也提议过机器采用的乘法组合,认为这种方法能够有效缓解梯度消失/爆炸问题。这也就意味着,机器选择的架构能够对发现新的神经架构大有裨益。

此外,机器还能教会人类为何某些神经网络的运行效果比较好。上图右边的架构有非常多的渠道,梯度可以向后流动,这也解释了为何 LSTM RNNs 的表现比标准 RNN 的性能要好。

「从长远看来,我们对于机器所设计的架构进行深入的分析和测试,这能够帮助我们重新定义原本自身对架构的看法。如果我们成功,这意味着将会启发新的神经网络的诞生,也能让一些非专家研究人员根据自己的需要创造神经网络,让机器学习造福每一个人。」

AutoML的案例应用

AutoML旨在创建可以由ML新手”开箱即用“的软件。

  • AutoWEAK:就是是一种可以同时选择机器学习算法和其对应超参数的方法;通过使用WEKA包,可以为各种数据集自动生成良好的模型。
  • auto-sklearn:则是一个自动化的机器学习工具包,甚至可以直接替换scikit-learn estimator模块。

测试结果

既然,AutoML说的这么厉害,笔者这里就使用auto-sklearn来测试下分类模型的效果。(使用MNIST数据集)
有关AutoML的具体使用教程,日后再放出,这里就先简单贴个结果,评测个模型准确率和性能。

首先,是笔者的LogisticRegression模型+特征工程:

LogisticRegression(PolynomialFeatures(input_matrix, degree=2, include_bias=False, interaction_only=False), C=0.5, dual=False, penalty=l2)
Accuracy score 0.988888888889

然后是auto-sklearn框架的自生成模型:

import autosklearn.classification
import sklearn.model_selection
import sklearn.datasets
import sklearn.metrics
X, y = sklearn.datasets.load_digits(return_X_y=True)
X_train, X_test, y_train, y_test = \
sklearn.model_selection.train_test_split(X, y, random_state=1)
automl = autosklearn.classification.AutoSklearnClassifier()
automl.fit(X_train, y_train)
y_hat = automl.predict(X_test)
print(“Accuracy score”, sklearn.metrics.accuracy_score(y_test, y_hat))
[Out] Accuracy score 0.993333333333

小结

看到结果,心情复杂。辛辛苦苦的特征工程果然还是干不过AutoML。不过,AutoML生成模型的耗时也是相当高的。
PS 说个题外话。微软一样的服务(不用写代码,不用调参数,会拖控件就能帮你训练深度学习模型)已经发布快一年了。
网站:https://www.customvision.ai/
新闻报道:https://thenextweb.com/dd/2017/05/10/microsofts-custom-vision-lets-build-computer-vision-ai-models-minutes/
只能说Google不愧是IT界的”广告公司“,这公关的对决高下立判。