王元涛的Blog

凯立德刷K币

2012-02-12T19:52:00.001+08:00

话说道高一尺魔高一丈。先从最弱的一个魔开始吧。

打开一个记事本，粘贴如下代码另存为kld.user.js，然后拖到Chrome或Firefox（Firefox需要安装GreaseMonkey插件）里。

然后在地址栏中输入：

http://www.kldjy.com/home.php#3

即可获得10个K币：因为访问了个人空间。虽说打开了12个Tab，但是家园的规则是每天由于访问个人空间获得K币的上限=10。

当然，要提前弄好自动登录。为了不影响正常使用家园空间首页，只对加了"#3"无害后缀的网址做批量弹Tab了。

// ==UserScript==

// @name KLD K coins

// @include http://www.kldjy.com/home.php#3

// @run-at document-end

// @description open 10 space

// @version 1.0

// ==/UserScript==

function main() {

for(i=0;i<$("portal_block_80_content").children[0].children[0].children.length;++i)

{

window.open($("portal_block_80_content").children[0].children[0].children[i].children[0].href, "_blank")

}

var inject = document.createElement("script");

inject.setAttribute("type", "text/javascript");

inject.appendChild(document.createTextNode("(" + main + ")()"));

document.body.appendChild(inject);

如果我的证明是正确的，就不需要任何其他形式的肯定

2012-02-11T23:50:00.001+08:00

想健康并且不浮躁地活，看来必须读一点书了。

2012年1月1日不是2012年的第一周

2012-02-03T21:56:00.001+08:00

> $ date -d "2012-01-01 00:00:00" +%U

> 01

> but in "man date"

> %U week number of year, with Sunday as first day of week (00..53)

> Is it a bug?

给gnu发了封信，得到了详尽的解答，长知识了。

ISO 8601 是一个神奇的标准，它定义了周数应该怎样数。

美国人按直觉定义每年1月1日所在的周为一年的第一周，而欧洲人制定的标准 ISO 8601 则不然。

首先， ISO 8601 定义每周的第一天为周日，其次，如果一年的第一周天数没有过半，则归入上一年的最后一周，第一个天数过半的周才真的算一年的第一周。

POSIX提供了支持这个标准的实现，也给出了多种其他符合直觉的选择，info date仔细读一下，深了去了。

查了一下中国国标，也是按这个ISO定义的。

幂级数

2012-01-18T21:37:00.001+08:00

假设有台印钞机每个月能印10元，不费纸不费电还终身保修，卖100元一台，要不要买一台呢？

假设月期望收益率是10%，那么还是不买为好，除非你忘了幂级数。

第一个月的收入折现：10/(1+0.1)=9.1，也就是说，第一个月月底的现金流入，在现在看来，价值只有9.1元。

同理，第N个月的工资折现：10/(1+0.1)^N

N趋于无穷的时候算部分和的极限=100

那应该怎么办呢？

把全部积蓄都拿来买这种印钞机，然后高价兜售给那些不懂幂级数的人XD

./

2012-01-14T14:36:00.001+08:00

有一天，我的C++程序在工作目录下建立了一个名为-1的目录，我想看看里面是什么鬼东西，不料
-bash-3.2$ ls -1/
ls: invalid option -- /
Try `ls --help' for more information.
分明是Tab出来的，怎么会这样。。。估计是ls的bug吧。cd进去看看？
-bash-3.2$ cd -1/
-bash: cd: -1: invalid option
cd: usage: cd [-L|-P] [dir]
我去。。。估计是没用了，删了吧
-bash-3.2$ rm -rf -1
rm: invalid option -- 1
Try `rm ./-1' to remove the file `-1'.
Try `rm --help' for more information.
原来如此，
rm -r ./-1
之后就天下太平了。

有一天，我的C++程序在工作目录下建立了一个名为~的目录，我一看
-bash-3.2$ tree
.
|-- bin
| `-- run
`-- ~
`-- workplace

3 directories, 1 file

没东西，那就删了吧
-bash-3.2$ rm -rf ~
rm: cannot remove directory `/home/ytwang': Permission denied
好生奇怪，我自己的程序跑出来的目录怎么会没有权限，难道程序里面用了sudo不成？。。。等等，为什么是/home。。。啊呀完了，赶快
-bash-3.2$ pwd
/home/ytwang/workplace
-bash-3.2$ ll
total 0
-bash-3.2$ cd
-bash-3.2$ ll
total 0
顿时心中拔凉拔凉的。。。

==结束语==
你永远不知道自己那充满了灵异事件的程序会创建出什么诡异的目录来，一旦他们被放出来，就将会像幽灵一样——你看不见他，而当你想砍死他的时候，自己就死掉了
要养成用./表示当前路径的好习惯
当你听说哪个同事不小心删除了某个重要文件的时候，不要惊讶，要同情

驾照换证攻略

2012-01-07T21:09:00.001+08:00

转眼驾照6年到期，要换本了。
第-1天：早上去照相馆照1寸照片，官网上给的规格说明很复杂，跟师傅说驾照用的，师傅都知道。非立等可取。15元。
第0天：早上去照相馆取相片。中午饭后在办公室打印《申请表》《身体条件证明》空白表，把自己能填的都填了，各贴一张照片，复印身份证正反面。
1点从朝阳门出发，打的到东直门医院，堵车（亏大了），11元，30分钟。
进门左手边有个门诊办公室，看一眼色盲的本本，然后大夫就边聊天便给你把表填好了，盖章，11元，2分钟。
出医院一路向南走，穿过海运仓小区，过东四十条的马路，再穿过军区总医院（好像这里也能体检），就到达了豆瓣胡同，15分钟。
在一条东西向的巨窄无比的小胡同里的路南，看到一个巨小无比的小门脸，上面写着"顺天府超市"，东城交通支队的门就在那旁边了，巨不起眼。
拿号，没什么人，30秒后就到我了。给他《申请表》《身体条件证明》，身份证及复印件，再一张一寸照片，旧驾照。
他在电脑前敲一阵，然后换个窗口给他10元，就看到一张滚烫的新驾照出炉了。5分钟。
走回朝阳门，也就2点多。

2012 resolution

2011-12-31T20:44:00.001+08:00

每年元旦前后都会制定来年计划，可是一年到头来，回头看看履行的情况总让人心灰意冷。
晚上老婆去天坛许愿，出门前问我来年有什么愿望。
我想了想"三条腿的理论"，说：长胖。

我家的挂表

2011-12-11T16:56:00.001+08:00

昨天在朝阳大悦城买的
超静音

新手提车小记

2011-11-30T15:08:00.001+08:00

用chli的话说，就这样一下回到了解放前。用工作以来的全部积蓄，买了辆小汽车。号快到期了，再加上反正也不懂，就随便买了一个。

虽然交款第二天这款车型就上了央视的每周质量报告一通批判，可我还是很开心。心想有了车以后，去看父母再也不用在公交车上浪费宝贵的时间了。

销售小李帮我把车从4S店的洗车房里开出来，放在了青年路的路边换我上去，嘱咐我说，别忘了松手刹。

离开驾校5年零9个月以来，第一次自己开车。之前1个月摸过两次朋友的车，不过都没开出去多远，一次在院子里转了转，一次在试驾路线上兜了半圈，而且都有司机在边上指挥着。

这次，说实话，很紧张。

我坐上驾驶座，等老婆坐上副驾驶座，都扣好安全带，就松手刹、给油、窜了出去。

这时销售敲车门，大声问，松手刹了吗？我踩住刹车回答，松了，然后又给油窜了出去。

刚窜出去又一脚刹车，让左边一辆残疾三轮过去。

等他过去之后，缓缓给油，慢慢开起来了。在青年路的试驾路线上开了一圈，试驾路线上几乎一辆车都没有，觉得太假了，于是回到了姚家园路，准备回家。

本来在手机上安装调试好了付费版的凯立德导航软件，下载了离线版北京地图，可是完全无心打开，只觉得周围都是车，特紧张。根本不敢分心去掏手机。

心想，嗨，我车没牌我怕谁。就顶着一张临时牌照一路。。。蜗行回家了。

按照常识最外侧车道用于右转弯或者出主路，那我就在最里侧走吧。没想到这条车道是左转弯掉头的。

汗死。

由东向西到了东三环，要右转，狠心往右并线要出主路，右侧车队死活不给空档，还挤进来一辆9XX公交挡住了我的半个屏幕。左右摇摆之际，后面喇叭声此起彼伏。

到了出口的位置，妈的老子还在直行道上打着蹦灯没拐出去呢，急了，一把转盘一脚油门挤了出去，然后一脚刹车踩死定在了9XX屁股后面。

原以为老婆在旁边吓得惊惶失措了，没想她还挺高兴的，说道："老公，你真棒！"

汗死。

到了桥下，发现其实刚才沿着主路直行到桥下也能右转，没必要拼了命地提早挤到辅路上来。

三环上走不起来。40%时间点着刹车往前蹭，30%时间滑行，5%的时间给点油，其他时间拉上手刹开风景。

当然，看风景的心思是没有的，只有在车里瑟瑟发抖的份：紧张+不知道空调热风怎么开。

安华桥掉头，两次右转270度=负180度那种，出主路又出不来了。打着蹦灯看着漫长的辅路车队，颇有沉舟侧畔千帆过的感觉。然而不是卧在船上看，而是要扭着脖子看。因为车身已经转过了30度，这次右侧后视镜视野完全不对，看不到辅路车辆，只能回头通过后排右侧车窗观察。差不多等了一个红灯的功夫，终于没车了，转上了安华桥。

终于到家了。可是。。。噩梦才刚刚开始！

最恐怖的路段，是进小区门之后。这是一个老小区，路只有三道宽，一道是规划停车位，停满了大大小小的车，另一侧地上画着禁止停车的标志，也停满了大大小小的车。

进门保安撕条递给我，我给油冲隔离墩，没控制好，留给他一个屁股。放了3个车窗才把自己的车窗放下来，伸出头去说，我停这儿去办车证，行吗。他说，这不能停，往前去，停前头那儿。

我再想往前开，发现已经不能开了，已经快顶到前面停那儿那车上了。想那就往后倒一点吧，后面进来一辆金杯，滴滴冲我按喇叭。

我勒个去。

我就在那里踩着刹车原地打轮，不敢碰油门。

金杯按耐不住，自己往后倒了一点，噌地窜了进去。

我一着急，刹车踩成了油门，往后跳出去老远。好在后面没人，赶快一脚刹车站住了。

看来倒车的时候不仅容易分不清左右，而且也容易分不清前后啊。

稳定一下情绪，踩着脚刹靠怠速驶入小区，边爬边寻觅车位。

这时对面从饭馆开出来一辆运货的小面包，迎面驶来。怎么办，只有一条可用车道啊。千军万马过独木桥么。狭路相逢勇者胜么。桥到船头自然直么。车到山前可路在哪儿呢。

缓缓地，他往左蹭蹭，我往右蹭蹭，好在相遇点两侧都是小车。我看到我的左侧后视镜距离面包车的左侧后视镜只有12厘米到15厘米，右侧后视镜距离停那儿那辆桑塔纳的右侧后视镜只有4厘米到9厘米。

过去之后我才纳过闷来，那叫后怕。天哪，我这是怎么过来的。。

继续在小区里逡巡。

逡巡许久，终于找到了一个理想的车位。停进去。右侧有矮墙，副驾门开不开。

老婆说，那我从你那边爬出去吧。

我说，别，我往回倒点。

老婆下车后，我又把车往边上靠了靠，停稳，轮回正，拉手刹，放P档，拔了钥匙，走人。

马尔代夫满月岛攻略

2011-11-20T15:18:00.001+08:00

11月15日
早上9点从T2航站楼出发，这个周二早高峰时间路面上的机场巴士不靠谱，于是7点半在北三环上打了个车，去三元桥换机场快轨。8点5分到机场柜台，被告知再有五分钟通道就关闭了，真悬！下午2点多到马累机场，出飞机差点被热死。老子还穿着秋裤，在赤道上忍了2个多小时才到房间。坐船到满月岛，前台交400$押金。喝了一种非常好喝的饮料，填了张表格。
先住沙滩屋，洗澡的地方是上半露天的，洗的慢了会晒伤。还有个阳台，阳台外有个水龙头可冲脚，因为从沙滩回来通常脚上都是沙子。傍晚到海滩走了走，订了第二天的海豚巡游。
16日
5点多起来看了个日出的尾巴，趁太阳不太猛，到处走走拍拍。金色的朝阳暖暖地打在沙滩上。
吃完早餐就回房睡觉。下午下海游了会儿泳，又睡了。
傍晚出海去看海豚巡游。没有见到海豚。日落很美，浪很大。
17日
搬家到水屋，有车接。水屋特别好，洗澡的地方有270°的视野看海景。从阳台可以直接下水。中饭后，我们就下海试了装备。
累了就回房睡觉，或者在阳台的躺椅上吹海风，吃苹果。
18日
又看了个日出的尾巴。上午在浅滩看到了小鲨鱼，还有上千条小鱼。
中午继续在阳台晒太阳，忽然想到一句词"海空凭栏阔"，不知出处。后来才知道是杜撰出来的，没有这样的说法。如果把"凭"和"栏"断开，把前者解释为"是因为"，那么这句词就有了含义：觉得海面空旷，是因为阳台宽敞。
下午在近海浮潜，拍水下小鱼。拍到了一些。
晚饭去吃了烛光晚餐。
19日
什么也没做，上午睡了个懒觉，中午交了房就坐船回机场了。

心得：
自始至终没觉得是出国，到处是中国人。还以为是到了三亚。
满月岛貌似由喜来登托管，国际化做的不错，几乎所有服务都有中文。那些不懂英语的富二代到了岛上也不会饿死。
自助餐貌似不用给小费。烛光晚餐和出海看海豚巡游直接记账就好，不用给现金。账单上已经包括了消费。
早餐相当丰盛，一坐下就有人来给你倒咖啡、橙汁。豆奶很好喝。各种面包点心粥。各种鸡蛋的做法：煎蛋、水煮蛋、蛋卷。各种果酱和水果。
浮潜装备因人而异，本人因为视力不好，带上面罩水下什么也看不见，只好换带度数的泳镜。气管的塑料味和海水的咸味让我恶心，在嘴里咬了一会儿就扔一边了。穿着救生衣，很难潜到水下太深的地方。大陆架很长，水也不深。脚蹼很好用，没有它基本上游的速度是0。光脚下水很危险，水下珊瑚轧脚。
水屋柱子上有一种奇怪的螺，不小心靠近柱子的时候碰到，手被剌出一个口子，直流血。
岛上只有一个地方可以免费上wifi，在Lobby旁边的网吧。
我们搭乘的是海航的飞机，往返各含两餐，一顿正餐有热米饭，一顿落地餐是三明治。回到北京是凌晨2点半，公共交通工具不要想了，直接打的回家。
自己带了一箱子吃的，因为酒店管4顿早餐，2顿晚餐，所以另外还有6顿饭要自己解决。
不用带水，房间里有瓶装水，还可以烧开水。
防晒霜不能省，我们每次下海前都抹，仍然会晒黑，不过不会晒伤。
不要为整天介睡觉而懊丧，心情好精神好是最重要的。这就是马尔代夫度假的意义。
如果平时苦日子过惯了，开始一两天会无法进入状态，适应不了这么慵懒并且奢靡的生活。带着一点批判去享受它，慢慢就好了。回到北京之后，你会忽然发现，哇塞自己的执行力好强。

测试MathJax

2011-11-12T21:57:00.001+08:00

下面的欧拉公式可以在网页上自动渲染。
$$e^{2\pi i}=1$$
这是因为我在Blog模板的<head>中加入了公式渲染引擎

<script type="text/javascript" src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>

还可以使用inline的公式渲染。例如$\frac{\partial^2 \phi}{\partial x \partial y} = \frac{\partial^2 \phi}{\partial y \partial x}$呵呵。
只需要加一行<script type="text/x-mathjax-config"> MathJax.Hub.Config({tex2jax: {inlineMath: [['$','$']]}});</script>

2-SAT

2011-11-11T23:18:00.001+08:00

总算release了，释放一下，找个简单的问题自娱自乐一把。

n个变量组成m个子句的句子，其中每个子句中恰有两个变量用或连接，每个变量可以替换成他自身的取反，子句之间用且连接。

给一个这样的句子s，求一组真值分配。

算法：

如果存在一个不成立的子句，就反转其中一个变量。

反复为之，直至超过循环的最大次数。

（听到一个声音说：这叫什么狗屁算法，这不就是在瞎蒙么！）

实现：

import random, re

sentence = lambda n, m: " and ".join(["(%sx%d or %sx%d)" % (

random.choice(["not ", ""]), random.randint(1, n),

random.choice(["not ", ""]), random.randint(1, n)) for i in xrange(m)])

def detect_self_contradict(sub):

(i, j) = re.match(".*?(\d+).*?(\d+)", sub).groups()

ret = i == j and sub.count("not") == 1

if ret: print i

return ret

def solve(s, a):

for _i in xrange(10000):

print _i, a

exec(";".join(["x%d=%s" % (i+1, a[i]) for i in xrange(len(a))]))

if eval(s): return a

for sub in s.split("and"):

if not eval(sub):

i = int(random.choice(re.match(".*?(\d+).*?(\d+)", sub).groups())) - 1

a[i] = not a[i]

break

return "failed to solve %s" % s

a=[True]*20

s=sentence(len(a), 24)

if any([detect_self_contradict(sub) for sub in s.split("and")]):

print "no way! %s" % s

exit()

print s

print solve(s,a)

一次试验：

python 2sat.py

(x19 or x6) and (not x1 or x20) and (x5 or not x20) and (x14 or not x20) and (x4 or x2) and (not x19 or x2) and (not x15 or not x6) and (x8 or not x15) and (x9 or not x17) and (not x11 or x4) and (not x16 or not x4) and (not x17 or not x11) and (x20 or x20) and (not x12 or not x16) and (x13 or not x10) and (x3 or not x17) and (x17 or x14) and (not x10 or not x6) and (not x13 or x19) and (x3 or x1) and (not x5 or x8) and (x15 or not x3) and (x1 or x1) and (x3 or x4)

0 [True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True]

1 [True, True, True, True, True, False, True, True, True, True, True, True, True, True, True, True, True, True, True, True]

2 [True, True, True, False, True, False, True, True, True, True, True, True, True, True, True, True, True, True, True, True]

3 [True, True, True, False, True, False, True, True, True, True, False, True, True, True, True, True, True, True, True, True]

4 [True, True, True, False, True, False, True, True, True, True, False, True, True, True, True, False, True, True, True, True]

[True, True, True, False, True, False, True, True, True, True, False, True, True, True, True, False, True, True, True, True]

我试验了很多次，发现凡是有解的时候，都不需要迭代很多次。

再分析：

如果存在一个满足句子s的真值分配A，那么令当前分配X[k]与A相同的变量个数为f(A, X[k])。

每次循环，至少有一半的概率f(A,X[k+1])=f(A, X[k]) + 1

为了证明这一点，假设一个当前未满足的子句X或Y，A以概率p取T,T，以1-p的概率取T,F或F,T。当前取值X[k]显然为F,F，那么f(A,X[k+1])=f(A, X[k]) + 1的概率依全概率公式得

p*1+(1-p)*0.5=0.5+0.5p>=0.5

证毕

f(A,X[k+1])=n时，我们就求得了一个真值分配。

这个试验告诉我们一个道理：如果你相信存在一个真理，并且在实践中发现每次尝试的试错概率不超过一半，并且每次尝试成功都可以向真理靠近一个常数距离，那么，你就大胆地去不断尝试吧！（听到一个声音说：这叫什么狗屁道理。。。）

1977

2011-10-24T21:48:00.001+08:00

1977年什么也没有发生。今天，2011-10-24，也不是什么特殊的日子。

但是当混用python和awk时，神奇的事情就这样发生了。

os.system("awk '{print %s}' a.txt" % str(datetime.date.today()))

输出

1977

坑爹啊。。。

Fwd: 分析笑话：数学系有三个班

2011-10-16T00:57:00.001+08:00

“数学系有三个班”是真的，可我从来没听过这样有趣的对话：

http://www.guokr.com/post/67116/

先分析一下原版的吧，说来惭愧，还是老婆给我讲明白的呢。

数学系一共3个班。

男生：你是3班的么？

女生：原来你是2班的啊！

男生：你是1班的啊！

初始：

男生知道自己是2班的，男生知道女生不是自己班的，所以男生知道女生是1班或3班的

女生知道自己是1班的，女生知道男生不是自己班的，所以女生知道男生是2班或3班的

第一句话：

男生问女生是不是3班的，相当于宣布自己不是3班的，因为如果男生是3班的，他就不用问这个问题了。这里假设这两个人认识自己班的所有人，并且彼此知道这一点。

所以女生知道男生是2班的了。

第二句话：

女生说“原来你是2班的啊”，相当于宣布自己之前不知道男生是2班的，现在得到足够的信息知道了。

第三句话：

男生是如何推断出女生是1班的呢？

假设女生是3班的，那么女生只能知道男生是1班或2班的。那么男生告诉女生自己不是3班的并没有给女生提供任何新的信息。那样的话，女生就不足以推断出男生是2班的。这和第二句话不符。

所以女生只能是1班的。

趁老婆还在玩iCloud没去睡，再分析一个鬼谷子版的吧

一天，鬼谷子随意从2-99中选取了两个数。他把这两个数的和告诉了庞涓，把这两个数

的乘积告诉了孙膑。但孙膑和庞涓彼此不知到对方得到的数。第二天，庞涓很有自信的

对孙膑说：虽然我不知到这两个数是什麽，但我知道你一定也不知道。随后，孙膑说：

那我知道了。庞涓说：那我也知道了。

这其实就是YaoQiZhi老师课上讲的Mr. S & Mr. M的故事，当时觉得好难啊，现在学了点python，觉得暴力破解法很直接，code也不过一屏。

from itertools import product

from collections import defaultdict

n=552

a0 = set([p for p in product(range(2,n), range(2,n)) if p[0]<p[1]])

print len(a0)

s=defaultdict(set)

m=defaultdict(set)

for (x,y) in a0:

s[x+y].add((x,y))

m[x*y].add((x,y))

print len(s)

print len(m)

#I don't know, you don't know

a2=set([(x,y) for (x,y) in a0 if len(s[x+y])>=2 and len(m[x*y])>=2])

print len(a2)

#I know you don't know

a3=set([(x,y) for (x,y) in a2 if all([len(m[p*q])>=2 for (p,q) in s[x+y]])])

print len(a3)

#you know

a4=set([(x,y) for (x,y) in a3 if len(a3 & m[x*y])==1])

print len(a4)

#I know

a5=set([(x,y) for (x,y) in a4 if len(a4 & s[x+y])==1])

print a5

ytwang$ time python v.py

150975

1097

77432

105205

4753

2058

set([(4, 13)])

real 0m18.002s

user 0m17.919s

sys 0m0.079s

计算复杂度的瓶颈在求a3，也就是“我知道你不知道”这句话蕴含了绝大多数信息量。

Fwd: 内省+变长参数的lambda表达式

2011-10-10T23:59:00.001+08:00

最近看某同事写的python代码写的很漂亮，不禁羡慕。

恰好另一个同事买了本python教程，便拿来学学。

写了个练习。

From Evernote:

内省+变长参数的lambda表达式

class Client:

def __init__(self, endpoint):

self.endpoint = endpoint

def wget(self, query):

url = "http://%s/%s" % (self.endpoint, query)

import urllib

#todo wget and return response

return url

def __getattr__(self, attr):

return lambda **kargs : self.wget(attr + "?" + "&".join(["%s=%s" % (k, kargs[k]) for k in kargs]))

if __name__ == '__main__':

a = Client("baidu.com")

print(a.q(keyword = "s"))

print(a.q(keyword = "s", charset = "UTF-8"))

结果

http://baidu.com/q?keyword=s

http://baidu.com/q?charset=UTF-8&keyword=s

滑雪板理论

2011-10-06T17:25:00.001+08:00

我假设你滑的是双板。

当你倒在雪地里，发现两个雪板都不见了，那是最开心的。你可以立刻起身，东张西望，去捡你的雪板。你可以先捡左脚的雪板，再捡右脚的雪板，也可以先捡右脚的雪板，再捡左脚的雪板。你可以在雪地上自由行走，比穿着雪板的时候自由得多。只要不被后面失速的其他人再次撞到就好。

当你倒在雪地里，发现左脚的雪板不见了，右脚的雪板还在脚上，是很郁闷的。你的右脚别在那里，不仅不能发力，而且不能提供合适摩擦力。单靠左腿的力量又不足以起身。你有两种选择，要么等人来拉你一把，当然这意味这你将命运交给了运气和泊松过程，要么你用撑杆把右脚的雪板脱掉，如果你的撑杆没扔出去的话。如果选择后者，相当于你把问题规约成了两个雪板都不见了的情形。你可能知道这个问题解起来更容易，也可能是歪打正着。

当你倒在雪地里，发现左脚得雪板还在脚上，右脚得雪板却不见了，也是很郁闷的。你右腿的力量可能强些，因为你踢足球，但是单凭右腿起身仍有些费力。你无非多一种选择：拼了命地站起来。虽然教练不是这样教你的。

当你倒在雪地里，发现两只雪板都TM还在脚上，你可以像个婴儿一样哭了。你觉得很亏，因为不仅要费力把两只雪板一只只脱掉，从这白茫茫的斜坡上站起来，还要把这双雪板摆成平行于等高线的角度，小心穿进去。你觉得委屈，因为这比掉了一只雪板还要费事，却在外人眼里看来还挺幸运的。

然而，倒在雪地里，还是很舒服的，至少它是一个稳态。

滑雪最不舒服的时候，是你发现自己已经超速，失去了对身体平衡和速度的控制的时候。

你喜欢一切尽在掌控，你后悔干嘛要跑到中级道上来吃饱了撑的没事找罪受。

你感到身体正在向左倾斜，同时向右旋转，你努力回忆教学视频上的内容，大脑却一片空白。

你感受到梯度下降，同时感受到无法忽略的水平位移。

你盼望着奇迹发生，身体赶快恢复平衡，严格按梯度下降，然后安全地在紧支集外面着陆。

但是力不从心，你知道你坚持不到平地了，你重心在下沉，心里盼望着赶快进入另一种稳态，哪怕摔得很难看。

你盼望着摔得两只雪板都不见了才好，如果不能，摔得两只雪板都在脚上也行，无非费点事，脱掉爬起来再穿上。

你想赶快结束这个痛苦的"要摔还没摔，没摔就要摔"的该死状态。

你开始破罐破摔，给自己捣乱，右脚外脚掌向里扣，想引入一个扰动，把自己弄倒。

然而速度太快了，这点扰动简直蚍蜉撼树，无济于事。

你重心继续下沉，两腿僵硬，死盯着雪板的前尖，越来越害怕，心里除了念着"wo cao wo cao"就已经没有了任何建设性意见。

就在这时，你发现你已经开始减速了。抬头看看。是的，恭喜你，你已经到平地了。

这时你腿一软，pia，摔那儿了。

这不是看辛普森。

上帝存在，帕斯卡说的

2011-10-04T10:35:00.001+08:00

“我们既不知道上帝的存在，也不知道上帝的本质。然而我们将倾向于哪一边呢？……，这里进行的是一场赌博，…… 让我们来权衡一下在上帝存在的赌注中的得失。让我们估计这两种可能性，如果你赢了，你赢得所有；如果你输了，你却一无所失。因此，你就不必迟疑去赌上帝的存在吧。”

这是帕斯卡在他的《思想录》中说的。

这其实是说一个条件期望。在相信上帝存在的条件下，期望收益是正的。在不相信上帝存在的条件下，期望收益是0。中国民间关于鬼神等无法实证的事物是否存在，流传着一句神似的说法：宁可信其有。

http://mathforum.org/isaac/problems/prob1.html

在那个年代，不奇怪，帕斯卡身边有个和他一样对古典概率论痴迷的神经质，比如费马。

一个周日的下午，帕斯卡和费马在巴黎的一个咖啡馆觉得无聊，于是玩抛硬币。正面费马得一分，背面帕斯卡得一分，谁先得到10分谁赢，输的买单。

15次后，费马得8分，帕斯卡得7分。这时他们玩烦了，决定提前结束这个无聊的游戏。他们不是中国人，不会争着去买单。也不是美国人，不会evenly share the bill。

他们太无聊了，决定用一种他们认为公平的方式来结束这个周末。

于是帕斯卡管邻桌的教授借来了笔和纸。帕斯卡和费马在几轮争论后，最终接受了费马的解法：

游戏将在至多4轮后结束，只需要将2^4种可能的序列枚举出来即可。h表示正面，t表示背面，*表示费马赢。

`h h h h *`	`h h h t *`	`h h t h *`	`h h t t *`
`h t h h *`	`h t h t *`	`h t t h *`	`h t t t`
`t h h h *`	`t h h t *`	`t h t h *`	`t h t t`
`t t h h *`	`t t h t`	`t t t h`	`t t t t`

显然这16种序列出现的可能性相同，所以11/16的概率是帕斯卡买单了。因此，公平起见，帕斯卡付账单的11/16，费马付剩下的。

帕斯卡当时不会想到，这其实是金融工程中经典的风险资产定价问题。

假设账单价格为1法郎，对于帕斯卡来说，他收益是一杯咖啡，成本有两种可能，一种是自己请客，也就是赌输了，支付1法郎，另一种是费马请客，自己一分钱不掏。前者的概率为11/16。帕斯卡可以选择一路赌下去，碰碰运气，不过他选择了无风险投资，避免自己全盘皆输，因为局势已经不利于自己了。费马也同意这个无风险的方案，因为他可以立即套现，什么也没有做就立即赚到1/2-5/16=3/16法郎，落袋为安。

为什么帕斯卡此时必须支付11/16法郎？

如果帕斯卡支付的少于11/16法郎，那么费马将选择继续赌下去，因为他可以算出来，在当前情形下，赌下去期望的花费更少。

如果帕斯卡支付的多于11/16法郎，那么帕斯卡还不如赌下去。

所以帕斯卡没有其他选择，只能支付他的条件期望。

当然，这里假设帕斯卡和费马都是风险中性的投资者。要是其中一个或者两个是风险厌恶型或风险偏好型，就另当别论了。

“我正在说的这句话是谎话”的三种等价说法

2011-09-02T22:08:00.001+08:00

一

什么是真理？

真理就是永远正确的。

什么是永远正确的？

没有什么是永远正确的。

二

我可以满足你三个愿望。

太好了，我的第一个愿望是：你可以满足我1000个愿望。

对不起，这是一个元愿望，我不能满足元愿望。现在你还有两个愿望。

好吧，我的第二个愿望是：你能够满足元愿望。

对不起，这是一个元元愿望，即使我能够满足元愿望，我也不能满足元元愿望。顺便提一句，我不能满足任意阶的元愿望，你懂我的意思，请珍惜你的最后一个愿望。

谢谢。让我来告诉你我的最后一个愿望：我希望我的愿望不被满足。

（http://www.eaglefantasy.com/archives/1209）

三

我：你为什么不高兴？

老婆：（摇头）

我：是因为男主角死掉了么？

老婆：（摇头）

我：是因为我不关心你么？

老婆：（点头）

我：我哪儿不关心你了？

老婆：（摇头）

我：是因为我老加班么？

老婆：（点头）

我：可是我现在不是正在关心你么？

老婆：（摇头）

我：你现在是在摇头么？

老婆：（顿了两秒钟，忽然转过来朝我就打）

修改注册表防止sftpdrive试用版过期

2011-08-28T19:28:00.001+08:00

sftpdirve只能试用45天，到期后卸载重装仍然提示过期不能使用。断言有残留注册表项记录了软件过期状态。

解决方案：

1、卸载sftpdrive

2、创建一个t.vbs的文本文件，粘贴如下内容并双击该文件运行

WScript.CreateObject("WScript.Shell").RegDelete("HKEY_LOCAL_MACHINE\SOFTWARE\SftpDrive\")
WScript.Echo("Done!")

3、重新安装sftpdrive

Enjoy~

期权定价

2011-08-14T15:45:00.001+08:00

读研的时候，很多同学都在研究期权定价，后来去了金融、保险、证券之类的地方，挣得不少，当然喝酒是免不了的。

最近迫于某些方面的压力，想换份工作，自然想到这些同学当年学的东西，便也弄了本书来看。还挺有趣的。

第一章中一个漂亮的结论：

在一个简单的模型中，期权的价格与股票上涨或下跌的概率无关。

假设：

银行存款利率确定并且已知为a>0，即今天存入1元，明天可以取出1＋a元

世界上仅有一只股票，收益率随机：今天买1元的股票，明天卖得1＋u元的概率为p，卖得1+d元的概率为1－p，d<u，0<p<1，别无它价

存款和股票的交易量可以为任意实数：连续、无界、可负（卖空）。

市场不允许没有初始投资的五风险利润－－无套利原则－－不存在空手套白狼的投资方案。

今天花c元买1份执行价为1＋s的看涨期权，d<s<u，明天得到u-s元的概率为p，什么也得不到的概率为1-p。

定理：

c＝x+y

其中(x, y)是如下方程组的解

x(1+a)+y(1+u)=u-s

x(1+a)+y(1+d)=0

证明：

首先方程组的解存在唯一，因为矩阵

1+a, 1+u

1+a, 1+d

的行列式＝(1+a)(d-u)<0

其次，

若c>x+y，则今日卖空1份期权，套现c元以(x, y)投资银行和股市，剩下c-x-y>0落入腰包。明日若股市上涨，则兑现期权付出u-s元，在银行和股市的资产x(1+a)+y(1+u)恰好可以填平这个窟窿，若明日股市下跌，则在银行和股市的资产x(1+a)+y(1+d)＝0亦无需兑现期权。总之无本净赚c-x-y。

若c<x+y，则今日以(-x, -y)卖空银行和股市套现以买入1份期权，剩下x+y-c>0落入腰包。若明日股市上涨，则期权收入u-s元，恰好用来还清在银行和股市的负债x(1+a)+y(1+u)，若明日股市下跌，则在银行和股市的资产x(1+a)+y(1+d)＝0亦无需兑现期权。总之无本净赚x+y-c。

自反+反向传递=>完全

2011-05-21T01:43:00.001+08:00

今天终于证明了这样一个定理：偏序是全序，当且仅当可反向传递。藉此可以理偏序和全序的本质差异。

回忆
偏序=传递+反对称+自反
全序=传递+反对称+完全

E是一个集合，≤是E上的一个二元关系，称
自反：a≤a，对于任意a∈E
完全：a≤b或b≤a，对于任意a, b∈E（显然完全蕴含自反）
反对称：若a≤b且b≤a，则a=b
反向传递：a≤b不成立且b≤c不成立，则a≤c也不成立，对于任意a, b, c∈E

求证，若≤自反和反向传递，则≤完全
证明：假设不然，即存在a, b∈E，使得a≤b不成立且b≤a不成立，则由反向传递，得a≤a不成立，与≤自反矛盾

求证：若≤完全、反对称、传递，则≤反向传递
证明：对于任意a, b∈E，若a≤b不成立且b≤c不成立，
由≤完全，有b≤a且c≤b，
由≤传递，有c≤a
假设a≤c成立，则
由≤反对称，有a=c，
代入b≤a且c≤b，得到b≤a≤b，
由≤b反对称，有a=b
代入a≤b不成立，得到a≤a不成立
与自反矛盾，故而假设不成立
■

也发美国商务签证面经

2011-04-15T22:05:00.001+08:00

公司给约了中午11点15的面签，三个同事一起去的，临出门前听说还要一张51mm×51mm的B1签证专用相片，惊慌失措，没有！
后来又听说如果网上提交的电子版照片足够清晰的话，没有纸板相片也可以，就直接去了。
上午十点从公司打车过去，不到半个小时。美国大使馆门前有很多人，所以下车后看哪里人多就好了。司机、咨询机构的服务员和路边的清洁工都知道。我们傻呵呵的还差点进了马来西亚驻京办。在没风的天气里，马来西亚国旗和美国国旗很难分辨。
排队5分钟后，院子大门口警卫会检查护照和确认页。
再排队10分钟进建筑物的门，建筑物门前的工作人员会每次放进去约10个人，门前有垃圾桶，可以扔路边被塞的小广告。
进去之后第一道关卡，出示护照和确认页。
往里走第二道关卡，安检，存手机，他会给你个牌。
再往里走第三关，检查材料：护照、确认页、邀请函。确认页个人信息旁边要用中文写名字、电码、公司地址、公司名称。
这关比较有趣，每个时刻会随机出现若干可以服务的窗口，其他窗口不服务，但是可以把你重定向到某个可服务的窗口。每个窗口服务6个人之后，会停止服务，然后把这6个人编成一队。注意如果你所在的窗口服务3、4个人（包括你）之后，这个窗口没人来了，而其他某个窗口排了队，你可以去叫几个人过来，（人肉优化），这样可以加速你所在的窗口编队，尽早进入下一关。队长会领到一个有颜色的牌子。完成编队的窗口会恢复服务，如果他不是去忙别的的话。这关的窗口都是中国人在服务。大约耗时15分钟。
第四关：按手印。跟队去按手印，左手三个手指、右手三个手指、两个大拇指，如果没按清楚他会提示你重新按。这个窗口是美国人服务，但是会讲中文。你不用说什么话，照他说的做就好了。排队等候大约15分钟。
第五关：面谈。排队等候大约70分钟。典型的单队列多窗口服务，效率已经很高了。因为到了中午吃饭时间，窗口可能会随时突然关闭，这时你将被重定向到其他窗口而不用像火车站一样重新排70分钟的队。问题内容因人而异，他可能会对你的公司业务感兴趣，问得很具体，我准备了一些面签常见问题，不过没想到他问：你们公司是如何帮助客户挣钱的？哇塞，就Business Model那个Page没好好看，太失败了。只好凭着印象胡扯了一通，不过他只在我打奔儿的时候盯着我看了一会儿，真正等我扯上了，他就一边点头一边忙别的去了。然后还问了家里的很多事情，你爸干什么的，在哪里，你妈干什么的，在哪里，你有没有移民倾向。问得很直接：你是否打算移民。直接回答否就好了。面试官是美国人，但是会讲中文，我前面的几个人一直说中文，也可以，但是如果是去留学就必须说英文。我因为之前背了一些面签常见问题都是用英文准备的，所以就赶上背过的问题就用英文，赶上没背过的问题就用中文，中英文夹杂着糊弄过去了。面试官似乎对此毫不介意。如果读过工程物理或空间科学的，可能会被Check，好像怎么解释也不行，被Check就是一审不过，回去等消息，可能过，也可能要重申。简历有可能会被看，不过我的没被看。51mm×51mm的相片这一路也没人找我要过。
出来拿牌子取手机。
第六关：去邮局办快递或自取。自取可以多个人一起取，比较便宜。可能是因为已经到了中午一点，几乎没用排队。
出来就是亮马桥地铁口，我们饿得眼冒金星没看见，又打了个的，回公司了。

2011推荐系统高峰论坛――参后感

2011-03-06T23:34:00.001+08:00

最近家事繁多，好久没更新Blog了。
今天参加
http://www.resysforum.org/
2011推荐系统高峰论坛，说来小有周折。消息刚放出来不久，同事zhangyifei就组织我们几个FreeWheeler报名了，结果被拒了。然后前几天我忽然收到resys内部人士邀请，让我去了。猜着可能是看在我是resys第一次线下活动的speaker之一的面子上吧。总之作为已经成为resys领域outer的我很是喜出望外，今天一大早，就背着小板凳去了。想必像我这样"走后门"进去的，顶多是张站票罢。不料竟然还有地方坐，让我很是Happy。在我之后进来的一些人，站在会场后面，让我很如坐针毡，很想把小板凳掏出来给他们，可是给谁不给谁呢，都不认识，再招他们互相谦让得打起架来就不好了，而且做小板凳上恐怕就看不见大屏幕了。算了，还是悄悄坐在我自己的位子上吧。

据谷文栋说，今天地大这场是resys第五次线下活动了，有1000多人报名，而第一次活动只有30来个人。第一次绝大多数是热心搞算法的草根或学生。这次，世界顶级学者、搞工程的、搞产品设计的、搞系统运维的，"啥都有"。

Koren介绍了Y!的一些最新研究，是今天最好的一个lecture。
分析了一下Netflix样本方差来源，超过一半无法用模型解释，可解释的成分中，"公认怎么样+个体打分松紧程度"占了一大部分。
推荐系统现在的发展，越来越多的和用户研究掺和在一起了。User Research原本是一个独立的功能机构，一提它，大多数人的第一反应是眼动仪和心理学专家。现在推荐系统揉进去，不仅多了一些统计学习的东西，而且日益工程化了。
用户研究，很直接的一个手段就是问卷调查。如果受某些限制只能问K个形如"你喜欢item i吗"这样问题，那么这K个item如何选取？Koren给出了一个cost function，使得给定的推荐算法在这组item上均方根误差最小，也就是说，当前模型能够较为准确刻画这组item，这样问用户，才能较为准确的理解用户的喜好。
另一个方法是决策树，它每次只问用户一个选择题，等用户选择完以后，根据用户的回答再问第二个问题，这样更为有效。
这两个办法都无法回避的问题是：用户可以回答"没听说过"，这样你就没办法知道用户的喜不喜欢，就白问了。所以不能轻易问太小众的item。
推荐系统应用到"问答社区"，就是把好问题推送给能高质量回答它的人。
另外，它提出置信度作为评价推荐系统的一个指标，可以改善推荐质量，例如可以对置信度低的预测点给予一定的"高斯模糊"。而置信度的计算，其实每一种模型都有给出，可惜人们往往无视他存在的意义。

淘宝作为spnsor有两场lecture，第一场比较概括，有三点引起了我的思考。
1、淘宝在推自动化数据运行，逐步替代人肉管理商品目录。随着广告市场的发展，我认为以后ad trafficker也有自动化的趋势，一些智能工具的价值将得到凸显。
2、数据可视化和交互是数据挖掘的重要手段。我认为，他不仅是数据探索的工具，更是结果表示的方法。
3、淘宝今年成立了研究院。没搜到，不知为何物，关注中。

淘宝第二场
1、推荐与广告的关系：掐。推荐系统是营销部做的，旨在提高更多long tail卖家的关注度，改善消费者的购物流程，优化整站生态。广告是广告部做的，是拿卖家钱，就要像狗一样给卖家效力，处处为给钱的人着想。
2、用户在不同阶段有不同的推荐需求，按时间顺序是发散、收敛、再发散的过程
四处闲逛阶段：各种大跨度的试探性推荐，揣测闲逛者的消费意图
精益求精阶段：发现用户开始货比三家，猎物相似了，那就要把推荐的范围逐渐缩小了，加速消费者下订单、支付的进程，说，嗨哥们儿别犹豫了，就它了！这时候再给他推荐一些天马行空的东西，没准人家一走神点走了，就不买这个了，去东逛逛西逛逛，一拍大腿啊呀三个小时过去了，人家说我得干点正事了，得，把浏览器关了。
意犹未尽阶段：这时候再开展你"啤酒=>尿布"的数据挖掘伎俩，推荐些哪怕看似不想关，却在统计数字上有关联的商品，即使推荐的莫名其妙，也无伤大雅，反正最坏情况就是就此打住，没准弄好了还能神奇的再下一城。
3、实时训练与预测。在燕莎、在双安、在大悦城，好的导购能通过消费者嘴角一个0.13秒的翕动，在0.013秒内判断其消费意图并推荐合适的商品，电子商务的推荐算法必须要向这样的导购一样。慢了，消费者就把tab关了，或者已经买完走了玩别的去了。

百度，大纲级
很适合我的口味，讲了很多我很有兴趣的学术界研究的课题，比如Learning to Rank。
演讲者高屋建瓴地将机器学习的问题归为三个：回归、分类和密度估计。听着很像统计学习的口吻，不知道工业界搞商业智能的兄弟们知不知道密度估计有什么应用，反正孤陋寡闻的我是闻所未闻。
五个成分：数据、机器表示、目标函数、cost function和evaluator
现实世界中的学习，几乎都是PU learning（partial learning），大量数据缺失，大量未标定数据。
一个我不太熟悉的模型：Prob Latent Topic。好像是十年前特别火的一个模型，读研时看文献综述的时候一坨一坨这玩意儿，但是在netflix price中表现很差，也可能是我当时实现的不好，又没机器去没日没夜的调参。总之，玩不转。三个方法：EM/Gibbs Sampling/梯度下降。
Learning to Rank三个方法：SVM、Boost和NN。可惜因为空间复杂度太高，没什么工业应用。
一个观点：ReSys以后将成为互联网应用的基础设施。哇塞，太牛了，听的让我热血沸腾，我以前只听说有人争论推荐系统究竟是个feature还是个product。基础设施是什么？在我脑海中就是磁盘阵列、路由器之类能被会计算作固定资产的东西。以后推荐系统也是，太牛了。
好的推荐系统的四要素：UI/UX占40%（难怪PM有饭吃），Data占30%，领域知识占20%，算法占10%。听的在坐的学生们心灰意冷，纷纷微博上墙感慨人生。
推荐系统研究的未来发展方向：领域专家系统、迁移学习和Temporal-bias。
还提到了Behavior Targeting：把恰当的广告推荐给对他最有用的用户。在我刚踏出校门的时候，我觉得这简直太美妙了。现在的我觉得，这想法太前卫了。

前Google研究员
将MapReduce在机器学习中的应用，比较工程。对于一个C++程序员来说，听的非常实在，没有一句话是飘在空中的，都是code。
学到一点：用BDB的B+树做外排，太毁硬盘了，因为要频繁读写一个Block。应该用内排加外部merge sort。

Hulu
这个之前聊过，很多都听过了。只记了一些新鲜的。
Hulu没有Live Sport。
内容分两种，On air show和lib show，On air show在Hulu上只能放最近5集，因为之后DVD上市，要是还在Hulu上能看到，谁还去买？
Seasonality是时间序列分析的重要话题，Hulu的Seasonality特征鲜明：暑期、圣诞是On air show的两个波谷。连续剧每周出一集。
Catch-up效应：show的日访问量关于days after release基本上是负指数曲线。
广告效果评估：Brand recall。奇怪FreeWheel研发这边怎么没人提过这个，估计美国那边在弄吧。
预测匿名流量的Demograph info，用于Demographic targeting。
Demographic targeting是广告主最容易接受的targeting方式，其次是Geo, daypart, category。。。
Behavior targeting的逻辑：请第三方来证明Hulu存在用户群A，具有购买商品P的习惯。然后Hulu通过某种方式，将用户群A具象化到一个形式化的targeting范围T里。向商品P的卖家或品牌广告商推荐：将广告targeting在T上吧！
在场有人质疑：品牌广告商需要支付更高的CPM来通过Behavior targeting保持现有用户群吗？以我对广告商的理解，是需要的，因为以他们的非理性程度和对人性深刻的把握，足以做出任何在常人看来是蠢事的事。

A spike

2011-02-20T19:40:00.001+08:00

马尔科夫链蒙特卡洛（二）

2011-01-11T22:56:00.001+08:00

马尔科夫链蒙特卡洛（二）

来，做一道书后习题：用MCMC生成服从Poisson分布的随机数。

习题嘛，瞎折腾，本来R就有现成的函数：rpois
既然MCMC，首先要构造一个转移规则。考虑射线上的随机游动。从0出发，必然走到1，从k出发(k>0)，一半的概率向左右一半的概率向右走。

计算接受概率：

其中

其中k>1。边界0和1处需要稍微特殊处理一下。注意到，概率密度函数求比值之后，我们再也不用算指数和阶乘了。剩下的事情，就是不停地走，走每一步时只需要回答：

1）往哪边走？---用随机游动的trial distribution：一半的概率向左右一半的概率向右走

2）真的要走到那里去吗？---使用上面算的接受概率（MCMC的核心）

下面是R源代码

p=function(s){
    runif(1) <= s
}
#pi_rate(k,lambda)=pi(k,lambda)/pi(k+1,lambda)
pi_rate=function(k,lambda){
    (k+1)/lambda
}
accept_walk_left_from=function(k,lambda){
    if (k==1)
        min(1, 2 / lambda)
    else
        min(1, pi_rate(k - 1, lambda))
}
accept_walk_right_from=function(k, lambda){
    if (k==0)
        min(1, lambda / 2)
    else
        min(1, 1 / pi_rate(k, lambda))
}
transition_from=function(k, lambda){
    if (k==0){
        if (p(accept_walk_right_from(k, lambda)))
            1
        else
            0
    }
    else{
        if (p(0.5)){
            if (p(accept_walk_right_from(k, lambda)))
                k + 1
            else
                k
        }else{
            if (p(accept_walk_left_from(k, lambda)))
                k - 1
            else
                k
        }
    }
}
mcmc_poisson=function(n,lambda){
    f=0
    for(i in 1:n){
        f[i+1]=transition_from(f[i], lambda)
    }
    f
}

samples=tail(mcmc_poisson(20000, 5), 10000)
hist(samples)

下面是一次实验结果的直方图，看起来和Poisson分布很像吧

马尔科夫链蒙特卡洛（一）

2011-01-08T21:27:00.001+08:00

马尔科夫链蒙特卡洛（一）

传统的马尔科夫链研究，是给定转移规则，关注平稳分布，而马尔科夫链蒙特卡洛研究中，我们知道想要的平稳分布，任务是构造有效的转移规则，使得它的平稳分布就是想要的那个。

记为目标平稳分布。Matropolis-Hastings算法的一般框架是：寻找一个函数T(x, y)满足T(x,y)>0当且仅当T(y,x)>0，

1）给定当前第t步的状态x[t]，依照分布T(x[t], y)随机采y

2）下一步状态x[t+1]以概率r(x[t],y)等于y（接受，转移），以剩余的概率等于x[t]（拒绝，留下）

其中,里面的分子是某个对称函数。

则第t步从状态x转移到y的概率

由分子的对称性得到

这说明上述算法的得到的转移轨迹是可逆的，或者说具备微观平衡的性质，下面证明，这个马尔科夫链以为不变分布。

也就是说，如果这一时刻的状态分布是，那么下一时刻的状态分布还是。

由经典马尔科夫链理论，既然它是非周期正常返的，而且以为不变分布，那么它必然收敛与，即成为平稳分布。

于是我们得到了采样方法：从某一步开始（即预热后），状态序列x[n].x[n+1],...就是我们需要的样本序列。

通常取

特别地，当T(x,y)也是对称函数时，它进一步退化成

如果你还记得重要性采样，是不是觉得这个接受概率的式子很眼熟？

目标分布不好算，但是两个状态的概率之比好算，就交给接受概率好了。这是蒙特卡洛的基本初衷：好算就算，不好算就狂试。

注意MCMC得到的样本的两个弱点：

1、方差大

2、不独立

遍历trie

2011-01-07T22:40:00.001+08:00

当时在网络所实习的时候，还没有学过数据结构，觉得trie是一种很神奇的东西。现在想来，tree作为用来查找的结构，要求元素构成全序即可，而trie还要求元素构成字典序。
http://en.wikipedia.org/wiki/Trie

通常trie是用来查找的，怎么遍历它呢？
假设假设某个国家实行实名制上网，你要负责设计一个IP(v4)到身份证号的字典，而且被要求只能用trie，而且是定义成这样的：

struct Trie{
Trie * left;
Trie * right;
}

根节点指向一个Trie数组的首地址。查询时指向左孩子，如果待查询IP的二进制表示最高位是0，否则指向右孩子。如果孩子结点地址不幸指向数组件外面，则表明那不是文件地址，而是身份证号+数组的字节数，这个人拥有的IP段是他从根节点走到这里沿途收集的0-1串作为IP二进制表示的前缀的所有IP。好吧，我应该事先声明，这个国家人口稀少，以至于一个人可以拥有多个全球互联网上的IP段，而且为了方便市民记忆，身份证号也被设计的也很短。
遍历：
调用p(&data[0], "")
其中
p(Trie* node, string m)
{
if node 指向超出数组边界：打印身份证号为((int)arr)-文件长度, 他的IP段(二进制表示)有m+(32-m.length)个'0' ~ m+(32-m.length)个'1'，return
p(node->left, m+‘0’)
p(node->right, m+'1')
}

谨以此文献给即将用尽的IPv4地址和囤积了大量地址空间的IP黄牛们。
强烈推荐一篇用时间序列方法研究IPv4何时耗尽的论文，十分有趣：http://www.potaroo.net/tools/ipv4/#r5

两个不幸等价的随机数生成器

2011-01-02T17:14:00.001+08:00

中午接到电话，研究生同寝的哥们儿回北京了。于是一起去平安里吃烤鱼，一边吃鱼，一边研究一个随机数生成器的问题，一个学生问他的。
随机数生成器A：1/2的概率返回1，1/2的概率返回0。
随机数生成器B：先以[a, 1-a]均匀分布生成一个p，然后以p的概率返回1，1-p的概率返回0。其中0<a<0.5。
问题是，若只能观测这两个随机数生成器的返回结果，如何区分二者？
他说，好像无法区分，因为计算发现二者返回值的分布函数相同，却不知如何证明。
我想：样本方差相同么？后验概率相同么？回家一算，发现很不幸，也都相同。
可能这个题当初是考察离散、连续混合情形下的全概率公式的使用吧，确实在推导过程中容易陷入符号上的混乱。
用分部积分公式可以将这个结论推广：如果p是连续型随机变量，分布函数为F(t)，那么上面两个随机数生成器等价当且仅当F(t)在[a, 1-a]上的积分=1/2-a

给我碗儿水喝

2011-01-02T00:47:00.001+08:00

今天元旦，晚上去工体听了崔健摇滚交响音乐会。听众已经不再年轻，但是能听出来听众的水平。能把摇滚和交响揉在一起的，恐怕只有崔健，一来崔健以后的中国摇滚缺少集体意识和时代底蕴，更没有作曲、指挥这种正统的音乐教育，二来是人脉。听众的注意力在崔健身上。我的注意力则在摇滚交响，所以我认为大可不必站起来，像听崔健个人演唱会一样――我最佩服的是，崔健和指挥背对背，怎么配合的那么好。
在《蓝色骨头》中，崔健重申了它那三条腿儿的理论，不是按优先级排序，而是"三角架有三条腿才稳定，少了任何一条都要不停的运动"。我便藉此为纲，总结并展望一下吧
"第一就是事业象我上面说的，能高高兴兴工作挣钱养活自己，有话就说有话就写而且要彻底，因为每次彻底之后才会出现美妙的空虚"
去年，计划在技术上深入学习几样东西（Linux Shell/C++/Design Pattern）并没有落实，倒是因为工作需要，启动了Monte Carlo方法的学习。不过后来也没用到，倒成了兴趣。饭否回归，是一件着实令我兴奋的事情，年底开始接触Android，拿饭否API练手，没有目标，只为了在过程中多学点东西。相反，发现工作中并不需要很Fancy的技术，公司总是结果导向。所以明年的技术重点在积累和巩固。拓展的项目仅保留2个：Monte Carlo, Android.
"第二就是身体一定要健康，因为身体要是不舒服，什么都是白给，所以我一周三次跑步加上一次游泳，在运动中想事儿是越想越起劲儿"
难怪崔叔叔都半张儿了，还能在台上那么使劲叫唤那么使劲蹦�，就仗着身子骨好。我要是他那岁数还能一口气唱下来《时代的晚上》间奏的那16个小节，就谢天谢地了。
今年本来说要增重，结果又TM瘦了。该死。明年必须把体重弄上去，不择手段！！！！年底开始了每周招呼同事一起去游泳，还不错，明年继续。
"第三当然就是一个爱情了。。。是不是我的工作太多了感情也变坏了，还是身体一独立个欲望就变野了"
想起有一次去科技园见朋友，聊天，本来答应一起吃晚饭，结果我临下楼接了家里一个电话，就赶紧走了。朋友表示理解：后院不能起火，要走就走吧。我问，你们老这么加班，家里人不急么？他说，时间长了就习惯了。我想，可能我这年头还短吧。有人说，不能总是一方迁就另一方，时间长了不行。也有人说，婚姻就是个错，长久的婚姻就是将错就错。不管怎么说，就算为了自己身体着想，我也要到点儿下班回家，到点儿上床睡觉。谁TM叫我加班熬夜我跟谁急！

评估在线广告的效果

2010-12-29T10:54:00.001+08:00

很久以前看到这样一篇论文
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，是否能假设等同于全体用户的集体偏好？这恐怕是一个难以检验的假设。

wirte C++ for MongoDB: smoking on Ubuntu

2010-12-23T23:03:00.001+08:00

这周利用晚上的时间在笔记本的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了，感觉远不像想象中那样轻量级

Proof of my CNF converter

2010-12-17T20:35:00.001+08:00

定义：
一个句子是一个变量或一个满足下列条件的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

convert to CNF/DNF from any simple sentence

2010-12-16T21:22:00.001+08:00

等着看《让子弹飞》，激动得写了个小程序，把任何不含“非”的简单句子转换成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-14T21:28:00.001+08:00

命题：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-12T20:35:00.001+08:00

在实验过程中发现，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-10T23:25:00.001+08:00

积和式近似

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

积和式定义
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
可以明显看出方差的减小

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

2010-12-01T22:40:00.001+08:00

命题：
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-10T22:00:00.001+08:00

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

序贯蒙特卡洛（一）

2010-11-06T21:53:00.001+08:00

遗传学的例子

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

考虑一组等位基因的演化历史。我们能观察到的，只有今天这组等位基因在某种群的分布。我们提出一个单倍体演化模型，只考虑繁殖和突变，一步概率转移矩阵由某个待定系数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-02T20:09:00.001+08:00

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

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

2010-10-25T20:37:00.001+08:00

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

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

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

没有光驱？死不了的

2010-10-23T22:49:00.001+08:00

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

基本的蒙特卡洛方法（三）

2010-10-14T23:31:00.001+08:00

基本的蒙特卡洛方法（三）

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

分层采样（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-07T16:08:00.001+08:00

多元时间序列分析一个特别醒目的概念，就是协稳：若干个时间序列，他们可能每个都不平稳，但是他们的某个线性组合是平稳的。
在计算广告学实践中，我们考虑一个网站的若干视频，每个视频的日访问量序列看作一个时间序列。很少有哪个视频的访问量是平稳序列，因为人们总是喜新厌旧，一个视频刚上线的一段时间内，看的人越来越多，一段时间后，就渐渐地没什么人看了。但是在成熟、稳定的市场（而非快速增长的市场），一个网站的视频日访问总量在较长一段时间内是平稳的，因为随着老视频陆续淡出，不断有新的视频上线，形成协稳的格局。
当然这只是美好的假设。现在让我们来看看这个假设究竟有多强。
假设有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是必要的。

修改App Engine上的cronjob

2010-10-03T21:50:00.001+08:00

以前用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

How to use VPN in Ubuntu 10.4

2010-10-02T21:15:00.001+08:00

参考
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-09-28T23:16:00.001+08:00

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

假设你开办了一个类似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]，代价是无法预知采纳的样本的个数。

Start to review a textbook with practice

2010-09-24T22:37:00.001+08:00

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-09-14T23:12:00.001+08:00

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

给定集合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”。再由传递性，c≤a。假设a≤c ，则由反对称性，得a=c。a≤b不成立，即c≤b不成立，而b≤c也不成立，与完全性矛盾。

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

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

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

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

URL size limit of Google Chart API

2010-09-13T22:43:00.001+08:00

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-08-18T13:26:00.001+08:00

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-07-27T20:45:00.001+08:00

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-05-08T16:26:00.001+08:00

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

GSL + MinGW + Eclipse

2010-03-14T22:34:00.000+08:00

想当年在学校的时候只用过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);

    }

[Reading] Optimal Online Assignment with Forecasts

2010-03-06T10:14:00.001+08:00

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

Space rail

2010-02-28T15:04:00.001+08:00

听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-02-19T11:27:00.001+08:00

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

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

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

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

Happy Chinese New Year!

2010-02-12T21:29:00.001+08:00

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

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

The limitation of the sampling algorithm: a simple case

2010-02-10T22:24:00.001+08:00

一个不能用采样算法的问题的例子
网站某天有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的无偏、相合估计？

The impact of resampling to RSE in sampling algorithms

2010-02-07T13:37:00.001+08:00

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

在古典概型中有放回和无放回差别很大。在随机采样算法中，对已经采过的样本，是否还可以再采？
问题：
假设有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了：）

Share Two Google Calendar App

2010-01-26T22:05:00.001+08:00

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 :)

Ubuntu 9.10 remember last boot choice

2010-01-24T11:37:00.001+08:00

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

Reading about Refactoring

2010-01-17T21:14:00.000+08:00

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 Resolution

2010-01-05T23:22:00.001+08:00

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

Free SMS

2009-12-27T12:15:00.002+08:00

功能：
通过互联网免费发送短信

原理：
Google Calendar可以通过运营商发送短信提醒
Google App Engine可以当虚拟主机运行简单的程序
Google App Engine上可以跑Google Calendar API的类库

详解：
假设你想给A发短信，而且不想花钱。
1、让A注册一个Google账号，如果A已经有Google账号，此步可以省略；
2、你自己用一个Google账号B，在Calendar里面创建一个新日历，比如名字叫aaa，分享给A；
3、让A打开他账号下的Calendar，如果还没设置SMS，要先点Set up your mobile phone to receive
notifications，选中国、输入手机号，点"发送验证码"，然后A的手机短信收到一个验证码，填到网页的验证码输入框中，点"Finishe
Setup"，然后页面显示状态变成了绿色的"Phone number successfully
validated"。关键--让A找到你分享给A的日历aaa，设置提醒，默认通过SMS，提前事件一分钟发送；
4、你去开通Google App Engine，在xxxxx.appspot.com部署一个网站，网站后台程序自动通过Calendar
API向账号B的日历aaa添加事件，事件的时间设置为当前时间+75秒，这样你基本上就可以实时给A发消息了。程序关键的一段代码(Java)：
private void send_msg(String to, String msg) throws Exception {
// login
CalendarService myService = new CalendarService("useless");
myService.setUserCredentials(USERNAME, PASSWORD);
// get URL
CalendarFeed resultFeed = myService
.getFeed(
new URL(
"http://www.google.com/calendar/feeds/default/owncalendars/full"),
CalendarFeed.class);
String post_url = null;
for (int i = 0; i < resultFeed.getEntries().size(); ++i) {
CalendarEntry entry = resultFeed.getEntries().get(i);
if (entry.getTitle().getPlainText().equals(to)) {
post_url = "http://www.google.com/calendar/feeds"
+ entry.getId().substring(
entry.getId().lastIndexOf("/"))
+ "/private/full";
break;
}
}
// create event
if (post_url == null)
throw new Exception(to + " not found");
CalendarEventEntry myEntry = new CalendarEventEntry();
myEntry.setTitle(new PlainTextConstruct(msg));
long now = DateTime.now().getValue() + (TIMEZONE * 3600 + DELAY) * 1000;
When eventTimes = new When();
eventTimes.setStartTime(new DateTime(now));
eventTimes.setEndTime(new DateTime(now + 60 * 1000));
myEntry.addTime(eventTimes);
myService.insert(new URL(post_url), myEntry);
}
然后外面简单包一层就可以用了
public class SimpleServlet extends HttpServlet {
private static final String USERNAME = "B@gmail.com";
private static final String PASSWORD = "B";
// now() + delay = the event time
private static final long DELAY = 60 + 15;//seconds
private static final int TIMEZONE = 8;//hours
public void doGet(HttpServletRequest req, HttpServletResponse resp)
throws IOException {
String to = req.getParameter("to");
String msg = req.getParameter("msg");
send_msg(to, msg);
}
}
（其中时区的问题你得自己动脑筋想一想）
配置web.xml
<servlet>
<servlet-name>api</servlet-name>
<servlet-class>SimpleServlet</servlet-class>
</servlet>
<servlet-mapping>
<servlet-name>api</servlet-name>
<url-pattern>/api</url-pattern>
</servlet-mapping>
5、通过Web发送短信
打开Firefox或UCWeb，输入网址，形如http://xxx.appspot.com/api?to=aaa&msg=balabala

限制：
msg不能超过57个字符（不同账号可能会有出入）
可能有一天app engine会撞墙
可能有一天你app engine的quota会用完
可能有一天Google会Ban掉你这样调用API的

开发：
个人认为用Java+Eclipse+Google Plug-in最方便，可以本地debug、一键部署
http://code.google.com/appengine/docs/java/tools/eclipse.html

应用：
1、私聊
用Web发送，用短信接收，所以只要能上网，就可以随时随地免费发短信，而且接受者可以不用上网。比飞信省流量，还不用装客户端。P.S.因为你可以登陆账号B，所以如果什么话说出去后又后悔了，可以在15秒内登陆Google
Calendar删除该Event即可。

2、当Twitter用
把自己的日历publish出去，让你那几个闺蜜来订阅。然后你api?to=yourself&msg=balabalalalala她们就收到短信了。

3、群聊
为了识别说话人的身份，Servlet程序需要多做一点工作。不然群里你一言我一语都不知道哪句话是谁说的就乱套了。
另外，怎样才能用一种简单的方式，不让自己的留言也发条短信给自己，我还没想好。

4、Monitor
比如你有一个系统，后台跑一个cronjob，让它
50 0 * * * wget "http://xxxx.appspot.com/api?to=your_team&msg=`grep
ERROR logfile | wc -l`" -O /dev/null

时钟纸牌

2009-10-17T21:52:00.001+08:00

《随机算法》一书中，时钟纸牌是这样一个程序：考虑13*4张标准扑克牌，随机排成13堆。先从第13堆中取一张，然后迭代地，每次从第x堆中取一张牌，其中x是上一次取的牌的面值，直到某次第x堆中没有牌为止。
命题：最后一次取的那张牌是13。
证明：假设最后一次取的牌不是13，是X，那么此时第X堆中没有牌，而程序开始时第X堆中有4张牌，所以在此之前已经从第X堆中取4张牌了，而每次从X堆中取牌，都必须要求上一次翻开的是X，所以在此之前已经翻开4张X了。一共就4张X，这次又翻出一张X，矛盾。

Principle of Shrink Method in collaborative filtering

2009-05-11T16:52:00.001+08:00

在Scalable Collaborative Filtering with Jointly Derived Neighborhood InterpolationWeights等很多篇流行的协同过滤算法的文献中，都提到了Shrink技巧。

Shrink技巧的理论基础是高斯假设下的贝叶斯参数估计，这部分内容是经典的，可以在贝叶斯方法、模式分类等的领域相关教科书中找到。我们将假设和推导过程整理、重述如下：

假设我们要估计参数，视为随机变量，关于的先验知识是它服从正态分布，即其概率密度函数为：

对的一组观测值视为独立同分布的样本，其总体为随机变量，满足条件正态分布。我们关心的是，有了一组观测后，我们应该如何用的信息修正先验知识。为此，计算后验概率密度函数：

其中c为归一化系数，倒数第二步将其整理为的二次型，是为了将其凑成正态密度函数的形式。上式告诉我们服从正态分布，记为，则

将上面两个式子中二项式对应项系数相等，得到下面两个方程：

从中很容易反解出我们关心的：

从中可以得到两点结论：

1、后验期望总可以表示成样本均值和先验期望的凸组合；

2、当观测足够多，即n充分大时，趋向于0，即我们确信非常接近其期望。而收敛到样本均值。

在

中取

记

即得到文献中的Shrink技巧的实用形式：

Implement RBM for netflix CF

2009-05-03T20:13:00.001+08:00

RBM的实现

RBM预测器（隐层节点数F）
参数	权值矩阵W，维数M×K×F(K=5) Bias矩阵b，维数M×K 原始显层节点向量v₀，元素取值于K，维数不超过M，与当前用户的评分个数相同，重构后的显层节点向量v₁，与v₀维数相同，m为v₀下标集到M的映射，即v_0j表示电影m_j的实际评分隐层节点0-1向量第一次采样h₀，第二次采样h₁，维数F 每个用户对应的隐层节点向量的期望Eh，维数U×F
训练算法	1)初始化W~N(0, 0.001),b_ik=电影i得k分的评分次数占电影i评分次数的比例迭代2)和3)直到收敛 2)对于每个用户u： 2.a)根据u的实际评分生成v₀ 2.b)h₀采样，生成随机数p来自U[0,1],若则h_0f=1，否则h_0f =0 2.c)轮盘赌法重构v₁ ：生成随机数p来自U[0,) 令其中k使得 2.d) h₁采样，生成随机数p来自U[0,1],若则h_1f=1，否则h_1f =0 2.e)更新W 若h_0f=1则若h_1f=1则 3)计算Eh
预测函数

实验结果

F=16,lrate=0.002
[1]1.0528 0.0508 1
[2]0.9617 -0.0018 135 0.0911
[3]0.9516 0.0024 273 0.0101
[4]0.9476 0.0007 411 0.0041
[5]0.9454 0.0008 550 0.0022
[6]0.9440 0.0011 689 0.0014
[7]0.9434 -0.0007 828 0.0006
...
[20]0.9404 -0.0003 2636 0.0002
...
[30]0.9400 -0.0024 4029 0.0001
...

各列含义：[迭代步数]RMSE(on ProbeSet) 平均偏差(on ProbeSet) 逝去时间(秒) RMSE改进量
发现前几步RMSE下降非常快。
注：使用Weight Decay可以略微改进结果
将2.e)改成

实验结果
F=16,lrate=0.002,decay=0.0001
...
[20]0.9393 -0.0011 2719 0.0005
...

The Stable Marriage Problem

2009-04-16T10:50:00.001+08:00

稳定婚姻问题

这是从社会生活中过分简化、抽象出来的一个组合数学问题。参考：
http://www.cs.columbia.edu/~evs/intro/stable/writeup.html
http://zh.wikipedia.org/wiki/稳定婚姻问题

严格描述如下：
定义：

考虑n个男人，用大写字母表示(X,Y,Z...)、n个女人，用小写字母表示(x,y,z...)
每个人有一个异性的偏好列表，例如X(x>y>z...)表示X喜欢x胜过喜欢y，喜欢y胜过喜欢z...（静态严格全序，即不存在相等或不可比较的情况，可传递，且不随时间改变）
一个婚姻表示为一个序对，例如Xx表示X与x结婚
一个婚姻组合为所有男人到所有女人的一一对应，表示为n个婚姻，恰好没有人单身
称一对男女(X, y)是不满的，若存在Y、x使得Xx、Yy、X(y>x)、y(X>Y)成立，也就是说，X和x结婚，Y和y结婚，但是X和y喜欢彼此胜过喜欢他们现有的配偶，所以X和y不满于现在的婚姻，存在重新组合的不稳定因素
称一个婚姻组合是稳定的，若其中不包含不满的男女对
一个稳定婚姻问题是指：给定n男n女和他们的偏好列表，问是否存在一个稳定的婚姻组合

定理：

任何一个稳定婚姻问题的实例，都是"是"实例。进一步地，存在一个多项式时间算法，求出稳定婚姻问题的解。

算法（被称作求婚算法）：
初始时所有男人和女人设为单身
每一步

任选一个当前单身的男人X
X向其列表中尚未标记为"拒绝"的、最喜欢的女人x求婚（其存在性在命题0中证明）
如果x单身，则Xx结婚，X从单身男人集合中移除
否则（存在Y，使得已经有Yx）如果x(X>Y)，则Xx，Y放回单身男人集合并从Y的列表中标记x为"拒绝"
否则X留在单身男人集合中并从X的列表中标记x为"拒绝"

直到单身男人集合为空集（可终止性在命题1中证明）

命题0：上述算法每一步的男人X，都能找到可以求婚的x
证明：反证法，假设不然，也就是说这时候X的列表中所有女人都已经拒绝X了。注意到按照算法的构造，只有结过婚的女人才能拒绝男人，而且结过婚的女人只能改嫁不能返回单身，因此这时候所有女人都非单身，从而所有男人都已经结婚了，这与X还单身矛盾。假设不成立，原命题得证。

命题1：上述算法至多n+n^2步终止
证明：每一步要么有一个女人第一次结婚，要么有一个男人被拒绝，注意到n个男人的列表总长度为n^n，因此算法至多n+n^2步终止

命题2：算法的正确性：上述算法终止时，得到一个稳定的婚姻组合
证明：反证法，假设不然，算法终止时存在(X,y)不满，即Xx、Yy、X(y>x)、y(X>Y)
由于X(y>x)，因此X向y求婚，发生在X向x求婚之前，而最终Xx，说明存在一个时刻，y拒绝或抛弃了X
若当时y拒绝了X，说明y的配偶Z在y心中比X好，即y(Z>X)，而最终Yy，说明要么Y就是Z，要么y(Y>Z)，从而都有y(Y>X)，与假设y(X>Y)矛盾
若当时y抛弃了X，说明那时候有一个Z向y求婚，而且y(Z>X)，而最终Yy，说明要么Y就是Z，要么y(Y>Z)，从而都有y(Y>X)，与假设y(X>Y)矛盾
得证。

应用：
这个算法北应用到大学联合招生、

Technical Details in RBM for CF

2009-04-03T15:08:00.001+08:00

Technical Details in RBM for CF
受限玻尔兹曼机的技术细节

参考
[1]http://www.netflixprize.com/community/viewtopic.php?id=1007
[2]http://todwang.blogspot.com/2009/03/understanding-rbm-for-cf.html
[3]http://www.iro.umontreal.ca/~lisa/twiki/bin/view.cgi/Public/DBNPseudoCode
[4]http://imonad.com/blog/2008/10/restricted-boltzmann-machine/
[5]http://www.ics.uci.edu/~welling/teaching/ICS273AFall06/NeuralNets273AFall06.ppt

Gibbs一步采样

在[2]中，我提到训练中要计算<v_ikh_f>_data和<v_ikh_f>_T，而按照[1]中的讨论，T通常取1，因为实践中发现T取大了也得不到更精确的解。因此，我们下面的讨论都是Gibbs一步采样。
这一步如何采，也有两种方式。[3]中<v_ikh_f>₁项中的hf使用了它的期望而不是采样值，而标准的Gibbs采样步骤这步也应该采样。[1]中的讨论中说，[3]的做法精度更好。而我实验发现，这两种方法效果差不多。[3]的具体做法在下一节中一起详细说明。

在线训练

Hinton那篇开创RBM for CF的文章中，用了Batch update weights的方法，也就是我[2]中解释的那样。而[1]讨论中发现，per-case-update（在线训练）也取得很好的效果。具体来说：
1、设trainning case v是一个M×K的0-1阵，表示一个用户的所有评分行为，vik=1表示该用户给电影i评分为k。
2、按照[2]中假设2采样得到隐层的F维0-1向量h，表示用户的特征。
3、按照[2]中假设1采样得到显层的M×K的0-1阵，记作v2，表示重构后的用户评分行为。
4、再按照[2]中假设2计算h的期望Eh，也是F维的向量，元素值取于区间[0,1]。
5、针对这一个trainning case v就可以更新W了：
W[i][k][f]+=LRATE * (v[i][k]*h[f]-v2[i][k]*Eh[f]),foreach i,k,f

利用稀疏性

在我的程序中，v使用了一种相对紧凑的数据结构，注意到v的列和总为1（如果此列不缺失的话）。
ur是一个二元组(mid, rate)的列表，长度和v的非缺失列个数sz相同，rate取值于K，它和v等价：v[ur[j].mid][ur[j].rate]=1恒成立，且v[i][k]=1当且仅当存在j<sz使得ur[j].mid=i,ur[j].rate=k。对应v2为结构ur2。
用这个结构改善上述第5步：
for j in sz
note u as ur[j], u2 as ur2[j]
foreach f
W[u.mid][u.rate][f] += LRATE * h[f];
W[u2.mid][u2.rate][f]-= LRATE * Eh[f];
这样一来，这一步的复杂度从MKF降到sz×F<<M×F了。

Bias设置与训练

[1]中讨论说hidden unit bias都设为零，也就是[2]中d_f设为0，visibel unit bias固定为先验概率分布的对数：b_ik=log(给电影i评k分的用户个数/给电影i有评分的用户数)。b_ik也可以一起在线训练：
bik+=<vik>data-<vik>T
利用稀疏性表示为：
for j in sz
note u as ur[j], u2 as ur2[j]
b[u.mid][u.rate]+=LRATE
b[u2.mid][u2.rate]-=LRATE
我实践发现，visibel unit bias固定或训练效果差不多。

Weight Decay和Momentum

这是常见的梯度下降技巧，[5]中有详细的解释。
Weight Decay可以理解为一个弹簧，把每个weight拉回0点，防止权变太大：

其中lambda是充分小的正数。
W[u.mid][u.rate][f] += LRATE * (h[f]-WEIGHTDECAY*W[u.mid][u.rate][f]);
Momentum是动量项，用于加速梯度下降，有点自适应步长的意思。

其中γ为小于1的正数。
相当于这一步计算的权变化量不仅对这一步迭代有贡献，而且对以后都有贡献，不过贡献随着迭代的次数指数衰减罢了。
注意到，这两个技巧主要是针对Batch update weights的，用到在线训练中的效果不好估计。

退火

这也可以用在所有迭代学习算法中，主要随着收敛步长减小的格式。[1]中给出了一个反比例函数的降温格式。
lrate = LRATE/(1 + iter / (3 * U));
其中LRATE是初始步长，iter是用过的trainning case的个数。可见，迭代30轮（pass），步长就已经下降一个数量级了。

连续受限玻尔兹曼机

参考[4]，注意到评分1,...,5原本是全序的，而不是K个无序的label，而发明RBM的初衷是用于模式分类，最后输出是一个离散label。因此引入连续受限玻尔兹曼机，是自然的。也就是显层向量v不是0-1的，而是连续取值的。这就要求更新[2]中的假设了，需要模型建立者重新假设条件概率分布，也就是重新定义整个贝叶斯网络。[4]引用的那篇Continuous RBM很有意思，用扩散网络建模，求解时用了随机微分方程。以后可以考虑做做这个方向。

结果

可能由于是随机算法，发现rmse随着迭代进行的下降过程有抖动，而不是严格单调的。
目前做到了probe rmse 0.9438 ：迭代46步，其算法参数：
100个hidden unit，在线训练、固定bias、没有退火，lrate=0.001，WEIGHTDECAY=0.0001，没有用动量项。

Understanding RBM for CF

2009-03-26T17:30:00.001+08:00

Understanding RBM for CF

理解协同过滤中的受限玻尔兹曼机

引论

协同过滤（Collaborative Filtering，简称CF）推荐系统中，古典的模型有近邻模型、隐性因子模型等。从2007年开始，受限玻尔兹曼机模型在协同过滤中开始流行起来。
参考文献
[1]Restricted Boltzmann Machines for Collaborative Filtering
[2]On Contrastive Divergence Learning
[3]Markov Chain Monte Carlo and Gibbs Sampling
受限玻尔兹曼机（Restricted Boltzmann Machines，简称RBM）是一种借鉴了模拟退火思想的人工神经网络模型，它在计算过程中又借鉴了蒙特卡洛马尔科夫链(MCMC)的思想。这篇Blog是[1]的阅读笔记，为了理解训练的算法，延伸阅读了[2]和[3]。

问题描述

在线系统中，有很多用户(User)给电影(Movie)评分的数据，协同过滤推荐系统利用这些数据，估计给定User对给定Movie的评分。评分值域为{1,...,K}，用户集合记作U={0,...,U-1}，电影集合记作M={0,...,M-1}，已知评分用有缺失值的矩阵R表示，，推荐系统使用某种协同过滤算法评分器估计(u,i)点处的评分记作，我们希望在给定测试集Q上均方根误差尽量小。

模型框架

存在这样一个对称网络可以表示用户给电影评分的机制，它包括显层和隐层两层节点。
显层共M×K个节点，对应M个电影和K个可能评分的卡氏积。每个u∈U的实际评分对应一组显层状态（Visible），一组显层状态用M×K的、允许缺失行的、0-1矩阵V表示，它刻画了u的实际行为。V=(v_ik)，若V的第i行数据缺失，则表示u没有给i评分，否则V的第i行总恰有一个非零元，v_ik=1表示u给i评分为k。
隐层共F个节点，对应F个隐藏的因素（Factor）。隐层表示为F维0-1向量h，非零元代表u具备的因素，它描述了u行为的原因。
任一显层节点(i,k)和任一隐层节点f之间用无向弧相连，弧的权值记作，表示两个节点联系的程度。权值矩阵W为M行K列F层的三维矩阵。不同的u，对应不同的V和h，但是他们共用一个W。根据显层状态分布假设P(V|h)，用隐层状态h和权值W决定显层状态的过程，是预测的过程。利用显层状态V反推权值W和隐层状态的过程，是训练的过程。作为Hopfiled网络，学习过程是以下两个步骤的交替：
1、激活：固定W，根据隐层状态分布假设P(h|V)，计算h的分布
2、反向传播误差修正W：固定h，极大化似然函数P(V)，计算它关于W各分量的偏导数，沿负梯度方向改变一充分小步长来修正W。

模型假设

固定一个u，
1、显层状态的条件概率分布：

其中

这里a_ik表示u给电影i的评k分的激发程度，它关于h和W_ik的向量内积成指数增加，也就是说，i具备u评k分的因素越多，被u评k分的可能性越大，而且这种正相关的影响呈指数增加。b_ik为电影i的得k分的内在因素，与u无关。
2、隐层状态的条件概率分布：

其中W^f为W的第f层，

d_f为因素f被观测数据影响的因子，它只与因素f自身有关，与u和i无关，
sigmoid(x)=1/(1+e^-x)
注意到，h_f=1的条件概率，数值上就是h_f的条件期望。
3、似然函数

其中c为归一化因子，E(V,h)为能量函数

疑问：
这个假设如何从前两个假设推出来？

训练

计算似然函数偏导数

其中
为训练过程遍历所有u∈U时，v_ikh_f非零的频率，具体来说，其中h(u)根据假设2采样得到；
是v_ikh_f非零的频数的数学期望，它的值由模型假定的分布决定，但计算复杂度太高，按照[2]介绍的Contrastive Divergence方法用近似；
通过Gibbs采样[3]迭代T次得到，T初始为1，随着训练T的值不断增加，[1]中说T太大也没必要。Gibbs采样是一种蒙特卡洛方法，意思有点像“左右互搏”，在这个问题中，就是先固定V0=V，根据假设2采样得到h1,然后根据假设1采样得到V1，这叫迭代了一步，迭代T步求平均，就可以近似它的期望了。具体来说
有了似然函数个方向偏导数，就可以“爬山”啦~

lrate为学习率，或者叫步长。

预测

训练好了权值矩阵后，就可以用它来预测(u,i)∈U×M了。给定u，先根据假设2算出h的条件期望，然后用它来代替假设1中的h，算出v_i1,...,v_iK的分布，然后计算r_ui期望作为预测值。具体来说：
step1

step2

step3

Netflix数据集的一些统计特征(2)

2009-03-19T23:06:00.001+08:00

很多网站的流量在一周之内变化的趋势类似，那就是工作日流量多，周末流量少。Netflix数据集也如此。无论是训练集还是测试集，评分日期的一周分布中，周六周日两天加起来不到20%，而周一周二尤其多，加起来超过了35%，说明人们上班的时候很闲可以上网，或者晚上喜欢放松看看电影，而周末则外出休闲去鸟。

另一个发现，是训练集和测试集的某些边缘分布很不同（测试集包括probe set和quiz set，分布相同）。
1、训练集的用户分布不均匀，少数用户非常活跃，10%最活跃的用户评分点占了所有测试集评分点的43.6%，而测试集的用户分布是均匀的，几乎所有用户被测试的概率相同，除了极少数最不活跃的用户。说明Netflix对所有用户一视同仁，不活跃的用户被认为与活跃的用户有相同的消费的能力。
2、测试集比训练集的时间更集中于离现在较近的时刻，测试集中最后一周的评分点占了20.4%，而训练集中最后一周的评分点仅占了1.5%。说明要求模型有一定预测的能力，尤其是预测当下用户偏好，而不是猜测用户过去偏好的能力，因为前者是有价值的。

run mysql with ruby in ubuntu 8.04

2009-03-04T10:36:00.001+08:00

用新立得安装ruby,mysql,mysql server的最新版
下载mysql-ruby-2.8.1，http://www.tmtm.org/en/mysql/ruby/
按照上面网址的说明，配置、编译、安装

中途遇到几个小麻烦，mkmf包在ubuntu 8.04默认的ruby 1.8中没有，mysql客户端的版本也有些问题。总之全都使用最新版就OK了。

Second-hand Hash

2009-02-25T17:23:00.001+08:00

二手散列这个名字是我起的，官方叫法是Cuckoo hashing，他是一个散列算法框架而不是一个散列函数。事实上，二手散列需要两个散列函数h1和h2，对应两个散列表T1和T2，散列函数从值域（domain of values）映射到键域（domain of keys），散列表把键应成值。判断值x是否存在于字典D中：lookup(x,D)=(T1[h1(x)]==x || T2[h2(x)]==x)。
插入的方式很有趣，如果有散列冲突，就把原有值的踢出去，提到另一个表里去，知道没有冲突为止。用伪代码来说：
insert x into D:
如果lookup(x,D)返回真，返回
重复maxloop次：
    若T1[h1(x)]为空，则T1[h1(x)]=x后返回
    交换T1[h1(x)]和x
    若T2[h2(x)]为空，则T2[h2(x)]=x后返回
    交换T2[h2(x)]和x
如果代码走到了这一行，调整T1,T2的大小和h1，h2的值域，重新insert x into D

这是一个完美散列，因为我们确定地知道，lookup至多探查两次，绝对O(1)。难点在insert的复杂度的分析。注意到在insert中使用了递归调用，因此必需首先证明，insert的做法能够以较大的概率在有限步内终止。然后还能证明，insert的均摊复杂度的期望也是常数。证明过程这里就不转载了。
------
另外有一种d手链式散列的框架，针对一种特殊的问题，也有很漂亮的结果，出自《概率与计算》一书。n个球放入n个箱子，每个球放的时候先从n个箱子中等概率地选d(>=2)个箱子，然后放入这d个箱子中球最少的那个箱子。可以证明，以1-o(1)的概率，球最多的箱子中球的个数（称为最大负荷）是
基于此，我们可以开发出确定性O(d)时间插入、依很大的概率在至多lnlnn/lnd时间查找的链式散列算法。
d=1时，没有选择的余地，书中引理5.12告诉我们：n个球独立等概率放入n个箱子，最大负荷以1-1/n的概率至少为lnn/lnlnn。这告诉我们，在散列表比较满的时候，“一手”链式散列的性能和二分查找相当。对比发现，d手链式散列对传统链式散列的改进还是挺大的。
------
应用中，二手散列或d手链式散列还有一个好处：在多核处理器上可以并行计算散列函数，进一步提高效率。
------
P.S.复杂度的概率分析中大量用到了切尔诺夫界技巧

Export all SMS messages as xml file from your Windows Mobile

2009-02-23T12:03:00.001+08:00

昨天觉得手机中病毒了，老莫名其妙拨号上网。准备重装系统，先备份短信。.NET Compact Framework类库中outlook类中的SMS若干方法形同虚设。于是找了一个第三方.NET的类库，MAPIdotnet。他构建在Windows API的MAPIlib基础上，专门用来访问Windows Mobile的短信的。他本身是for Smartphone的，稍微改一下就能在我的Pocket PC上跑了。
            //export all
            FileStream fs = new FileStream("smsExport"+DateTime.Now.ToString("yy_MM_d_H_mm")+".xml", FileMode.OpenOrCreate);
            StreamWriter sw = new StreamWriter(fs, Encoding.UTF8);
            //
            sw.Write("<?xml version=\"1.0\" encoding=\"UTF-8\" ?><sms>");
            for (int k = 0, K = this.stores.Length; k < K; k++)
            {
                IMAPIMsgStore store = this.stores[k];
                sw.Write("<store name=\""+store.ToString()+"\">");
                IMAPIFolder folder = store.RootFolder.OpenFolder();
                IMAPIFolderID[] subFolders = folder.GetSubFolders((int)folder.NumSubFolders);
                foreach (IMAPIFolderID fId in subFolders)
                {
                    IMAPIFolder f = fId.OpenFolder();
                    sw.Write("<folder name=\"" + f.ToString() + "\">");
                    IMAPIMessage[] messages = f.GetNextMessages(f.NumSubItems);
                    for (int i = 0, length = messages.Length; i < length; i++)
                    {
                        IMAPIMessage msg = messages[i];
                        msg.PopulateProperties(EMessageProperties.DeliveryTime | EMessageProperties.Sender | EMessageProperties.Subject);
                        //store folder subfolder msg
                        sw.Write("<msg when=\"" + msg.LocalDeliveryTime.ToString("H:mm d/MM/yy") + "\">");
                        sw.Write("<Sender name=\"" + msg.Sender.Name + "\"><![CDATA[" + msg.Sender.FullAddress + "]]></Sender>");
                        sw.Write("<Recipients>");
                        for (int j = 0; j < msg.Recipients.Length; j++)
                        {
                            IMAPIContact recipient = msg.Recipients[j];
                            sw.Write("<Recipient name=\"" + recipient.Name + "\">" + recipient.FullAddress + "</Recipient>");
                        }
                        sw.Write("</Recipients>");
                        sw.Write("<Subject><![CDATA["+msg.Subject+"]]></Subject>");
                     //   sw.Write("<Body><![CDATA[" + msg.Body + "]]></Body>");
                        sw.Write("</msg>");
                    }
                    sw.Write("</folder>");
                }
                sw.Write("</store>");
            }
            sw.Write("</sms>");
            //
            sw.Close();
            fs.Close();

Chernoff bound

2009-02-21T19:04:00.001+08:00

Chernoff bound

《概率与计算》一书中介绍了一个叫作“切尔诺夫界”的不等式，在概率论中，可以用来估计随机变量尾分布的性质，也就是随机变量在远离其期望值的某区域中取值的概率。切尔诺夫界可以用于分析算法的概率意义下的性能，具体来说，能给出算法运行时间超过某一阈值的概率的上界，或者解的相对误差超过某一阈值的概率的上界。

Wikipedia上有切尔诺夫界的若干中形式的公式。我在这里不加证明地引用其中一个比较简单的公式，然后给出一个它在算法分析上应用。

切尔懦夫定理：X₁,...,X_n是独立随机变量，X_i服从0-1分布，，则对任意的ε，0<ε<1，有：。

考虑这样一个判定问题：G=(V,E)为无向图，D={(s_i,t_i),i=1,...,n,}为一个顶点对的集合，问是否存在{p_i,i=1,...,n}，其中p_i为连接s_i到t_i的路，使得任何一条边e属于E都没有被两条路p_i、p_j覆盖。
对应解一个优化问题：给一组路{p_i,i=1,...,n}，对于任一e属于E，定义cong(e)=|{p_i|e属于p_i,i=1,...,n}|，即通过e的路的条数。求{p_i,i=1,...,n}使得C=max{cong(e),e属于E}达到最小，即让被通过次数最多的边，通过次数最少。

这个问题是NP难的，不过我们可以设计这样一个算法：线性规划松弛，然后随机近似。通过切尔诺夫界能够给出这个算法的近似率的一个估计。

令P_i为从s_i到t_i的所有路的集合，记f(p,i)=1{p属于P_i}，也就是表示p是否是一条从s_i到t_i的路。优化目标是min C。约束条件有下面两个：
首先任何一对点s_i,t_i恰有一条路p通过，
，对于任意i，
其次每一条边的冲突数不超过C，
，对任意e属于E。
我们把f(p,i)松弛到[0,1]，用线性规划的方式求解。然后对每个i=1,...,n，用轮盘赌的方式以概率f(p_i,i)选p_i属于Pi。

上面就是线性规划松弛加随机近似的算法，下面用概率的方法分析这个算法。

定义随机变量表示路选中的路p_i是否通过边e，i=1,...,n，e属于E。则：

定义，表示边e的被通过的次数。则：

由切尔诺夫界定理，立得：

旁听《电视新闻研究》1

2009-02-17T23:45:00.001+08:00

从这周开始，在北大旁听《电视新闻研究》，阿忆讲的。其实早就对电视媒体发生了兴趣，尤其是互联网对电视媒体的变革。要想读懂互联网电视，我认为首先应该读懂电视，然后才是读懂互联网。新闻，作为内容的一种，既出现于传统平面媒体、电视媒体，也出现于互联网后的新媒体。新闻，是媒体之所以存在的最原始也是最具有生命力的理由。
这个系列Blog是随记，不仅仅是课堂笔记，还记录了这门课上其他所见、所想。

2009年2月17日
"今天下午，我和丽冒着2009年的第一场雪骑车去北大找硕硕旁听阿忆的《电视新闻研究》。"
按照新闻学的分类，上面这句话是不折不扣的"硬新闻"。今天这堂课，主要就讲了一个"软新闻"的概念，它与"硬新闻"相对，有很多版本的定义，我理解的软新闻，就是讲故事，而硬新闻，则是中国人原来印象中的新闻。
"软新闻"不是电视媒体发明的，而是平面杂志发明的。1968年以前，电视新闻主要是硬新闻，因为它在时效性和感染力方面有着得天独厚的媒介优势，不需要用内容上的革新就能赢得广泛受众。而受到杂志上频见的"软新闻"的启发，CBS的《Cronkite晚间新闻》的执行制片人Don Hewitt在1968年引入了个性化新闻的理念，将"软新闻"带进电视媒体。值得注意的时间节点有两个，一个是1968年，当年《60 minutes》栏目就上线了，说明他的团队执行力强，另一个是在栏目不受好评的低谷中坚持摸索了7年，在1975年才成为相当于中国春晚一样的收视率极高的节目，说明他们能厚积薄发，耐得住寂寞。
这让我想到互联网电视，它相对于广播电视，和广播电视相对于平面媒体有很多相似之处。首先它们都是技术的革新带来媒介上的突破，使得新闻传播的更实时、更鲜活，影响力更大。其次，他们都在抢占传统媒体市场的同时，学习传统媒体中有价值的元素，例如电视从平面学来了"软新闻"，而互联网也从电视学来了"文字直播"、"参与式新闻"等等。新事物的不会凭空产生，他们必须在旧事物中汲取养分，才能成长。互联网电视也一样。这也是我为什么想做互联网电视，却热心于学习传统电视媒体的原因。
《60 minutes》的Hewitt改变了新闻理念。Now, News is the story you never heard。好新闻的标准，据说是"昨晚播出的节目成为第二天早上人们的话题"。不过，我更认可另一种也许更保守、更古典的说法：新闻是历史的第一稿。仔细想想，二者并不矛盾，只是表意倾向不同罢了。
关于电视广告，有这样一个案例，每期《60 minutes》分为3个独立新闻故事，中间插播两个广告。即使这节目中抨击了广告主，广告主很不爽，却还要在这里打广告，因为其他地方的回报无法和这里相提并论。这让我想起以前经常听做互联网电视的朋友提起的一句话：内容为王。

也许是学校不同，也许是专业不同，北大新闻系的同学给我感觉放的更开。上课时，阿忆一度困扰于笔记本电脑接到音箱上没有声音，很快第一排主动站起来一个男生一个女生，走上讲台帮阿忆解决问题。女生打电话，和物业的人员沟通，按照提示尝试讲台上的各种开关，沟通很得体很有效率，男生则和女生配合，帮着拉开讲台的柜门、查看各种开关的说明等，很快找到了正确的开关，解决了问题：能听到阿忆笔记本电脑的声音了。他们也许上过阿忆的课，也许在坐很多都是他俩的同学罢，所以感觉特别自然，就像中学的一个班上，同学帮老师去拿粉笔或帮老师擦黑板一样，给人感觉很亲切很温馨。在清华理工科的课堂上，我从没见过这种场面。遇到笔记本电脑声音放不出来，要么是所有同学都对老师的无助无动于衷或爱莫能助，要么是于无声处跳出来一个惊世骇俗的大牛，上去一按，电脑就出声了。

正心诚意读书不二

2009-02-15T23:56:00.001+08:00

最近经历了一些大悲大喜，发现自己内心越发平静了。
正心诚意，出自《大学》，刚上大学时特别特别喜欢的一句话，烙在了骨子里。虽然多年来学习成绩一直很一般，但是却自以为学得很洒脱，很愿意砸时间去享受自由的思考的乐趣。"古之欲明明德于天下者，先治其国；欲治其国者，先齐其家；欲齐其家者，先修其身；欲修其身者，先正其心；欲正其心者，先诚其意；欲诚其意者，先致其知，致知在格物。"格物，我以为指的是学好物理和数学，或者说了解自然和理性。这一串"责任链"中最重要的，我认为是正心诚意。自古以来，许多治国者没有很好的根基，风光一时甚至一世，但最终被人民群众和历史否定，无法达到"明明德于天下"的儒家理想。也许是我的唯物历史观太教条了罢。。。
读书不二，出自曾文正公的自律条款。"敬、静坐、早起、读书不二、读史、谨言、养气、保身、日知所亡、月无亡不能、作字、夜不出门"。我特别喜欢其中"读书不二"这句，以为是不要同时读两本书的意思。很多时候，我们到图书馆会搞一大堆书，东翻翻西看看结果哪一本都没看透。这样不好。最近我在潜心读一本书《Algorithms for hard problems》，去图书馆只带这一本书、一支笔、一个笔记本（纸的）。只读一本书时，心灵会觉得特别纯净。虽然这本书很难啃，而且对眼下的毕设、对即将面对的工作没有多少直接帮助，却提供了很多我认为很有价值的算法设计与分析的思想，读起来很享受。专注，是一种做学问的态度，也是一种做人的态度，适合一小撮人。他们是畏惧社会吗？看你怎么看了。
正心诚意，读书不二。继续，作我的快乐学术男。

赫鲁晓夫在哪里

2009-02-06T19:01:00.001+08:00

1965年2月14日，苏共第二十次代表大会在莫斯科召开。2月24日深夜，赫鲁晓夫忽然做了《秘密报告》，系统揭露和批评斯大林的重大错误，要求肃清个人崇拜的流毒。顿时在国内外引起强烈反响。很多人就有疑问，你既然知道他的错误，为什么在斯大林生前不提出来呢？后来，在党的代表会上，赫鲁晓夫又谈此问题时，有人从听众席传来一张纸条，上面写着：当时你在哪里？
赫鲁晓夫不愧是政治家，看看他是如何机智应对的。
他拿起纸条，大声念出了上面的内容，然后向台下喊道："请写这张纸条的人，马上从座位上站起来，并走到台上。"台下鸦雀无声。赫鲁晓夫又重复了一遍，台下依然一片死寂。于是，赫鲁晓夫淡淡地说："好吧，就让我告诉你，当时我就坐在你现在坐的那个位置上"。

这个故事来自一本博弈论的科普读物，当时正在讲多人博弈的人质困境。

Simulating N-Body Problem

2009-01-20T10:26:00.001+08:00

理解多体问题的数值模拟

多体问题
物理上，多体问题（N-Body Problem）是指：3维空间中给n个质点在某一时刻的位置和速度，已知他们的质量，假定他们彼此之间只受万有引力作用，求他们后续任意时刻的位置和速度。数学上，多体问题是一个高阶的常微分方程初值问题，具体说来：
设位置函数向量x=(x₁,...,x_n)，，质量m_i
满足方程

其分量形式为

和初值条件

求x(t)
当n=2时（又称二体问题），有解析解，解曲线是平面上的二次曲线（《常微分方程教程》P132）。
n>2时目前没有解析解，虽然理论上得到了很多很漂亮的定性分析结果，但是实践中还要归诸数值模拟，也就是用差分代替微分。

多体问题的数值模拟
Part1：差分格式
参考http://en.wikipedia.org/wiki/Verlet_integration
数值模拟最朴素的想法是欧拉迭代，在每一个时刻，根据当前位置计算加速度，更新下一时刻的速度，根据当前速度更新下一时刻位置：（设时间步长为h）
a(t) = f(x(t))
v(t + h) = v(t) + a(t)h
x(t + h) = x(t) + v(t)h
这个方法误差很大，即使假设速度v没有产生误差，也很容易分析出（通过Taylor展开到二阶项），迭代一步产生的误差为O(h²)。

经典力学中的数值模拟常用的是Verlet方法，它利用两侧Taylor展开消去速度项，构造误差阶数更小的差分格式：
首先Taylor展开到4阶

注意到

并将-h带入上面的h，得

两式相加，消去1阶项，顺带把3阶项也消掉了（赚到了，呵呵）。得到：

定义差分格式为

显然一步误差O(h⁴)
然后来分析误差扩散。
（注：通常Big O Notation是不能做减法的，这里我们不加说明地假设它可以做减法，使用Lagrange余项的Taylor展开也许能解释清楚，我们有深究它。）
来归纳证明，k步误差有如下二次扩散关系

k=1、2时不难直接验证。假设直到k-1都成立，来看k时

其中假设f连续可微，x取值在一个紧集上，就有f是Lipschitz连续的，从而两个h²f项做差是高阶小，直接扔掉。从而上式成立。
因此，若模拟总时间T=kh固定，则总误差O(h²)

Part2：受力计算的分治算法
参考http://www.cs.berkeley.edu/~demmel/cs267/lecture26/lecture26.html
计算每个质点受的合力，按照定义，是其他质点对它的万有引力的矢量和。n个点，每个点计算(n-1)个分力，这样迭代一步的复杂度是O(n²)，不能忍。
Barnes-Hut算法使用分治策略，将复杂度降低到O(nlogn)，它的思想是，将较远的若干个点当做一个点，一起计算来自他们的分力。举例来说，假设x,y到点集A的质心的距离很远，A包含a个点，则按照定义计算x,y受A的力需要计算2a个分力，而若把A当做一个对象处理，预先计算它的质心和总质量，则完成同样的工作，仅需要a次计算。
算法用到了四叉树的结构（这是对于2维情形，3维使用八叉树）。
递归地将矩形分成四等分，落入第i象限的点属于当前节点的第i子节点，直到当前节点只有一个节点为止，依此从平面上的点集构造四叉树，如图——

通过适当选取划分轴，可以保证四叉树均衡，也称自适应四叉树，这一步复杂度O(nlogn)
得到四叉树后，后续遍历一次得到每个节点的质心和总质量。最后，计算每个质点的受力：
从根节点开始考虑每个子树，如果质点到该子树的根节点的质心足够远，则累加质点受到该子树（作为整体对应的质心和总质量）的万有引力，否则，质点受到该子树的合力等于该子树的每个子树的引力矢量和，而计算子树的引力又递归调用当前算法。由四叉树对应的几何空间划分，很容易看出，计算一个质点的受力，要遍历树上的节点个数不超过树的高度。若四叉树均衡，则复杂度O(nlogn)，从而总复杂度也如此。

总结：
1、利用Taylor展开构造有效差分格式可以减少计算误差
2、利用分治思想改进算法可以降低计算复杂度

题外话：
计算数学看算法就看三件事：收敛性、复杂性、稳定性。这个例子中讨论了前两个，第三个稳定性主要是针对良好的问题的，也就是说，如果真实解对初值的依赖是连续的，则要求算法的数值解对初值的依赖也是连续的。但是多体问题是一个典型的混沌系统，问题本身是病态的，讨论它的算法稳定性也就意义不大了。
一直以来，我就特别看好物理的sense好的人，像上面这个分治算法的设计，物理背景就很强（berkeley的讲义上用仙女座星系做四叉树构造思想的导言），或者说，对问题域有透彻了解的人，往往能发明出新的、好的方法。

Keypoints for Basic Functional Analysis

2009-01-13T09:09:00.001+08:00

Keypoint for Basic Functional Analysis
基础泛函分析知识点整理
度量空间

聚点：x∈G'，若x的任何去心邻域交G非空

闭集的等价定义：
1、余集是开集
2、闭包为自身
3、若其中一序列收敛，则极限点属于集合

闭包的等价定义
1、G的闭包=G∪G'
2、G的闭包是包含G的最小闭集

可分空间的子空间可分
对X'可分，则X可分
C[0,1]可分
l∞不可分

完备子集闭
完备空间的闭子空间完备
自反，则完备
离散度量空间总完备
lp、l∞、c0、c、C[a,b]完备

M是紧子集，若M的任何序列都有收敛到M中的子列
紧空间完备
紧子集是有界闭集
紧集中的闭集是紧的
点到紧集的距离总可达到

压缩映射：T:X->X成为压缩的，若存在0<a<1，d(Tx,Ty)<=ad(x,y)对于任意x,y∈X
不动点：x为T不动点，若Tx=x
不动点定理：X完备，存在m>0使得Tm压缩，则T有唯一不动点

连续映射：
开集的原像是开集、闭集的原像是闭集
把紧集映成紧集
与极限号可交换

线性空间

G是Hamel基，若G无关且spanG=X
有限维线性空间代数自反

实线性空间X，Z为子空间，p为X上次线性泛函，f属于Z*，f(x)<=p(x)对于任意x∈Z，则存在g属于X*使得g|Z=f，g(x)<=p(x)对于任意x∈X
线性空间X，Z为子空间，p为X上半范，|f(x)|<=p(x)对于任意x∈Z，则存在g属于X*使得g|Z=f，|g(x)|<=p(x)对于任意x∈X

赋范空间

度量具有平移不变性、比例缩放、球可线性运算

(e_n)为Schauder基，若对于任意x∈X，存在唯一λ_n使得x=∑λ_ne_n
X有Schauder基，则X可分
Schauder基无关，则dimX=∞

Reiz引理：X是赋范空间，Y是Z真子空间，Y是闭集。对于任意0<θ<1，存在z∈Z，||z||=1，使得ρ(z,Y)>θ
X是有限维赋范空间，(e_i)是X的基，则存在a>0，使得a||x||₁<=||x||对于任意x成立
有限维赋范空间是Banach空间
有限维赋范空间中有界闭集是紧集
赋范空间X，X有限维当且仅当单位闭球是紧的
赋范空间中，点到有限维子空间的距离总可达到

X严格凸，点到闭子空间的距离至多可达一处

X,Y赋范空间，T:X->Y线性，则三者等价：T连续，T在0处连续，T有界算子

Z为X子空间，f∈Z'，则存在g∈X'，g|Z=f，|g|=|f|
对于非零元x，存在f∈X'，|f|=1，f(x)=|x|
|x|=max{f(x)|f∈X'，|f|=1}
Y真闭，x=X-Y，则存在f∈X'，|f|=1，f|Y=0，f(x)=ρ(x,Y)>0

T属于B(X,Y)，定义共轭算子T×:Y'->X'，T×f定义为fT，f∈Y'

自反，则完备
H自反，lp自反
l1、l∞、c0、C[a,b]不自反

Banach空间

赋范空间完备，当且仅当任何绝对收敛的级数都收敛
X,Y赋范空间，B(X,Y)完备当且仅当Y完备

X完备，Tn∈B(X,Y)：

对于任意x∈X，Tnx有界，则Tn有界
Tn弱收敛到T，则T有界
T∈B(X,Y)，Tn强收敛到T，当且仅当Tn有界且Tn在X的一个完全集上处处收敛到T

X,Y完备：

T∈B(X,Y)，若T满，则T为开映射，若T双，则T逆连续
D(T)闭，T闭，则T∈B(X,Y)

X,Y赋范空间，T∈B(D(T),Y)：

D(T)闭，则T闭
T闭，Y完备，则D(T)闭

Hilbert空间

范数由内积给出，当且仅当平行四边形等式成立|x+y|²+|x-y|²=2|x|²+2|y|²

内积空间中，点到完备凸集的距离唯一达到，特别到闭子空间对
内积空间中，点到子空间的距离可达到

任何集合的正交补都是闭子空间
X=MN，若任意x∈X，存在唯一m∈M,n∈N，使得x=m+n
X为Hilbert空间，M为闭子空间，则X=MM正交补
M完全当且仅当M正交补为{0}

Bessal不等式：(ei)为标准正交列，则∑|<x,ei>|2<=||x||2

完全标准正交集，可数则Schauder基，有限则Hamel基
Hilbert空间总有完全标准正交基
伴随算子：H1,H2为Hilbert空间，T∈B(H1,H2)，T*∈B(H2,H1)使得<T*y,x>=<y,Tx>

点到集合的距离
有限维，紧，内积，凸

广播加密算法调研报告

2009-01-12T12:08:00.001+08:00

广播加密算法调研报告

一、广播加密问题
1、问题的提法
一个系统有一个广播源和n个用户，订阅者S为用户集的一个子集，大小为k。广播源将消息m加密后广播给所有用户，使得S中的用户解密得到m，非S中的用户无法解密得到m。其中如何加密解密，并使得密钥和密文长度尽量短，就是广播加密问题。相当于多把钥匙一把锁，锁可以制定，想让谁能开谁就能开，而且想让谁能开谁才能开。
注意到，不可以给每个订阅者单独加密，这样虽然可以把广播加密问题归结为已经解决的点对点加密问题，但是实际中不能接受：应为那样广播源加密的计算开销将随着k的增大而线性增大。这是广播加密问题和传统的加密问题的区别。
广播加密问题的形式化提法：
通过Setup(用户总数 n)函数得到一个广播源的公钥PK，每个用户i的私钥di，i=1,..,n。广播源针对订阅者集合S加密消息m，即执行加密函数Encrypt(广播明文m, 订阅者集合S, 公钥 PK)得到密文Hdr，通过网络广播Hdr给所有用户。所有用户都有相同的解密函数Decrypt(密文Hdr,私钥d)，用户i接收到Hdr后计算Decrypt(Hdr, di)，Decrypt(Hdr, di)=m当且仅当i属于S。注意到，初始化公钥和私钥的算法与订阅者S无关，实际中订阅者会经常变化。我们希望开发出其中Setup、Encrypt、Decrypt的算法，满足一系列安全要求，并且PK、di关于n不太长，给定m后Hdr不太长。
关于广播正文m的一点注解：实际中广播内容是一段比较长的视频流stream，但是用此非对称加密算法加密stream可能无法满足实时性需要，因此，我们用计算复杂度比较低的某种对称加密算法加密stream，设对称密钥为K，然后用广播加密算法（非对称加密）加密此对称密钥K，然后广播非对称加密后的K和用K对称加密过的stream即可满足实时性的需要。因此我们下面讨论的广播正文m，实际就是这里的对称密钥K，因此，m只需要是某个随机的整数即可，与待广播的视频流stream无关。

2、广播加密问题的若干安全要求：
防破解：防止非S中的用户通过Hdr得到m。
盗版追踪：如果S中的某个i将di出卖给非S中的j，如何检测到并确定i的行为。
防恶意插播：如果另一个非法广播源在广播信道中广播内容，接收端如何识别，也就是如何验证广播源的身份。

3、接收端无状态
[todo]

二、广播加密算法
1、多项式插值
91年Berkovits的《How to broadcast a secret¹》提出这样一个基于拉格朗日多项式插值的广播加密算法，它利用了这样一个事实：一个n阶多项式能被平面上（横坐标互不相同的）n+1个点确定。
算法构造的思路如下（a special but typical case）：
Setup(用户总数 n)
随机生成平面上横坐标互不相同的n个点，分配(xi,yi)作为用户i的私钥di。不妨设所有xi都不为0。
不用公钥。
Encrypt(广播明文m, 订阅者集合S)，S的大小为k
计算一个k阶多项式P，过下列k+1个点：(0,m), (xi,yi) ,i属于S
然后广播密文Hdr为一个长度为k的串/数组：(xi,yi)，i不属于U，且(xi,yi)在P上，不妨设所有xi都不为0。
Decrypt(密文Hdr,私钥di)
如果用户i属于S，则i知道k+1个在P上的点（k个来自Hdr，还有一个是自己的私钥），从而确定P并能够求出P在0点的值，也就是密文m
否则i只知道k个在P上的点，自己的私钥(xi,yi)不在P上，无法求出P在0点的值。

实际加密的时候还要引入j个扰乱的点，可能是为了安全，为了计算方便，作者给出了矩阵描述的等价算法。

2、子集覆盖
3、配对
05年Boneh等人的《Collusion Resistant Broadcast Encryption With Short Ciphertexts and Private Keys²》利用配对技术构造了这样一个广播加密算法，它利用了双线性DH假设。具体来说，给定误差限，给定p阶双线性群G和它上面的双线性映射e（相关定义和性质见附录），若h,g随机取值于G中，随机取值于中，则不存在这样一个（确定性的？且多项式时间的？）算法A使得
，其中
注意到，如果能在G中从g和计算出，则显然可计算出因为其中已知。在G中从g和计算出就是经典的离散对数计算，因此双线性DH假设比离散对数假设更强，如果离散对数多项式可解则双线性DH假设不再成立！（介于二者难度之间的还有一个DH假设，见附录）
算法构造的思路如下（a special but typical case）:
Setup(用户总数 n)
随机取值于中，令
公钥PK：
私钥di：
Encrypt(广播明文m, 订阅者集合S)
t随机取值于中，不妨设，广播密文为一个二元数组Hdr=(C0,C1)，其中

Decrypt(密文Hdr,私钥di)

三、研究者
1、Naor
2、Boneh
四、应用
1、DRM
2、CA for IPTV
附、Mathematical Background

¹http://www.springerlink.com/index/XXBDHG4QH59RGCH4.pdf

²http://crypto.stanford.edu/~dabo/papers/broadcast.pdf

王元涛的Blog

凯立德刷K币

如果我的证明是正确的，就不需要任何其他形式的肯定

2012年1月1日不是2012年的第一周

幂级数

./

驾照换证攻略

2012 resolution

我家的挂表

新手提车小记

马尔代夫满月岛攻略

测试MathJax

2-SAT

1977

Fwd: 分析笑话：数学系有三个班

Fwd: 内省+变长参数的lambda表达式

From Evernote:

内省+变长参数的lambda表达式

滑雪板理论

上帝存在，帕斯卡说的

“我正在说的这句话是谎话”的三种等价说法

修改注册表防止sftpdrive试用版过期

期权定价

自反+反向传递=>完全

也发美国商务签证面经

2011推荐系统高峰论坛――参后感

A spike

马尔科夫链蒙特卡洛（二）

马尔科夫链蒙特卡洛（二）

马尔科夫链蒙特卡洛（一）

马尔科夫链蒙特卡洛（一）

遍历trie

两个不幸等价的随机数生成器

给我碗儿水喝

评估在线广告的效果

wirte C++ for MongoDB: smoking on Ubuntu

Proof of my CNF converter

convert to CNF/DNF from any simple sentence

序贯蒙特卡洛（三）

序贯蒙特卡洛（二）代码

序贯蒙特卡洛（二）

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

如果饭否回来了

序贯蒙特卡洛（一）

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

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

没有光驱？死不了的

基本的蒙特卡洛方法（三）​

基本的蒙特卡洛方法（三）​

分层采样（Stratified sampling）

控制变量法

协稳的假设有多强？

修改App Engine上的cronjob

How to use VPN in Ubuntu 10.4

基本的蒙特卡洛技术（二）

Start to review a textbook with practice

全序集的本质

URL size limit of Google Chart API

不能在偏序集上排序

不等式移项

挑鸡蛋

GSL + MinGW + Eclipse

[Reading] Optimal Online Assignment with Forecasts

Space rail

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

Happy Chinese New Year!

The limitation of the sampling algorithm: a simple case

The impact of resampling to RSE in sampling algorithms

Share Two Google Calendar App

Ubuntu 9.10 remember last boot choice

Reading about Refactoring

2010 Resolution

Free SMS

时钟纸牌

Principle of Shrink Method in collaborative filtering

Implement RBM for netflix CF

RBM的实现

实验结果

The Stable Marriage Problem

Technical Details in RBM for CF

基本的蒙特卡洛方法（三）

基本的蒙特卡洛方法（三）