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是不可能的