2010年12月29日星期三

评估在线广告的效果

很久以前看到这样一篇论文
Evaluating Online Ad Campaigns in a Pipeline: Causal Models At Scale
一直想看进去却没有整块时间。
终于在今年年底开年会的时候,在密云的一个安静的角落里看进去了。趁着还有点印象,写篇Blog。

假设广告商在网上投放广告的目的是提升品牌影响力,就像医院给病人用药的目的是让病人缓解病症。
证明广告有效,就像证明药品有效一样,需要"临床"数据。

考虑这样一个实验:
你有100个胃炎病人,程度不一,生死掌握在你手上。你有一种新药,随机给其中50个人吃,结果这些大多数都活了,剩下的50个人差不多都死了。于是你说,斯达舒,管用。

一个怀疑的声音说:人活并非因为药有效,而是因为心理作用:医院给我吃药,我有救了。
为此改进的实验:给不吃新药的人吃安慰剂:一种外观和新药相同的胶囊。

类比到此结束,你不能给targeting中的用户看一个没有任何品牌符号的广告(吃安慰剂),因为没有广告商愿意为了你的对照组而浪费Inventory。
更糟糕的是,你不能随便找一个没看过此广告的用户当作对照,因为有可能他压根儿就屏蔽了所有广告(好吧,我承认屏蔽视频广告不那么容易)。
所以,对照组的用户必须要看过其他广告才行,以表明他没有屏蔽广告。

即便如此,也还不够。另外一个不能类比的地方在于,吃什么药,是医院说了算,被实验的病人没有选择,吃完药病症如何自己也左右不了,生理现象。而看广告不一样,决定要不要看广告的是用户,决定看完广告去干什么的还是用户。你除了观测、分析,能控制的很少。文章中,选了15k个看了给定广告的用户,只统计他们在看这个广告之前的行为,同时选了70k个没看这个广告的用户,发现:前一组用户访问该品牌网站的次数的中位数,是后一组
用户访问该品牌网站的次数的中位数的4倍多。
如果因果联系必须是先行后继的,那么这样两组用户不能互相对照,因为他们活跃程度不同,至少在看广告之前就对此品牌的兴趣度不同。

有两个方案:
接受拒绝:扔掉对照组中的一部分不活跃用户,使得剩下的用户的活跃度和考察组在看广告之前的活跃度相同。实验发现,这样很可能失败因为对照组的用户太不活跃了,以致于几乎都被拒绝了。
调权重:通过简单乘上概率的倒数来降低对照组中不活跃用户的权重。虽然是无偏的,但是方差很大,以致于通常也不可用。
文章给出了一个鲁棒的估计,在调权重上花了很大功夫,也是最有趣的部分。实验结果发现,其实广告效果不是那么明显,想观测出来都难,更甭提评估了。也许这才体现了这篇文章的价值。

这种Research只有Google这样的公司能做,因为它有很多Tool Bar安装的用户,可以采集用户的很多行为数据。
我的怀疑:安装Google Tool Bar,是否能假设等同于全体用户的集体偏好?这恐怕是一个难以检验的假设。

2010年12月23日星期四

wirte C++ for MongoDB: smoking on Ubuntu

这周利用晚上的时间在笔记本的ubuntu上装了mongodb,把官网上的一小段C++ code点亮了。
参考
http://www.mongodb.org/pages/viewpage.action?pageId=133415
http://www.mongodb.org/display/DOCS/Building+for+Linux
怨念:
1)它依赖boost等,得先安装一大堆乱七八糟的库
2)源代码庞大,从github上全clone出来需要很久,然后scons又要build很久,都build完整个目录将近800MB,晕死
3)编译C++ code时候好像绕不开-L参数,不然link的时候报错,分明安装的时候所有东西都已经存到系统路径了,奇怪
g++ tut.cpp -o tut -lmongoclient -lboost_thread-mt -lboost_filesystem-mt -lboost_program_options-mt -L/home/ytwang/lib/mongo
我猜是因为libmongoclient.a在/home/ytwang/lib/mongo
而且引用应该写成这样
#include <mongo/client/dbclient.h>
4)运行时需要后台有守护进程mongod,像mysqld,而不是像sqlite那样纯粹基于文件
5)生成的binary也比较臃肿,刚点亮就快10MB了,感觉远不像想象中那样轻量级

2010年12月17日星期五

Proof of my CNF converter

定义:
一个句子是一个变量或一个满足下列条件的list,至少有3个元素,第一个元素是连词AND或OR,其余元素是句子(递归定义),这些句子称作子句。没有NOT这种连词。
称一个句子是CNF的,若它是一个变量或满足下列条件的list。第一个元素是AND,其余元素是变量或满足下列条件的list。第一个元素是OR,其余元素是变量。(非递归定义)
对偶地可以定义DNF。
称一个句子是melt的,若任何一个子句的连词都不同于这个句子的连词。显然CNF和DNF是melt的
称两个句子等价,若他们用到的变量相同,且对于任何一种变量分配,他们表达式求值后同真同假。
任何句子可以转化成等价的melt形式,只需要将多余的连词合并即可。

命题:
如果句子x1,...,xn都是melt的,则melt(Op,x1,...,xn)=[Op]+trim(x1,Op)+...+trim(xn,Op)也是melt的,其中Op是AND或OR
其中,trim(x,Op)=x[1:]若Op是x的连词,否则trim(x,Op)=[x]
证明:显然

命题:
如果句子x1,...,xn都是CNF的,
那么melt(AND,x1,...,xn)也是CNF的
证明:显然

命题:
如果句子x1,...,xn都是CNF的,将他们的子句(n个列向量)的卡氏积写成一串n阶向量y1,...,ym,|yk|=n,
那么AND,OR(y1),...,OR(ym)也是CNF的
证明:反复应用交换律即可。

基于上述命题,可知昨天的递归算法正确。

按照定义验证也很直接,只需代入全部可能的变量分配即可。python可以帮你方便地对表达式求值。
import itertools

def is_leaf(x):
    return type(x)==type("")

def to_str(x):
    return is_leaf(x) and x or "(" + (" "+x[0]+" ").join(map(lambda e : to_str(e), x[1:])) + ")"

def _nfy(x, O1, O2):
    if is_leaf(x):   return x
    if x[0]==O1:  return melt([O1] + map(lambda e : _nfy(e, O1, O2), x[1:]))
    return [O1] + [[O2]+list(c) for c in itertools.product(*map(lambda e : is_leaf(e) and [e] or _nfy(e, O1, O2)[1:], x[1:]))]       

def cnfy(x): return _nfy(x, "and", "or")

def dnfy(x): return _nfy(x, "or", "and")

def melt(x):
    if is_leaf(x): return x
    ret = x[:1]
    for e in x[1:]:
        if not is_leaf(e) and e[0]==ret[0]:   ret+=melt(e)[1:]
        else:                                 ret.append(melt(e))
    return ret

def terms(x):
    ret = set()
    q=[x]
    while len(q)>0:
        head = q.pop(0)
        if is_leaf(head):
            ret.add(head)
        else:
            q+=head[1:]
    return ret

def eq(x, y):
    xterms = terms(x)
    yterms = terms(y)
    if xterms != yterms:
        return False
    xstr=to_str(x)
    ystr=to_str(y)
    termlist = list(xterms)
    for assignment in itertools.product(*([("True", "False")]*len(termlist))):
        x1=xstr
        y1=ystr
        for i, a in enumerate(assignment):
            x1=x1.replace(termlist[i], a)
            y1=y1.replace(termlist[i], a)
        if eval(x1) != eval(y1):
            return False
    return True
       
def main():
    x = ["or", "x1", ["and", ["or", "x2", ["and", "x3", "x4", ["or", ["and", "x5", "x6"], "x7", ["or", ["and", "x8", "x9"], "x10"]]]], "x11"]]
    print terms(x)
    print eq(x, x),
    print eq(x, cnfy(x)),
    print eq(x, dnfy(x))
   
main()

P.S.
修正了_nfy函数最后一行:用不着melt

2010年12月16日星期四

convert to CNF/DNF from any simple sentence

等着看《让子弹飞》,激动得写了个小程序,把任何不含“非”的简单句子转换成CNF或DNF范式
其实只做了一件事:把AND和OR按照交换律不断展开
import itertools

def is_leaf(x):
    return type(x)==type("")

def to_str(x):
    return is_leaf(x) and x or "(" + (" "+x[0]+" ").join(map(lambda e : to_str(e), x[1:])) + ")"

def _nfy(x, O1, O2):
    if is_leaf(x):   return x
    if x[0]==O1:  return melt([O1] + map(lambda e : _nfy(e, O1, O2), x[1:]))
    return melt([O1] + [[O2]+list(c) for c in itertools.product(*map(lambda e : is_leaf(e) and [e] or _nfy(e, O1, O2)[1:], x[1:]))])
       
def cnfy(x): return _nfy(x, "AND", "OR")

def dnfy(x): return _nfy(x, "OR", "AND")

def melt(x):
    if is_leaf(x): return x
    ret = x[:1]
    for e in x[1:]:
        if not is_leaf(e) and e[0]==ret[0]:   ret+=melt(e)[1:]
        else:                                 ret.append(melt(e))
    return ret

def main():
    x = ["OR", "a", ["AND", ["OR", "b", ["AND", "c", "d", ["OR", ["AND", "e", "f"], "g", ["OR", ["AND", "h", "i"], "j"]]]], "k"]]
    print to_str(x)
    print to_str(cnfy(x))
    print to_str(dnfy(x))
   
main()

运行结果
(a OR ((b OR (c AND d AND ((e AND f) OR g OR ((h AND i) OR j)))) AND k))
((a OR b OR c) AND (a OR b OR d) AND (a OR b OR e OR g OR h OR j) AND (a OR b OR e OR g OR i OR j) AND (a OR b OR f OR g OR h OR j) AND (a OR b OR f OR g OR i OR j) AND (a OR k))
(a OR (b AND k) OR (c AND d AND e AND f AND k) OR (c AND d AND g AND k) OR (c AND d AND h AND i AND k) OR (c AND d AND j AND k))

参考
First-Order Logic Resolution Theorem Prover In Haskell

更有趣的问题是形式化验证算法的正确性,下回分解

2010年12月14日星期二

序贯蒙特卡洛(三)

命题:M2和M3都是无偏估计,E(w[1]...w[n])=Zn,其中w[k]是第k步产生的调整系数,Zn是满足约束的路径数
证明:Zn=sum{1(x) , x是所有哈密尔顿路径}其中1(x)=1当且仅当x满足约束。
1(x)=w(x)p(x)
p(x)为按算法规则产生这个路径x的概率,链式分解:p(x)=p(x[0],x[1],...,x[n])=p(x[0])p(x[1]|x[0])...p(x[n]|x[n-1])
w(x)=0当x不满足约束,否则w(x)=w[1]...w[n],w[k]=1/p(x[k]|x[k-1])
从而得到右端就是期望E(w[1]...w[n])
点评:神奇之处在于概率的链式分解给出了权重的连乘表达,同时给出了迭代构造试探路径的方法。好的构造,即方差小的构造,会以较大概率得到满足约束的路径(接受的样本),这也是Important sampling的初衷。

2010年12月12日星期日

序贯蒙特卡洛(二)代码

在实验过程中发现,M2并没有比M3、M4差太多。
M3和M4等价是因为转置不改变矩阵的积和式。物理上可以理解:每一步(一个时间)找一个地方走,和每个地方找一步(一个时间)走,其实都是一种置换。
M1虽然计算复杂度低,但是方差实在太大了,没法要。

import math
import random
import itertools
import pylab
import numpy

def weighted_choice(weights):
    cum = numpy.cumsum(weights)
    return numpy.searchsorted(cum, numpy.random.rand() * cum[-1]), cum[-1]

def restrict(A,p):
    return all([A[i][p[i]] for i in xrange(len(A))])

def make_zero_one_matrix_with_given_fillrate(dimension, fillrate):
    return numpy.random.random((dimension, dimension)) < fillrate

def perm(A):
    return sum([restrict(A, p) for p in itertools.permutations(range(len(A)))])

def eperm1(A, m):
    c=0
    p = range(len(A))
    for _ in xrange(m):
        random.shuffle(p)
        c += restrict(A, p)
    return math.factorial(len(A))*float(c)/m

def eperm2(A, m):
    Backup=A.copy()
    w=numpy.ones(m)
    for trip in xrange(m):
        A=Backup.copy()
        for step in xrange(len(A)):
            candidates = filter(lambda k : A[step][k], xrange(len(A)))
            if len(candidates) == 0:
                w[trip] = 0
                break
            k = random.choice(candidates)
            A[:,k]=False
            w[trip] *= len(candidates)
    return sum(w)/m

def eperm3(A, m):
    Backup=A.copy()
    w=numpy.ones(m)
    for trip in xrange(m):
        A=Backup.copy()
        csum=A.sum(0)
        for step in xrange(len(A)):
            for k, s in enumerate(csum):
                if s == 1:
                    A[:,k]=False
                    break #you have no other choice but goto k, w is kept
            else:
                candidates = filter(lambda k : A[step][k], xrange(len(A)))
                if len(candidates) == 0:
                    w[trip] = 0
                    break
                p = 1.0 / (csum[candidates] - 1)
                k, s = weighted_choice(p)
                w[trip] *= s/p[k]
                A[:,candidates[k]]=False
    return sum(w)/m

def eperm4(A, m):
    Backup=A.copy()
    w=numpy.ones(m)
    for trip in xrange(m):
        A=Backup.copy()
        rsum=A.sum(1)
        for pos in xrange(len(A)):
            for step, s in enumerate(rsum):
                if s == 1:
                    A[step,:]=False
                    break
            else:
                candidates = filter(lambda i : A[i][pos], xrange(len(A)))
                if len(candidates) == 0:
                    w[trip] = 0
                    break
                p = 1.0 / (rsum[candidates] - 1)
                k, s = weighted_choice(p)
                w[trip] *= s/p[k]
                A[candidates[k], :]=False
    return sum(w)/m
   

def main():
    n=15
    K=10
    m=100
    A=make_zero_one_matrix_with_given_fillrate(n,0.7)
    print(A)
    print(n,math.factorial(n),m)
#    print(perm(A))
    x=numpy.zeros((K,4))
    for i in range(K):
        x[i][0]=eperm1(A.copy(), m)
        x[i][1]=eperm2(A.copy(), m)
        x[i][2]=eperm3(A.copy(), m)
        x[i][3]=eperm4(A.copy(), m)
    print(pylab.average(x,0))
    print(pylab.var(x,0))
#    pylab.figure()
#    pylab.boxplot(x)
#    pylab.show()
   
main()

2010年12月10日星期五

序贯蒙特卡洛(二)



积和式近似

前面两个例子都没有展开,这次写具体些。
白老师曾说,计算数学其实分两块,一块研究“计数”,就是“数数儿”,另一块研究算数,就是解方程。积和式属于前者的范畴。

积和式定义
http://en.wikipedia.org/wiki/Permanent

A是一个n阶0-1方阵,积和式perm(A)=#{A(p)=1, p是(1,...,n)的一个排列},其中A(p)=A[1][p[1]]*...*A[n][p[n]]
物理意义:p是一条长度为n的哈密尔顿路径,A是约束矩阵,A[i][j]=1表示允许路径的第i步走到j。A(p)=1表示哈密尔顿路径p满足约束A。
perm(A)就是满足约束A的哈密尔顿路径的个数。
(上个版本说和哈密尔顿路径问题等价是错误的)

下面考虑计算perm(A)的蒙特卡洛方法。

M1:大海捞针
生成N个(1,...,n)的排列p_1,...,p_N
perm(A)的估计=n!*(A(p_1)+...+A(p_N))/N

M2:走一步说一步
一步步构造哈密尔顿路径,每次等概率地走到还没走过的且能走的地方,直到无路可走。
尝试m次:
w[0]=1
走第k步时,如果没路了,则w=0,尝试失败
否则w[k]=w[k-1]*没走过且可走的地方的个数
最后w[n]作为这次尝试的权重w
perm(A)的估计:m次权重w的平均值

M3:倾向于先走不容易走到的地方
思路是如果一个地方少有机会能走,则应该先走。
一步步构造哈密尔顿路径,每次走到还没走过的且能走的地方k的概率,反比于k的剩余的可走到的机会数,直到无路可走。如果一个地方必须这次去,则概率1地走向哪里。
尝试m次:
w[0]=1
走第k步时,如果没路了,则w=0,尝试失败
否则w[k]=w[k-1]*第k步选择所选落脚点的概率的倒数
最后w[n]作为这次尝试的权重w
perm(A)的估计:m次权重w的平均值

实验
A=
[[ 0.  1.  1.  1.  1.  1.  1.]
 [ 1.  1.  0.  0.  1.  1.  1.]
 [ 1.  1.  1.  1.  1.  1.  1.]
 [ 1.  1.  1.  1.  0.  1.  1.]
 [ 1.  0.  1.  1.  1.  0.  0.]
 [ 1.  1.  1.  1.  1.  1.  1.]
 [ 1.  1.  0.  1.  1.  1.  1.]]
真值perm(A)=1434
三种方法估计的均值(样本数20)
[ 1415.          1448.91071429  1434.89212938]
下图是M1,M2,M3的Boxplot
可以明显看出方差的减小


2010年12月1日星期三

n次根号下p[1]..p[s]是无理数

命题:
n次根号下p[1]..p[s]是无理数,其中p[i]是素数,n>1
证明:
假设存在自然数a,b使得(p[1]..p[s])^(1/n)=a/b
则对a,b做素因子分解
a=p[1]^a[1]..p[s]^a[s]
b=p[1]^b[1]..p[s]^b[s]
其中a[i],b[i]为非负整数
则p[i]^(1/n)=p[i]^(a[i]-b[i])
即1=n(a[i]-b[i]),这对于n>1是不可能的

2010年11月10日星期三

如果饭否回来了

如果饭否回来了,我将表示情绪稳定。

2010年11月6日星期六

序贯蒙特卡洛(一)




遗传学的例子

大三那年自从上了白老师的一个讨论班之后,就对生物学中的很多问题很感兴趣。这 本教材中也有不少生物方面的例子,让我时不时眼前一亮。

考虑一组等位基因的演化历史。我们能观察到的,只有今天这组等位基因在某种群的分 布。我们提出一个单倍体演化模型,只考虑繁殖和突变,一步概率转移矩阵由某个待定系数t唯一决定。用极大似然方法估计这个参数,也就是极大化历史演化成今 天这个样子的概率。

按照定义,要计算

P(H0)=sum{P(H)|H 的最后一代是H0}

其中H0是 今天的分布,H是历史演化过程。这在计算上是不可行的,因为H的空间太大了。自然想到蒙特卡洛。

任意给定参数t,根据贝叶斯公式,可以知 道P(H-1|H0),从而每一步,根据今天的分布,按条件概率采样尝试回到昨天的分布,也就是P(H-1|H0)的估计。反复为之,可以得到 P(H|H0)的估计,也就是在今天这样分布的条件下,某历史演化过程存在的概率。
从而得到P(H0)的估计:P(H)/P(H|H0)的多次实 现的均值。




聚合物的例子

考虑一个随机方案来构造一个长度为N的、在平面网格点上的、不能自交的 路径,使得任何两条满足条件的路径出现的概率相同。不失一般性,假设从原点出发第一步向右走。
方案A:
每一步从{左转、右转、直行}三个 中等概率选择一个,如果自交了,则拒绝当前路径从头再来。否则走到长度为N接受它。
方案B:
每一步如果没有能走的地方,则拒绝当前路径从 头再来,否则从能走的1或2或3个地方中等概率选择一个,所谓能走的地方,指的是{左转、右转、直行}的三个点中之前没有到达过的那些。

B 看起来更有效率,然而不幸的是它不是无偏的:它更倾向于走出一条压紧、褶皱的路径。例如,对于N=4,走出[右,右,右,右]的概率就要小于[右,上, 左,左]
为了使得B均匀分布,需要为它生成的每条路经赋一个修正权重,而这个修正权重的计算,正是蒙特卡洛方法的精髓。

2010年11月2日星期二

“给我照张相吧,挤死了”

其实我没想照他。在按下快门的刹那,他喊了这么一句。令我印象深刻。
可惜我不是记者,我只能发到自己的技术博客上。
嗯。
"这哪儿算挤",人们说。

2010年10月25日星期一

FreeWheel isn't sitting in the middle, it is making money

上周跟着FreeWheel的男女老少去清华做校园招聘,CEO想方设法让同学们理解FreeWheel是怎么赚钱的,可是同学们依旧一脸茫然,如同他们走进教室时一样,如同两年前的坐在这里的那些人一样。
PPT上有些图,如同所有公司宣传一样,FW被画在中间,上下左右是各种类型的客户,仿佛天地因FW的存在而运转,仿佛钱从四面八方向FW而来。这种八卦图看多了,同学不禁会想,谁都说自己是中心,到底谁是中心呢?
在我看来,其实大家是在一个N维单位球面上。一共有N个公司。每个公司都是这个球面上的一个点,因此任何两个点之间的距离都相同。但是,由于每个公司都觉得世界是平的,而且看到所有其他公司到自己的距离都相同,所以觉得自己必然巍然中央,所有其他公司落在以自己为中心,半径为1的圆环上。

同学们也觉得世界是平的,而且点特别多,所以觉得应该至多存在一个点,使得所有其他点到这个点的距离相同。这个点就是中心。

争谁最大最牛没意义,争谁是宇宙中心也没意义,挣钱才是真的。

2010年10月23日星期六

没有光驱?死不了的

过去我一直觉得重装系统就必须有光驱,而且得有张光盘带光启,非常神奇。这周买了个小上网本,没有光驱,没有预装操作系统。表示惊讶之后,只好在另外一台能上网的电脑上反复百度,终于用U盘搞定了。参考:http://apps.hi.baidu.com/share/detail/866549

2010年10月14日星期四

基本的蒙特卡洛方法(三)​

基本的蒙特卡洛方法(三)​

根据蒙特卡洛方法得到的样本总是无偏的,即样本期望等于总体期望。实践中,期望相同还不够,还要求方差尽量小。因此方差减小技术是蒙特卡洛方法极重要的一部分。

分层采样(Stratified sampling)

充分利用积分的可加性。

仍然考虑计算广告学中的例子,考察样本空间E=<V, C>,其中V是视频集合,C是城市集合。一个点e=(v, c)表示一个来自视频v且来自城市c的广告请求。

给定E的子集E1,求E1中的广告请求个数。

如何广告请求顺序写在了一个文件里,那我们只好随机采。

如果来自不同城市的请求分别写在不同的文件里,那我们可以每个城市随机采,从而有望减小方差。

方差减小的充要条件可以模糊地形容为:同一个城市的用户看视频的偏好相近。

不证明了,举个小规模的例子:

两个城市,100个视频,每个城市每个视频恰好有一个广告请求,共200个。考虑E1={城市1的视频1到90,城市2的视频1到10}。我们一眼就看出来E1中广告请求有90+10=100个。让我们采两个样本,来估计之。

在全空间上随机采:从200中等概率、有重复地采2个,X个属于E1,则100X是无偏估计。经计算,X的方差是0.5。

分层采:从每个城市各采1个,城市1中采中属于E1的个数是X1,城市2中采中属于E1的个数是X2,则100X=100(X1+X2)是无偏估计。经计算,X的方差是0.18。

控制变量法

也叫借鸡下蛋法,呵呵。

这部分的分析很简洁漂亮,情不自禁地说两句。

假设我们关心随机变量X的期望EX。已知另一个随机变量C的期望EC,而且X和C比较相关。用待定系数法,令随机变量X(b)=X-b(C-EC)。(原书中写成了+b有误!)

显然EX(b)=EX,因此,无论系数b如何取,求EX(b)即可。关键在于,适当选取b使得X(b)的方差尽量小。

varX(b)=varX-2bcov(X,C)+b2varC

中学生都会,这是关于b的抛物线,当b=cov(X,C)/varC时最小,最小值为


神了,这看上去就是白给。其实有很多代价:首先要知道EC,其次还要能算b,这太浪漫了。

举个浪漫的例子吧。考虑100个视频请求,一个广告C要投放在1~60,另一个广告X要投放在1~50。假设C的情况我们了解了,X是新来的,要通过采样估计它的访问量。

先按公式算出来b=5/6。从100个请求中随机采1个,属于广告X投放范围的有X个,属于广告C投放范围的有C个,则X(b)=X-5/6*C+0.5。

X

C

X(b)

概率

0

0

0.5

0.4

0

1

-1/3

0.1

1

0

1.5

0

1
1
2/3
0.5


varX(b)=1/3*0.4+1/9*0.1+1/2*4/9=1/12

varX=1/4






2010年10月7日星期四

协稳的假设有多强?

多元时间序列分析一个特别醒目的概念,就是协稳:若干个时间序列,他们可能每个都不平稳,但是他们的某个线性组合是平稳的。
在计算广告学实践中,我们考虑一个网站的若干视频,每个视频的日访问量序列看作一个时间序列。很少有哪个视频的访问量是平稳序列,因为人们总是喜新厌旧,一个视频刚上线的一段时间内,看的人越来越多,一段时间后,就渐渐地没什么人看了。但是在成熟、稳定的市场(而非快速增长的市场),一个网站的视频日访问总量在较长一段时间内是平稳的,因为随着老视频陆续淡出,不断有新的视频上线,形成协稳的格局。
当然这只是美好的假设。现在让我们来看看这个假设究竟有多强。
假设有N个视频,用N维时间序列X(t)表示。M个网站,用M维时间序列b(t)表示。假设视频和网站的多对多关系用M*N的0-1矩阵A表示:即A[i][j]=1当且仅当视频j出现在网站i上。于是有下列线性方程组:A*X(t)=b(t)
模型假设说每个网站的流量都平稳,也就是A*EX(t)=Eb(t)=b。考虑方程组AX=b的非负数解,下面简称解。
当方程组无解时,说明"每个网站的流量都平稳"的假设不成立。
当方程组解存在唯一时,X(t)也平稳,这是退化的情形,即所有视频的流量都平稳,这在现实世界很难发生。
当方程组解存在且不唯一时,才是模型假设的经典情形,即若干非平稳序列的线性组合平稳,由线性代数知识可知,此时A的行秩小于N是必要的。

2010年10月3日星期日

修改App Engine上的cronjob

以前用Windows的时候用Java+Eclipse写了一个免费收短信的工具,每天两条,一个天气预报,一个每日英语。
最近觉得每天发的那个每日英语内容没什么意思,决定取消掉。可是已经换了电脑,Eclipse一键部署的功能没有了。
全新的Ubuntu下,下载了一个Python的App Engine SDK,唯一的目的就是停掉发每日英语的cronjob。
原本担心:
1. 没有代码树如何提交?
2. 当时部署的是Java的环境,现在如何用Python工具修改?
实践发现:
只需要建一个目录a,下面放两个文件:
空白的应用程序配置文件
app.yaml
application: freewheel-sms
version: 2
runtime: python
api_version: 1

handlers:
- url: /
  static_dir: /

和新的cronjob配置文件
cron.yaml
cron:
- description: daily_weather
  url: /routine
  schedule: every day 21:02
  timezone: Asia/Shanghai

然后执行
python appcfg.py update_cron a

即可
这个部署脚本特别单纯,不管代码diff,也不管原来部署的运行时语言,只管更新cronjob

2010年10月2日星期六

How to use VPN in Ubuntu 10.4

参考
http://www.mcmaster.ca/uts/network/vpn/vpnclient_linux.html
http://www.lamnk.com/blog/vpn/how-to-install-cisco-vpn-client-on-ubuntu-jaunty-jackalope-and-karmic-koala-64-bit/
步骤
第一步:下载并解压
http://www.mcmaster.ca//uts/network/software/vpnclient-linux-x86_64-4.8.02.0030-k9.tar.gz
第二步:打补丁de关键三句话

wget http://lamnk.com/download/vpnclient-linux-2.6.31-final.diff

patch < ./vpnclient-linux-2.6.31-final.diff

sudo sed -i 's/const\ struct\ net_device_ops\ \*netdev_ops;/struct\ net_device_ops\ \*netdev_ops;/' `find /usr/src -name netdevice.h`

第三步:安装
sudo ./vpn_install
sudo /etc/init.d/vpnclient_init start


第四步:配置并运行
把你公司的ABCD.pcf复制到
/etc/opt/cisco-vpnclient/Profiles
然后就可以运行了
vpnclient connect ABCD

我的测试环境:
Linux version 2.6.32-21-generic (buildd@rothera) (gcc version 4.4.3 (Ubuntu 4.4.3-4ubuntu5) ) #32-Ubuntu SMP Fri Apr 16 08:10:02 UTC 2010


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,显示调试信息,例如告诉你为什么某一组参数无法渲染,非常好用

2010年8月18日星期三

不能在偏序集上排序


bool myfunction (int i,int j) { return (i==4) && (j ==2); }



int main() {

   {

       {

       std::vector<int> x(4);

       for(int i=0;i<4;++i)x[i] = (i + 1);

       for(int i=0;i<4;++i)printf("%d,", x[i]);

       printf("\n");

       std::stable_sort(x.begin(), x.end(), myfunction);

       for(int i=0;i<4;++i)printf("%d,", x[i]);

       printf("\n");

       }

       {

       std::vector<int> x(4);

       for(int i=0;i<4;++i)x[i] = (i + 1);

       x[0]=2;

       x[1]=4;

       for(int i=0;i<4;++i)printf("%d,", x[i]);

       printf("\n");

       std::stable_sort(x.begin(), x.end(), myfunction);

       for(int i=0;i<4;++i)printf("%d,", x[i]);

       printf("\n");

       }

       return 0;

   }



Result



1,2,3,4,

1,2,3,4,

2,4,3,4,

4,2,3,4,





http://middleangle.com/rif/derifatives/Home/13/sorts-and-partial-orderings

[12:32:59] wangyuantao: 构造了一个反例,不能排序。序列1,2,3,4,只定义4<2,按照归并排序,结果还是1,2,3,4

[12:33:07] 翟倩(Zhai Qian): 这是谁的blog啊

[12:33:18] wangyuantao: 不认识google出来的

[12:47:14] wangyuantao: 关于偏序集上的排序和TopN算法,这篇文章被引用的比较多

http://www.siam.org/proceedings/soda/2009/SODA09_044_daskalakisc.pdf

光看了摘要,貌似最好的结果也不保证最坏情况下的复杂度比O(n^2)好。所以,就按照目前的实现两两比吧

发送自 HTC



2010年7月27日星期二

不等式移项

        unsigned int x=6;

        int y=8;

        printf("(x>y)=%d\n", x>y);

        printf("(x-y>0)=%d\n", x-y>0);

        printf("(y-x<0)=%d\n", y-x<0);

今天实验后才发现,这三个表达式两两不等价,C语言的类型转换真是博大精深啊。
第一个:编译警告,有符号和无符号整数比大小有危险。
第二个:错的一塌糊涂。x-y返回一个无符号整数-2,竟然大于0了。
第三个:对。

比较稳妥的方法是显式类型转换。

2010年5月8日星期六

挑鸡蛋

下午去菜市场买鸡蛋,挑大的买,忽然脑袋里蹦出来这样一个问题:
假设鸡蛋的大小服从[a,b]均匀分布,那么如果我任意取一个,则期望大小为(a+b)/2。如果我任取两个,要较大的,期望是多少呢?感觉应该比(a+b)/2大。
回来用期望的定义算了一下,是a/3+2b/3。直观上这个结论也非常好接受:具有某种对称性。推广到N个也成立。对离散情形也有类似的结论。

2010年3月14日星期日

GSL + MinGW + Eclipse

想当年在学校的时候只用过Matlab,没有用C++做过科学计算,现在觉得很遗憾。周末在家琢磨GSL,懒得切到Linux,又不敢往公司的开发环境上乱装东西,所以很没追求的又用最熟悉的MinGW+Eclipse了。

参考
http://www.diybl.com/course/3_program/c++/cppsl/2008326/107387.html

步骤
1.下载并安装
http://gnuwin32.sourceforge.net/downlinks/gsl.php
2. 安装库
复制C:\Program Files\GnuWin32\bin下的libgsl.dll和libgslcblas.dll到C:\MinGW\bin
复制C:\Program Files\GnuWin32\lib下的libgsl.a 和 libgslcblas.a到C:\MinGW\lib
3. 在C++ Project中添加include path
"C:\Program Files\GnuWin32\include"
4. 运行实例
#include <gsl/gsl_sf_bessel.h>
int main()
    {
           double x = 5.0;
           double y = gsl_sf_bessel_J0 (x);
           printf ("J0(%g) = %.18e\n", x, y);

    }

2010年3月6日星期六

[Reading] Optimal Online Assignment with Forecasts

正在读Yahoo! Research的一篇旧文。文章中,将在线广告投放的问题,理解成一种assignment problem。一个顶点集是advertiser,另一个顶点集是user。
模型中有两点关键:
1、online,user集中的点依次到达,要求依次assign给advertiser集中的点。
2、forecast,允许对offline数据进行采样,生成"compact allocation plan",用于指导online assignment。
约束条件:
每个user不能提供太多的流量(User Experience?)
每个advertiser的需求有不等式约束,可能是流量下限(至少达到一定的投放效果),也可能是流量上限(不能超预算)。
目标函数:
生没看懂。。。
文章的数学味很浓,我已经把凸优化、对偶空间那套理论忘光了,读起来相当吃力。。。嗯,正好督促我复习一下

2010年2月28日星期日

Space rail

听zeta同学推荐,订了一个space rail,貌似很好玩的样子,相当期待。在YY的时候,忽然发现这个东东貌似很适合搞软件的团队做team building~
需求:造出来的过山车要能转起来(functional),而且样子要酷(polished)
设计:利用力学原理(theory),设计轨道(architecture)造型(product)
构建:每个人组装(implement)一个模块(module),确保小球在给定进入速度范围内可以通过自己的模块,并且离开速度满足给定要求(unit test),然后集成(integration)
管理:分配底座、支架、轨道等零件(resource)给各个玩家(engineer),协调各个模块组装的进度等
迭代:当过山车做装好后(release),每当有了新装置或设想(new feature),稍加改动(refactoring)就可以安装到成品中,而不需要全部拆散重来(redesign/rewrite),并且递归地,改动后的过山车仍然具有此属性(continual improvement)。

2010年2月19日星期五

春节。上海。西塘。杭州。

以为,杭州是一个充满了怨念的城市,可我还是去了。
大年初二凌晨5点,和LP踏上了南灰的旅途。落地是上海时间10:30,淅淅沥沥下着雨,冷得棉袄里面好像装着冰。下午从南站往嘉善。
晚上赶到西塘,被嘉善的出租车黑了。不能怪公交车收的早,也不能怪正品黑车杀人不见血,只能怪晚辈我人生地不熟不敢跟人家犯蛮,成全了那披着出租车绿皮的黑车。
晚上住在杨老板家,杨老板人很好,给我们热粽子吃。终于见识了南方的冬天没有暖气的冷是啥滋味。
初三,天晴了,西塘上午的阳光真的非常非常亚克西。只想坐在水边的石台上一直发呆,看着游人走着石板路,看着爱人吃着煮青豆,看着老太太不急不忙地递过来臭豆腐、粉蒸肉,惬意地只想睡了算了。LP食欲和玩心超好,逢食必吃,逢店必进,逢景必拍照,活蹦乱跳。
下午在杨老板的帮助下换到了提前发往杭州的汽车票,于是两个小P孩拍PP走了。一觉醒来,就到了莫干山路pod,传说中的家庭房,超cute的小房间。晚饭杭大路尘衣食府,尘依小炒和扒羊肉相当不错(偷偷拍照),片儿川没有下次了。
初四,彻底在武林门/广场迷失了方向,因为杭州修轨道一号线,公交车已经不再既定路线上行动了。灵隐,九溪,都没去成,因为找不到Y1,Y2,很崩溃。只想摔了手机,卸了Google
Map。中午灰溜溜地去河坊街吃王润兴,吃到一半犯病了,不知道是不是因为太难吃了。比08年来杭州那次犯得还厉害。下午继续在河坊街逛,买了丝巾和被面,反正不会砍价,索性没有砍,宰就宰吧,反正杭州我已经来四次了,没有下次了。
初四晚上,一咬牙一瞪眼,回到北京了。

春节前:
阳光明媚。上班浪费。

春节:
怨念沪杭。怀念西塘。

春节后:
龙体欠安。提前下班。

2010年2月12日星期五

Happy Chinese New Year!

富饶之城,赛乐堡,全城热恋,放鸽子

春节快到了,娱乐活动多了起来。
周四午饭和501的同事一起吃,那些家伙比传说中的客气。饭后响应weiwei的号召去玩富饶之城。我是第一次玩,还是挺有意思的。桌游的魅力在于现场互动,大家都相当生动。Free有一局不幸地当了我的替罪羊被刺了,然后从此一蹶不振,Snow是永远的军阀,挣钱就是拆台,口号是我赢不了谁也别想赢,Dandan为了争第二费劲了心思,我一直很想在博物馆下面藏一个大学,因为觉得那样特别逗,不过一上来就被军阀摧了。赢家当然是weiwei,在正确的时机建造了八级浮屠,不愧是有经验的玩家。
下午公司里人越来越少,每个人起身离去都要和剩下的所有人告别说"明年见",搞得我越来越没心情coding了。下班后跑到赛乐堡和朋友K歌,含自助晚餐。因为下班前吃了一包饼干,不饿,于是我说,"你们先吃,我给你们唱歌!",然后开开原唱,到处找调。找到所有人都吃饱了我也没找到调。5个人水平都一般,个别人个别曲目个别段落有感觉有旋律,大多时候有感觉没旋律或没感觉。然而非常High,不知不觉4个小时过去了,还有4屏的曲目没轮上,就很伤心的走了。
周五叫嚣着去上班,然而早上起床后发现阳光明媚,心想上班浪费,于是和MM去看全城热恋了。虽然很肥皂很商业,但是还是有几处很有趣的。最逗的是字幕,在显示人物发短信的时候,居然可以支持退格!还有就是显示制作群的时候,文字可以配合音乐抖啊抖的。两只小鱼的动画也很搞。
下午回家干活,处理了服务器上的一个小事件,又读了篇PM的Design。然后上Skype,发现DZF在线上,一个昨天也号称今天要去上班的童鞋,结果一问,啊哈,也没去公司,再一问XXX,也没去。。。瀑布汗,今天503貌似就剩前台姐姐了。。。心中顿时升起一坨坨的歉意。
明儿就三十儿了,大后天就去南方玩了,保佑我们的Server keep alive吧!

2010年2月10日星期三

The limitation of the sampling algorithm: a simple case

一个不能用采样算法的问题的例子
网站某天有N个request,每个request都包含一个session_id,现在希望统计y=平均每个session有几次request。若M为distinct的session的个数(通常M<N),则y=N/M。
现在从N个request中随机取n个,记这其中distinct的session的个数为m,令随机变量X=n/m,则EX!=y
反例:当n=1时,m恒=1,X恒=1。
问题:是否存在一个采样算法,使得算法的输出是y的无偏、相合估计?

2010年2月7日星期日

The impact of resampling to RSE in sampling algorithms

采样算法中,是否可重复采样对相对标准差的影响

在古典概型中有放回和无放回差别很大。在随机采样算法中,对已经采过的样本,是否还可以再采?
问题:
假设有N个零件,次品率为p(即有pN个次品)。现在已知N,要求通过采样率为q(检查qN次)算法来估计p。

有放回算法:
依次生成qN个独立同分布U[0,N)的随机数xi,检查第xi个零件是否为次品,记随机变量X为检查出的次品的个数。来求X的分布。
令Xi=1表示第xi个零件为次品,则X=X1+...+X_qN
P(Xi=1)=p,EXi=p,根据期望线性性,EX=EX1+EX_qN=pqN,DX=pqN(1-p)
相对标准差的平方RSE2为DX/(EX)^2=(1-p)/(pqN)

无放回算法:
一次从N个零件中随机抽取qN个零件,记其中次品个数为X,则X服从超几何分布
EX=pqN,DX=p(1-p)q(1-q)N,RSE2=(1-p)(1-q)/(pqN)

对比发现,两个算法的RSE2就差个(1-q),即相同采样率q下,无放回更精确。当然,对于千分之一的采样率,二者差别可以忽略,但是如果要采50%的样,最好用无放回算法了。这也符合物理直观。
但是不管有没有放回,RSE都关于p,q,N单调递减,因此,固定p(p来自问题域),可以通过调节qN来寻求性能与精度的平衡――
当N较小时,用比较大的q来维持精度,反正样本数qN也不大。
当N较大时,不需要很大的q,因为反正RSE2的分母上已经有个N了:)

2010年1月26日星期二

Share Two Google Calendar App

Calendars
1) Daily English
ID
gfn0kpsechk859laf5effn31uo@group.calendar.google.com
preview
http://www.google.com/calendar/embed?src=gfn0kpsechk859laf5effn31uo%40group.calendar.google.com&ctz=Asia/Shanghai


2) Beijing Weather
ID
oiio53f5pflo33a3nd328gfics@group.calendar.google.com
preview
http://www.google.com/calendar/embed?src=oiio53f5pflo33a3nd328gfics%40group.calendar.google.com&ctz=Asia/Shanghai

How to use
step 1
go to your Google Calendar and paste the ID above into the text box title with "Add a friend's calendar" at left;
step 2 (optional)
Add SMS as a Notification for the calendar, then you can receive Daily English/Beijing Weather everyday!

Copyright
Daily English piped from dict.cn
Beijing Weather fetched from weather.com.cn

More
Tell me what you want, I'm pleasure to pipe it into your SMS :)

2010年1月24日星期日

Ubuntu 9.10 remember last boot choice

http://ubuntuforums.org/showthread.php?t=1195275

step 1
sudo gedit "/etc/default/grub"
change
"GRUB_DEFAULT=0"
to
"GRUB_DEFAULT=saved"

step 2
sudo update-grub

step 3
reboot

2010年1月17日星期日

Reading about Refactoring

I'm reading a book about refactoring these days, and found some good points useful in my work.
1. Definition
Refactor do not change any observable behavior of software
//But new features are needed, they are mixed in time, but should not mixed in mind. They are two things!
//Programmers should never add any new feature without PM design while refactoring.
2. Aim
Improving the design
//Performance
//Data Consistency
//Reuse the logic
Make the code easy to understand
//Not easy, because the logic is really complicated sometimes. We should try.
Find bugs
//Clarifying the assumptions of each unit!
Move on fast
//Decayed design is a huge burden. Without it, we can move fast!
3. When and when not
Three strikes and you refactor.
//I used to refactor at two strike, and always find it unnecessary in the end.
Right time: adding new features, fixing bugs, doing a code review
Wrong time: Need to rewrite, near to deadline
//Stop refactoring on branch since staging release
4. Limitation
DB, interface
//change them carefully. Data migration is always dangerous. Other team depend on your interface, so if you change it, they have to change!
Design
//The feature design is sometimes lead to bad structure of program. Engineers are obliged to reduce the bad design.
5. Bad smell
Long parameters list
//Erlang program is easy to be writen into these style
Switch statement instead of polymorph
//Erlang has no switch, it is good! I love polymorph
Lazy class
//Someone loves it, but I never use it
Use member fields as global variable
//I hate it very much, although I do it sometimes
//Member fields are state of an instance, rather than a basket holding the mess you do not know where to put it.
Data classes
//All the time I wish all small data struct can be thrown into protobuf, but sometimes they need quick lookup or other functions. It is difficult to trade off.
6. Test
When you get a bug report, start by writing a unit test that expose the bug.
//I'm not sure if all the refactoring can be covered by unit test. I love TDD, however, I find it hard to maintain the unit tests because the inner interfaces changed frequently. Unit tests turn from my proud into my burden that I have to change accordingly. So I think writing too much unit test is not a good idea. Today I agree that unit tests have two major purpose: 1) Code coverage, 2) reported bugs coverage

2010年1月5日星期二

2010 Resolution

Many resolutions around, even a generator was found:p

My wish list this year

Healthy
Weekly sweat
Gain weight
Tech
Expert in Linux Shell/Python/C++
Get start with Design Pattern
Sci
Keep on following Recommendation Algorithms
Get start with Microeconomics
Sandbox
Dig into App Engine/GData
Popularize my Blog