2011年1月11日星期二

马尔科夫链蒙特卡洛(二)

马尔科夫链蒙特卡洛(二)

来,做一道书后习题:用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年1月8日星期六

马尔科夫链蒙特卡洛(一)

马尔科夫链蒙特卡洛(一)

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

为目标平稳分布。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、不独立





2011年1月7日星期五

遍历trie

当时在网络所实习的时候,还没有学过数据结构,觉得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年1月2日星期日

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

中午接到电话,研究生同寝的哥们儿回北京了。于是一起去平安里吃烤鱼,一边吃鱼,一边研究一个随机数生成器的问题,一个学生问他的。
随机数生成器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

给我碗儿水喝

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