2010年9月28日星期二

基本的蒙特卡洛技术(二)

将计算广告学中的一个问题形式化:

假设你开办了一个类似Hulu的网站,有n个视频,上周视频k的访问量是x[k],x[k]>0,通过某种时间序列模型,你预测下周视频k的访问量是y[k],y[k]>=0。

记x=sum{x[k],k=1..n}​, y=sum{y[k],k=1..n}

假设你现在有上周的全部x个视频请求,从中采样m个,第i个代表w[i]个,使得它们代表了下周的视频请求,即E(sum{w[i]h(i,k), i=1..m})=y[k]对任意k=1..n成立,其中

示性函数h(i,k)=1当且仅当样本i请求视频k。​

1. 加权采样

算法:

有放回、等概率地从上周x个视频请求中随机选取m个,如果选中的视频请求i来自视频k,则令w[i]=c*y[k]/x[k],其中c是归一化系数,等于x/m。

无偏性:

注意到E(w[i]h(i,k))=w[i]P(h(i,k)=1)=x/m*y[k]/x[k]*x[k]/x=y[k]/m

停时:

T显然概率一地为m

2. 拒绝方法

算法:

重复下列步骤直到接受m个请求为止,有放回、等概率地从上周x个视频请求中随机选取1个,假设它来自视频k,则以概率r[k]接受之,其中r[k]=c*y[k]/x[k],其中c=1/max{y[k]/x[k], k=1..n},接受的请求的权重w都是y/m

无偏性:

只需证明P(h(i,k)=1)=y[k]/y

令事件S(i,k)表示从历史采中的请求i来自视频k,S表示sample,事件A(i)表示接受采中的请求i,A表示Accept。我们已知P(S(i,k))=x[k]/x, P(A(i)|S(i,k))=r[k]

注意到h(i,k)=1是说,在接受的请求i的条件下,i来自视频k,因此转化成计算条件概率

P(h(i,k)=1)=P(S(i,k)|A(i))=P(A(i)|S(i,k))*P(S(i,k))/sum{P(A(i)|S(i,k))*P(S(i,k)), k=1..n}

注:​P(h(i,k)=1)=y[k]/y表明,拒绝方法可以达到按照目标分布采样的目的

停时:

即接受m个请求,需要尝试T次采样。E(T)=m/P(A(i))。接受概率P(A(i))=y/x/max{y[k]/x[k], k=1..n},也就是说,要想接受概率不太低,那么历史分布和目标分布必须一致接近。

3. 拒绝方法的拉斯维加斯版本

算法:

重复下列步骤直到视频k上的请求采够y[k]/w个,w=y/m。有放回、等概率地从上周x个视频请求中随机选取1个,假设它来自视频k,接受它当且仅当视频k上的请求还没采够y[k]/w个。接受的请求权重都是w。

注:拉斯维加斯版本的好处是不用计算接受概率,即不需要知道x[k]的情况,这一点在“历史分布”难以计算的情况下尤其有用。

正确性:

显然算法停止时,P(sum{w[i]h(i,k), i=1..m}=y[k])=1。

停时:

同前

4. 混合加权拒绝方法

加权方法完全靠样本权重将历史分布拉到目标分布,以牺牲权重相近为代价,换得最有停时性能。

拒绝方法完全靠接受拒绝控制了分布,以牺牲停时为代价,换得相近的权重。

有些场合,例如利用样本做某种预算控制的模拟,要求权重差异不能太大。这种混合方法目的是在停时和权重差异之间寻求平衡。

算法:

条件接受概率r[k]=min{1, c*y[k]/x[k]},c>1/max{y[k]/x[k], k=1..n}。

权重可以反算出来,当h(i,k)=1时,w[i]=y[k]/m*q/q[k],其中q=sum{q[k]},q[k]=r[k]*x[k]/x,这里q其实就是接受概率。

停时

注意到c越大,接受概率越大,停时越小,直到退化为纯加权方法。

5. 停时有界的混合加权拒绝方法

实际应用,由于SLA等原因,我们要求算法概率1地保证停时有界,即存在M,使得P(T<=M)=1。

这时我们松弛掉“恰好采纳m个样本”的要求。这个方法可以看做方法4的拉斯维加斯版本。

算法:

按照历史分布采M个,M>m,权重都是y/m,采够了的视频,拒绝更多样本,最后没采够的视频,通过调整其上每个样本权重,使之满足y[k]。

注:

这个算法不仅概率1保证停时有界,而且概率1保证满足等式约束sum{w[i]h(i,k), i=1..m}=y[k],代价是无法预知采纳的样本的个数。

2010年9月24日星期五

Start to review a textbook with practice

2 years ago, when I was a graduate student, the lecture of Monte Carlo Method interested me because I needed it for the RBM model in my thesis, although the examples provided by Prof Wang in the class were all about finance which I don't interest in.

The book I prepare to review in Chinese (not translate only) is

Monte Carlo Strategies in Scientific Computing

I want to mix in my practice in Computational Advertising where possible.

====标志开始的分割线====

评论+实践计划:

对应原书章节

主题

时间

2​

基本的蒙特卡洛技术

2010年9月、10月

3、4​

序贯蒙特卡洛

2010年11月、12月

5​

马尔科夫链蒙特卡洛

2011年1月、2月

6、7、8​

条件分布的采样

2011年3月、4月

9​

分子动力学与混合蒙特卡洛

2011年5月

10、11​

蒙特卡洛采样的性能优化

2011年6月、7月

12、13​

收敛性分析

2011年8月、9月



基本的蒙特卡洛技术(一)

蒙特卡洛常被理解为随机模拟算法,它的第一个问题就是如何生成满足某种分布的随机数,当没有办法直接从目标分布函数抽样,而只能从另一个分布函数抽样的时候。

1、​“拒绝方法”是基本的方法。它不保证每采一下都能采中一个能要的随机数,但是可以保证采中的要的数满足目标分布函数。

2、“加权采样”其实是退而求其次,它并没有去真的生成满足目标分布的随机数,而是假设该分布函数用来计算某个可积函数的期望,构造一个(随机点、权重)序列,使得该可积函数在这些点上的值的加权平均是原期望的无偏估计。

(未完待续)

2010年9月14日星期二

全序集的本质

上次在和同事讨论问题的时候发现,偏序集不能排序。后来我就很好奇为什么不能,可以排序的集合上应该至少定义一个满足什么条件的二元关系?
先回顾一下偏序集和全序集的定义:

给定集合S,“≤”是S上的二元关系,若“≤”满足:

   1. 自反性:∀a∈S,有a≤a;
   2. 反对称性:∀a,b∈S,a≤b且b≤a,则a=b;
   3. 传递性:∀a,b,c∈S,a≤b且b≤c,则a≤c;

则称“≤”是S上的非严格偏序或自反偏序。通常简称偏序

给定集合S,“≤”是S上的二元关系,若“≤”满足:
   1. 完全性 a ≤ b 或 b ≤ a
   2. 反对称性 如果 a ≤ b 且 b ≤ a 则 a = b
   3. 传递性如果 a ≤ b 且 b ≤ c 则 a ≤ c

则称“≤”是S上的全序

显然完全性蕴含了自反性,因此全序都是偏序。

命题1:如果“≤”集合S上的偏序,且满足4. 反向传递性:∀a,b,c∈S,a≤b不成立且b≤c不成立,则a≤c不成立”,那么“≤”集合S上的全序

证明:只需要证明“≤” 具有完全性。假设a ≤ b 或 b ≤ a 不成立,即a≤b不成立且b≤a不成立,由反向传递性,a≤a不成立,与自反性矛盾。

命题2:如果“≤”集合S上的全序,那么“≤” 满足反向传递性

证明:假设a≤b不成立且b≤c不成立由完全性,“a≤b不成立”蕴含了“b ≤ a ,“b≤c不成立”蕴含了“c≤b。再由传递性,ca。假设a≤c ,则由反对称性,得a=c。a≤b不成立,即c≤b不成立,而b≤c也不成立,与完全性矛盾。

结合上述命题,可知:一个偏序集合是全序集合,当且仅当具有反向传递性。也就是说,反向传递性在某种程度上是全序集合的本质。

经典排序算法中,如何隐含地使用了指定序的反向传递性来保证算法正确性呢?

注意到,排序算法总是要求传入一个"<",但是保证的是数组按照""排好序。其中ab定义为b<a不成立。
以快速排序算法为例,假设pivot为b,所有<b的元素放在b的左边,其他元素放在b的右边,不失一般性(递归地),只需证明"b左边的任一元素a""b右边的任一元素c"。

假设存在a和c,使得ac不成立。由a在b的左边,知a<b成立,即ba不成立。由反向传递性,bc不成立,即c<b,这与c在b的右边矛盾。























2010年9月13日星期一

URL size limit of Google Chart API

Google 帮助这样说
The longest URL that Google accepts in a chart GET request is 2048 characters in length, after URL-encoding (e.g., | becomes %7C). For POST, this limit is 16K.

按照帮助中的建议,将坐标数据压缩后,可以画更复杂(线条更稠密)的图了。
最简单的文本编码前缀是
chd=t:
后面接1,2,3.1415,2.17这样的字符串

压缩编码可以用下面这个函数实现单向转换:

def encoding_chart_data(l, upper_bound):#l is a list of uint
    enc_map = 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789'
    return "".join(map(lambda x:enc_map[int(round((len(enc_map)-1) * 1.0 * x / upper_bound))],l))

编码后的前缀是
chd=s:

调Google Char API的时候,很多时候感觉都像是在调一个矢量绘图的API,脑子里时刻要充满Shape/Scale的概念,而不是X,Y的欧氏坐标。很亲切。

有两种比较另类的Chart很吸引我:
http://code.google.com/apis/chart/docs/gallery/graphviz.html
可以将一个“图”可视化(GraphViz语法)
http://code.google.com/apis/chart/docs/gallery/formulas.html
渲染一个公式(TeX语法)

还提供了两个工具,相当有用:
http://code.google.com/apis/chart/docs/chart_wizard.html
向导,帮助你快速熟悉API的各项参数
http://code.google.com/apis/chart/docs/chart_playground.html
在线调试工具,帮助你缩进API URL,显示调试信息,例如告诉你为什么某一组参数无法渲染,非常好用