分类时一定只进行分类吗?

分类时一定只进行分类吗?

产生这个疑问主要两个来源于两个切身体会。

来源一

最近做一个项目过程中,遇到一个问题。由于实验室条件的缘故,对于某指标的检测只能在某个水平上检测出来,比如指标值为50以上才能有显示,于是实验人员会根据检测结果(一次检测会有多个位点的值测出来)来标注结果中是否存在缺失,他们的目的是想确定下在不同的水平下,缺失概率是多少?检测水平多少比较经济、合适?因为在50的时候,他们还需要去掉许多背景噪音,工作量大,想不做那么严格控制。对于这个问题,细抠起来想建好模其实还比较麻烦。不过我的疑问不在这里,而是在于由于不同实验室实验仪器和操作水平不同,阈值水平不一样,即使建模所得的结果,也是针对不同实验室得出来的结论。而造成这个结果不同的主要原因就是因为人为划分的阈值水平不同,导致类别结果不同。换句话说,我分类、打的标签不同,你建模的结果就会大不一样。同一批数据,我选择阈值不同,缺失的结果就不一样,所以我觉得挺难去确定一个适用于绝大多数实验室的阈值。
Continue reading

PKU暑期高维统计学习心得(II)

前言

距上一篇时间颇长,不过继续Jiashun老师的讲课心得。上一篇谈到稀疏、弱信号的一种处理框架——Higher Criticism,在分类、聚类等领域可以有比较好的应用。具体如何应用,此处不详谈,大家可以看看他第二节课的PPT以及该篇论文。在第二节课结束时,他提了一个结论:

Surprisingly, penalization methods (e.g., the L0-penalization method) are not optimal for rare/weak signals, even in very simple settings and even with the tuning parameters ideally set。

也就是说在稀疏、弱信号下,由L0衍生出来的方法并不是最优的,比较容易出问题。虽然我依稀记得某些论文模拟显示信噪比过低时候不少penalty方法结果并不太好,不过Jiashun老师的这个结论还是让我比较吃惊,毕竟被很正经的提出来了,而且他还有相对的解决方案!着实让我很感兴趣。
Continue reading

PKU暑期高维统计学习心得[I]

印象

为其两个周的北大关于高维统计的暑期课程即将告一段落,我回来奔跑了两周,身体略感疲惫,现在总算可以休息下,然后停下来消化下讲过的内容。

这次来讲课的老师学术能力都很强,都是四大paper等身的青年学者。老师们讲课的风格不一,最好玩的当属Tiefeng Jiang老师,他讲起课来就像说东北二人转,段子一个接一个,东北味的口音让我第一节课毫不犯困。而且深入浅出,随机矩阵这种比较数学的研究领域,也被他讲的比较好理解。不过后面由于有事情,以及之后的内容过于数学化,我就没有再跟下去了。Zhu Ji老师讲的很细致,不过内容偏简单了,听了两节课后我也没有跟下去。Cun-Hui Zhang老师做的很理论,深厚的数理分析功底,以及对高维问题理解的深刻让我感觉很敬畏,不敢靠近。不过对他后面做的scaled lasso和LPDE的结果很感兴趣,想用来做点检验的试验,不过邮件找老师要代码现在还没有回复,略感伤心,看来只能过几天自己写了。Yang Feng老师很年轻,在Fan老师那边做了很多非常好的工作,不过由于之前我看了不少Fan老师的东西,对他的讲的思路相对比较熟悉,也就没有太用心听而刷微博、做项目去了,真是一大罪过啊!

整个课程中对统计所持的观点和态度,我最欣赏的是Hui Zou和Jiashun Jin老师。
Continue reading

距离测度学习

上个学期快结尾时,我和小南决定找点东西做做。左看看右看看,觅了个距离测度学习来研究研究,这个在matlab中已经有了相应的包,只是功能并不是全,于是我们觉得把这个还算有点意义的东西写成R包,算作一个短期的学习。

不过人总是有惰性的,上个学期末我写了个距离测度学习的总结草稿,本打算寒假时候开始R包工作,但是谁知寒假毫无动力,堕落的只剩下吃喝睡觉。伴随着这个新学期的开始,毕业论文现在还没有眉目,距离测度学习这个东西小南却已经在github上开了个头,哪有半途而废之理。于是经过一天的梳理,我把这个总结的文档又重新完善了一下,虽然没有任何实例,但是对写R包还是很有帮助。现在赶紧挂出来,也算给自己一个督促吧。毕竟还有个毕业论文等着我呢。

距离测度学习文档总结

爱的声音

今天是母亲节。我是记不住这个节日的,可能本身这个节日在并不算做传统节日的缘故吧。长这么大,我一共陪母亲度过两个母亲节。一个是高三的时候,一个高四的时候。

那两年高考的最后几个月母亲过来陪读,很辛苦,很煎熬。但那是我第一次有了那么多时间与母亲交流,也是从那时候起,我才开始懂事,开始关心母亲,懂得一个小小的表达母亲都会很感动。每天的散步闲谈,让我深深体会到了父母对我的无条件爱,母亲说的有句话我一直刻在在心底“:只要你好,不管有没有出息,过的健康快乐,我和你爸就什么都好,都安心。”每次我一想起这句话,生活学习中的压力顿时烟消云散,母亲诚挚的话语,让我懂得了生活意义。虽然我无数次的质疑生活的意义,但是我从未放弃母亲告诉我的两样东西:爱与理想。

还有个场景我一回想就会热泪盈眶:高考前夕,烈日炎炎下,一群家长焦急等待自己孩子下课,手里都提着自己亲手做的可口饭菜,母亲因为身材娇小,而且怕热,每次都不得不站在人群后面,然后踮着脚,仰着头,远远的眺望,满脸的焦急与期盼。一想起这个场景,我心里都发酸,母亲为我倾注了那么多爱,而我却回报的如此至少,内心充满了愧疚与感动。

一晃,三年就过去了,如此之快。记得第一次我给母亲送了百合花,愿她健康快乐;第二次我送给她一只小乌龟,祝愿她长命百岁。我是一个不善于表达情感的人,满肚子的话,到了嘴边都会被憋回来。那两次送礼物的时候,我只是傻傻的笑,连句祝福不敢说,男孩子的倔强让我说不出那些肉麻的话。至今我也没对母亲说一句“妈,我爱你!”我知道,没有这些话母亲都很感动了,如果说了,怕我们两个都会哭出来,都是天蝎座的人,感情很丰富,很敏感,很脆弱,受不住这样的爱的表达。

三年的每个这个星期,虽然没法再送个我精心准备的礼物,但我都会或早或晚的打电话回去问候下母亲,今天也是一样。不过但是今天看到这个网页时候,我想到了一种表达爱的方式。我要一种比较特别方式送给母亲我一直不敢说的那几个字。不多言说,附图附代码。

1
2
3
4
5
6
7
8
9
10
11
12
library(seewave)
library(tuneR)
woaini = readMP3("我爱你.mp3")
woaini1=resamp(woaini, g = 2000)
muqinjie = readMP3("母亲节快乐.mp3")
muqinjie1 = resamp(muqinjie, g = 2000)
iloveyou = readMP3("i_love_you.mp3")
iloveyou1 = resamp(iloveyou, g = 2000)
op = par(mfrow = c(3,1),mar = c(4.5,4,2,2))
oscillo(woaini1, f = 2000, from = 3, to = 6,cexlab = 0.75,colwave = "red",title = "妈-我爱你")
oscillo(muqinjie1, f = 2000, from = 2, to = 5,cexlab = 0.75,colwave="gold",title = "妈-母亲节快乐")
oscillo(iloveyou1, f = 2000, from = 0, to = 4,cexlab = 0.75,colwave="lightblue",title = "i love you so much")

图解:

  • 最上是我用家乡话(信阳罗山话)说的:“妈,我爱你!”(罗山话中,“你”发音类似于“恩”,故振幅小)
  • 中间一幅依然用家乡话:“妈,母亲节快乐!”(录得不好,貌似有杂音)
  • 最下面的那副是英文:“I love you so much!”

其实从图的走势还是可以看出我羞涩的内心,爱的表达,对我来说,竟然真不是一件容易的事情,但是对许多男生来说,又何尝不是呢?男人有着内心的倔强,都懂得~

看过的“推荐系统”文献

酱油掉“树蛙”比赛后,深感遗憾。准备删掉这群“推荐系统”的文献前,还是做个记录,起码比赛前的那段时间看的那些推荐系统的文献,让我知道了推荐系统还是挺有趣,思想平易近人。更重要的是,让我等从未真实玩过”树蛙“的文盲,终于深刻体会到多挖挖数据比做一个模型重要得多。

  1. 我们最终敲定的模型,参考文献《Pairwise Preference Regression for Cold-start
    Recommendation》
    。只可惜当时用R写的代码慢的让人自杀,不得不放弃,于是用R整理好数据后,让同学代写C程序(大一时候学的C已经忘得差不多了,真是深感惭愧)。但是谁知道一直有个BUG没有解决,直接导致我们最后参赛结果:0提交!择人不善啊~~以后还是找个靠谱的熟悉矩阵的同学来写C代码吧。模型思想:根据推荐者与被推荐者的features与得分矩阵做拟合,得到推荐与被推荐者features之间权矩阵,然后预测得分。直接了当!模型评价:elegant!思想简约,有数学显式解,对推荐系统的冷启动问题效果应该不错。但是不知道对这次婚介的数据做出来效果是否好,真希望我们能够把那程序完善,去测试一次,了结心愿。
  2. 次选模型论文《Hydra: A Hybrid Recommender System》。也很不错,将需要的features信息都能够较好的连接起来,也是通过这个文章,知道了有个对等网络这玩意,也还不错。缺点就是在原基础上扩展的矩阵会很大,特别是对于这次婚介的推荐,SVD做起来消耗估计也会很大,不会太快,效果可能也比较一般。所以对于特征的选取还是得靠前期的数据了解,知道哪些人的特征在此数据中可以省下很多工作。跟这篇文章相关的文章很多,提高精度也罢,改进模型也罢,反正看着挺好,就是不知道用起来如何...sigh~~
  3. 《A Singular Value Decomposition Approach Recommendation Systems》这篇文章还是讲SVD的,但是感觉看了那么多论文,感觉这篇文章前面的Introduction介绍推荐系统比较系统全面,如果对推荐系统问题有了较好的定位,我觉得不管做哪些问题,思路都会清晰很多,把握住了问题的关键点,这样结果才会有好的突破。当然这篇文章也详细论述了SVD用在推荐系统的实际意义。
  4. 最后再提一个我和小南都认为很不错的论文《Learning User Preferences in Online Dating》,思想很简单,高中生也能不过很能够解决问题。当然这篇文章有个致命的缺陷:数据获取,要能够获得双方发MSG时候的信息。这个实际中估计是是不可能的,所以这篇文章只能提供一个比较好的思路,不能解决实际问题,仅供把玩~~
  5. 至于协同过滤的东西,看了挺多,觉得这个东西重要的不是数学,而是他的思想,比如从USER变换到ITEM就可以获得一个很大的突破,真的是让人感觉很畅快。

比赛早已落下帷幕,小南说他已经把PAPER全扔进了垃圾桶,我也想删掉这群文献了。虽然结果让人很XX,但是没有关系,起码都获得了很多,这点只好深埋心里,笑而不语也。看排行榜龙争虎斗也是快事一件!

纠结的问卷设计

转眼就是三个月过去了。

记得还是年前那会,我和小南在“八哥”小店吃饭,两人商量着美赛的事情,我一时兴起,想了个无聊的点子:互相帮对方设计下未来1年半的路途,立下字据,1年后看完成的怎么样。于是操起桌上的菜单,刷刷龙飞凤舞起来。paper,MCM,信安比赛,paper,旅游,TOFFEL,GRE,paper,毕业旅行...一路写的眼冒金光,仿佛世界已被我们主宰,然后顺风顺水的过完了整个大学...那时,两个年轻人真是幻想加激情,笑声都轰天震地。

今晚,我翻开钱包,那张纸虽然破烂,但是还是被我好好留存着,看着上面满满的计划,不禁内牛满面。想起小南那篇<successsful participant>的文章。狗屁般得论文与绚烂的排版形成的鲜明的对比,让我满心的愧疚与伤感。再回头看看数据挖掘比赛,更是让人无语凝噎。只好作罢,留作下次来谈谈对推荐系统的心得了。今晚,先聊一个自作孽的抽样调查。

其实这个本是开学前就该完成的课程设计,但是当初雄心勃勃,想写出一篇“唱作”俱佳的好文章来,于是一腔热血的跑向中南的湘雅医院,找了个心理学研究生准备做个“人格与人际交往”方面的研究。后面也定好了位,研究“大学生人格因素与人人网人际交往以及现实人际交往的关系”。本以为可以借着现在的社交网络研究热来给心理学也点点火,本以为心理学测量这东西水不深,本以为一切都是顺水推舟的事情,但是“一入侯门深似海,从此萧郎是路人”。从此我对心理测量敬而远之,还是看看实用心理学罢了。

当然心理测量并不是不好,而是我通过这次做问卷设计才知,原来一份好的心理学量表需要耗费如此多的功力。按那个湘雅研究生的说法,湘雅可以算是国内心理测量方面的发源地,但是做心理测量用的量表基本还是从国外翻译过来,最多做一点点修改。一份好的心理学量表是要经过无数次的检验才能达到很好的信度和效度。这次我们调查用了两个问卷,一个问卷是“卡塔尔16pf人格量表”,一共187道题,经过几十年的修正后已经非常成熟,很权威。而另一份量表的我们自己做,正是这个量表的设计,让我知道我的抽样学的很差,再次深刻的理解到抽样重要的是各种方法的思想和应用范围,而不是那些重复相似的数学公式,更重要的是,抽样只是一个辅助工具,就像推荐系统一样,算法终究不是大头,重要的知识。对于心理学量表设计,如果没有一定经验和对心理学一定的理解,是没法做出一份像样的问卷的,就如我们这次做这个问卷一样,才返工了两遍,就把被试者做烦了,还是认识的好朋友呀!如果是一般的受试者,估计不给钱不给个礼物,真是难以修改重测,而且修改也是困难重重,像我们这种心理学门外汉,凭借着一点文献和心理学书籍,哪有可能几遍修改就能做好。烦躁激动纠结集于一身,于是乎问卷设计以及论文...后续略去几千字吧...

其实通过这次做问卷设计,我认识到问卷搜集到的数据其实可以是很可信的。只是许多没有相关专业背景和统计学背景的人把问卷胡乱一通的用,导致人们一般认为得来的问卷结果是受试者敷衍的回答罢了。对于心理学量表来说,这些问卷的设计很多很精妙,除非如果你自认为你能考虑这些问题的设计会比心理学家更加全面,但是一般没有受过专业心理学训练的人是不可能做到的。一份成熟的量表,我们考虑的问题,都已经被心理学家考虑到了。当然问卷是一个双方互动的过程,很多测试问题我们感觉很浅显,都知道他问的方向,但是此时我们需要做的就是诚实做答,这样一份好的问卷才能收到一个好的效果,不然问卷就没有任何意义,一个人也永远别想从专业心理学角度了解自我,发现自我。

现在的社会人们更加机灵了,不愿意这么明显的在问卷中暴露自己,怕被骗或者被利用,这种保护心理是可以理解的,但是殊不知每天我们那些无意识的网络行为其实也被人利用了呢?无意识的暴露和有意识的暴露都不重要,因为都可破解,重要的是我们愿不愿意配合各种调查的进行。

问卷设计和问卷调查是一门艺术,或者说很多文科性质的东西也是艺术,因为他们研究的对象很多都是人或者人的行为,难以捉摸这才是魅力所在。不是所有的东西都可以被定量化,如果都程式化研究,也就没啥意思了,人不是电脑,不是数学方程式呀!

不过话说到最后,我还是喜欢定量化研究,因为我喜欢数学,简洁,具有形式美。

误差论->回归(1)

###竟然发现这篇辛辛苦苦写的文章不见了,1年前我的理解的还是挺正确的 by 2011.9###

最近一直不能很好的理解回归中最小二乘的意义,首先就被一个傻傻的问题困惑:假设试验中没有系统性的偏差或者人为的失误,而误差是不可控的,为何要对误差求其平方和的最小值?这个最小到底意义何在?

虽然书上回归分析所得的结论都是很自然,但由于这个方法的前提(历史)我不十分清楚,内心甚是不安,总有种被人灌了酒,模模糊糊接受了结论之感。于是立马查阅相关史料和名家所写的各种回归分析的书籍,梳理下对回归分析的认识。

一. 误差论的过往
其实可以这么说,回归分析是线性模型中误差分析的推论而已,只是之前没有一个完整的数理统计框架,没有一个对回归分析意义的清晰解释罢了。因此,要想很清楚的理解回归,必须要从误差论开始,从最小二乘法的历史来源开始。

最小二乘法最先由数学家勒让德在测量子午线长度时候,为了解决方程个数n大于未知数k的矛盾方程组发明的方法(这里还有一段与高斯的发明最小二乘优先权的故事,勒让德是1805年发现,理论上还是比高斯的1809年早)。这个式子是

 x_0+x_1\theta_1+\dots+x_k\theta_k=0


现在我们按代数的观点看的话,如果要解出各个\theta_i的精确值,很可能是不存在的,历史上像欧拉、拉普拉斯这样的大学者也尝试过解决这个问题,可能是他们更习惯于求解提法严谨的数学问题,因此对于这种实用性的数据处理问题,一直没有什么建树。直到勒让德发明最小二乘法,这个问题才算得到了不错的解决。

勒让德的思想是考虑误差在整体上的平衡,即不使误差过分集中在几个方程内,采取

 \sum_{i=1}^n=(x_0i+x_1i\theta1+\dots+x_ki\theta_i)^2=min


的原则去求解\theta_i
其实勒让德思想看起来是挺自然的,但是我却仍然困惑,最小二乘得到的结果顶多是个最好的估计罢了,而误差是不可控的,你怎么能使其平方后相加能够变的小呢?所以此时的最小二乘法只是个算法,并没有体现什么更深刻的思想。因此我们还是得沿着历史的轨迹,去寻求更完美的解决方案——误差与正态分布。

误差论的研究最初起源于天文学,而天文学之于数理统计发展,丹麦统计史学家哈尔德在其著作指出:“天文学自古代至18世界是应用数学中最发达的领域,观测和数学天文学,给出了建模的及数据拟合的最初例子。在这个议一下,天文学家是最初一代的数理统计学家……天文学的问题逐渐引导到算术平均,一集参数模型中的种种估计方法,以最小二乘法为顶峰。“

这里既肯定了最小二乘法的显赫地位,同时也指出这种地位的确立,极大程度上取决于误差理论的建立,如果没有这样一个理论,最小二乘法只是一个算法,没有了与统计分析相联系的纽带,于是就会依然回到之前那个问题——对不可控的误差求最小平方和的意义何在?

误差理论研究的基本问题是指:随机测量误差服从怎样的概率规律,既有怎样的概率分布?

其实我们根据现代的概率论观点,由拉普拉斯中心极限定理,很容易知道误差的分布是符合正态分布的。当然现在看来是十分简单的,但谁都知道,那个时候还没有严谨的数理体系,早期天文学家根本不知道中心极限定律是啥,现在看起来简单的问题,从前可能就是一个难题,可见科学的发展是多么的不易。因此了解误差理论的发展,对回归的认识非常有益。

随机误差的开拓者——伽利略从第谷的观测数据,提出观测误差的3个条件:
1.所有的观测值都可以有误差,其来源可归因与观测者、仪器工具及观测条件。
2.观测误差对称的分布在0 的两侧(排除系统误差)
3.小误差出现的比大误差更频繁。

综合这几条,伽利略所设想的误差分布,就是一个关于0的对称分布,概率密度f(x)随|x|增大而递减。这个原则性提法,成为日后学者们研究这个问题的出发点。当然由于当时概率论发展水平所限,伽利略们没有取得什么进展。真正的误差理论工作开始于Simpson,而正是他提出用平均值来估计真实值。

他的尝试是在误差满足某种分布下,计算平均误差的分布,从而证明在某种概率的意义上,平均误差小于个别误差。虽然他的工作并未触及建立一般的误差概率理论的问题,但是有一点很重要,就是撇开位置的真值不论而把注意力放在误差上。这样就把问题给简化了。

随后拉普拉斯也进行了研究,与Simpson不同,不是先假定一种误差分布后设法证明平均值的优良性,而是直接涉及误差论的基本问题,即误差分布是怎样的,有了分布后,再如何根据未知量\theta的多次测量结果估计\theta。当然这是一条正确的思路。可惜的是,拉普拉斯没有想到用他的中心极限定律,而当时也没有估计参数的矩法和极大似然发,于是他用了贝叶斯式的推理方式去求解,但是由于他的误差概率密度分布(现在我们常见的拉普拉斯分布)不对,因此后面的计算非常复杂,后面他虽然仍在这个问题上耕作了几年,但是仍没有得到理想的结果。

这时终于要轮到高斯登场,高斯正态误差分布也即将横空出世。此处简单介绍下高斯求解的过程。
首先他用到了最大似然估计法(当然费歇尔之前是没有这个叫法的,可以说高斯首创了这个方法,当然不够一般而已)。设真值为\theta,n个独立测量值X_1,\dots,X_n,高斯取概率

 L(\theta)=L(\theta;X_1,\dots,X_n)=f(X_1-\theta)\dots f(X_n-\theta)

然后取使是概率最大的\hat{\theta}作为\theta的估计,同时先承认算术平均值\hat{X}为应取的估计,然后高斯证明,这样的取法只有误差的概率分布为正态分布N(0,\sigma^2)时成立。具体证明略去,有兴趣的可以自己根据伽利略所提的误差概率性质,再结合一些函数知识便可以证得。然后再根据概率最大原则便可得出最小二乘估计式,至此,高斯便赋予了最小二乘了更深刻的意义,把它与数理统计结合了起来。

当然那个时代高斯正态分布还不能体现其深远的影响,直到20世纪正态小样本理论充分发展起来以后才充分显现。但此刻我反观高斯正态分布,再看看现代的数理统计结构,由高斯导出的理论简直充满了我的眼球,很多地方都可以有他的身影,那回归分析简直就是高斯重生了!所以我不得不感叹,高斯真真真是个天才!

但是崇拜之余,这并不表明了高斯创建的误差理论已经十全十美,细细品味,高斯的论证有循环论证的味道:由于算术平均的优良,推出误差必须服从正态分布;反过来又由正态分布推出算术平均及最小二乘估计的优良性。而算术平均到底说来,总是没有一个自行成立的理由,只是从古到今,人们除了平均下,没有太好的选择,而以它作为理论中一个预设出发点,终有其不足之处。直到后来拉普拉斯终于开窍,用了自己的中心极限定律轻松证明了各种误差之和组成的随机变量近似正态分布后,误差论终于完美建立,随后高斯又证明了一系列最小二乘估计的优良性质,其中最著名的就是高斯-马尔科夫定理,而这个定理的高明之处在于,说明了没有啥估计还能比用最小二乘的得到的参数估计方差还小了。换句话说,我就是稳定可靠,没有比我更好的了。

至此,我们可以看到误差理论的发展已经能够达到了很深的水平,如果高斯他们知道高尔顿的工作,那么真的不再是高尔顿发现回归,而是高斯了。当然话也不能这么说,回归需要对事物间的关系有更深层次的认识,也许从数学理论上来看,回归分析不过是误差理论的简单推论罢了,但是从思想上,它是无法以此类推。科学史的发展就是如此,理论数学家已经把数学弄得出神入化,数学的东西一般都不缺,缺的是一个创新运用的思想。

总结:误差理论的建立,使最小二乘法不再仅仅是一个解决矛盾线性方程组的算法,而且让他有了更深刻的统计意义。这也是为什么我不太接受勒让德的最小二乘,而接受高斯的最小二乘的原因,最大似然更能触及概率的本质。而高斯建立的正态分布的误差理论的全部影响需要直到正态小样本理论充分发展起来后充分显现出来,同时他也让后面的回归有了更深刻的背景。按我的理解,回归分析其实就是各种影响因素叠加后,在宏观方面上对几个因素关系的分析,说到底,他只是求条件期望后平均值意义上的误差分析。这也是为什么一旦这些随机误差的假设条件不满足时,回归方法需要改进。因为高斯的最小二乘法是建立在随机误差正态分布上的,如果不满足随机,不满足正态性,不满足独立性,最小二乘法很容易失效,回归效果会很差了。因此我们会用岭回归、LASSO、M估计等方法来解决,而至于回归的深刻含义,留作下一节继续。

P.S. 今天老师讲解计量时突然讲不清楚 残差和为何要用\mathcal{X}^2检验,我感到很是悲哀,因为这正是误差正态假设的前提下很自然的结果。

数据和理性

奥运420亿美金、60国庆开销…不明,反正上了百 亿、世博4000亿RMB、09年全国新开工高速公路计划总投资约 7000 亿元,沪杭磁悬浮上1000亿只为节约十分钟…然后我党给 水深火热中的西南五省下拨抗旱经费—…1.55亿元。

这是我从一个同学的博客上抠下来。其实我也没啥意思,只是最近不管是寝室内还是寝室外,总听到不少人在骂dang不人道,不太管西南老百姓的死活。而一同学天天粘着猫扑上,更是大骂社会不公。

当然我没有那么愤青,听着这些话,总觉得网络语言中有太多的戾气和不理性。人们对数据的比较始终充满了主观性,特别是数据的比较,缺乏统计人的理性。其实我也不是为dang做什么辩驳,只是觉得上面这种数据的比较除了激起不明事理的人的愤怒感,没有太大的作用(当然会促使dang拨更多的钱给西南,但是难免其中又有贪污了)。国家的财政预支都是上一年就做好的了,如果咱国把应急预险的预算做到和世博一个价位的,不知道有多少人会跳起来大叫“疯了!疯了!”国家有国家的打算,个人有个人的看法,如果我们总是从自身出发看问题,缺乏大局的眼光,那么我们总是会被别人的思想牵着走。

有些人觉得国家搞个奥运会,弄个世博会,再折腾个磁悬浮,花了那么多钱,仿佛是关心形象工程,而对于真正的民生,比如西南大旱却吝啬如高老头。对这我不想发表太多的言论,一个中心两个基本点,主要矛盾与次要矛盾。国家要发展,强健民心的工程必然是不可少的,而对于自然灾害,国家肯定有自己的打算,目前据说国家已经拨了20几个亿给西南了,这对于缓解居民的基本生活困难帮助应该不小了,但是这肯定无法从根本上解决问题,钱不能让老天马上下大雨,根本问题是如何加大水利建设和环境保护工程,增强抗旱的能力。所以20多个亿了,数目还是不小的,毕竟对于一次自然灾害来说(当然无法与汶川大地震相比较的)。

不要看到和世博那比起来这些钱是那么的寒碜,然后就怒气冲冲的找国家说理。一同学说的好,“西南大旱又不是建磁悬浮引起的。”一群人抓着个数据就像抓着个宝样,这不是个好兆头,特别是一群不懂数据、不明事理的人(当然我也不明事理),这样只会让自己愤青升级,一天到晚张着个嘴骂这骂那,而什么实事都不做。问题是要来解决的,关键看你用什么手段解决。咱要理性点,是好是坏,咱也不是瞎子,咱也是明了的。
所以这样看来,在众人中普及点统计知识是迫在眉睫呀,特别是一些愤青们和一些“领导上”们。所谓心中有牛屎,看什么都是牛屎。数据也可以是理性的,关键是我们有没有理性了。唔,路漫漫其修远兮,吾将上下而求索之...

提高Monte Carlo模拟定积分计算的精度

前几天在COS主站上看到一篇谈Monte Carlo求定积分的方法的文章。记得去年我学matlab时,求定积分便用到了这两种方法,而在上概率论课时,老师也略微涉及到相关的内容。当时的感觉时,定积分这种一般用确定型算法计算的东西,竟然还可以用随机的方法来实现,确是一种非常巧妙的思想,(Monte, Carlo是摩纳哥的一个赌城,二战时Von, Neumanm一伙人研究原子弹时弄出来的方法。)于是被震撼,也没有再深究其中更详细的内容。

但是昨天一看到这篇文章,便想起寒假时候看的那本《计算方法》中关于Monte, Carlo平均值方法的论述,而最近在上计算方法的课,正纠结与一堆无聊的误差计算中,于是便再次翻开此书看看,竟发现了一些一直没太注意的很重要的章节。

具体如下:

  • 平均值方法是根据弱大数定律证明得来的。构造如此

    I(f)\approx 1/N\sum_{i=1}^{N}f(X_i)=I_N(f)


    均方误差

    E|e_N|^2=E(I_N(f)-I(f))^2=\frac{1}{N}Var(f(X))


    而由Schwartz不等式知

    E|e_N|\leq \sqrt{E|e_N|^2}=\sqrt{Var(f)/N}

    ,从中便可知 MC 法具有半阶收敛性,且收敛阶不依赖问题的维数,这样在计算高维积分的时候,比较优越。但是收敛速度还是很慢的,是一个比较坏的收敛速度。因此在有其他高效率、高精度的确定型算法时,就不要采用MC法,那是在没有更好的算法时更好的抉择,比如超高维积分。
  • 由均方误差可知,要提高 MC 法的精确度,那么久必须得增大N,减小VAR(f),也就是说,要增加样本量,同时改善抽样方法的选择,这时便有很多方法选择,书中提到了四种方法,我就谈谈我的理解和试验。
    1. 分层抽样。核心就是将原来的 [0,1] 区间 M 等分后在各个小区间上分别再以均匀分布抽样 N/M 个点,然后以

      I_N=1/N \sum_{k=1}^{M}\sum_{i=1}^{N_k}f(X_i^{(k)})


      理论证明知这样方差减少,具体证明式子较多就不再写,否则就是一本教材的翻版了。
      模拟\int_{0}^{1}exp(-x^2/2)/\sqrt{2pi}dx,具体实验结果如下 

      1
      2
      3
      4
      5
      6
      7
      
      M&lt;-100
      N&lt;-1000
      A=matrix(as.numeric(gl(M,N/M))-1,ncol=N/M,byrow=T)
      B=replicate(M,runif(N/M))
      x=1/M*t(B)+A/M
      y=exp(-x^2/2)/sqrt(2*pi)
      I=mean(y)

      结果为0.3413287,可知只取了1000个样本点的情况下,精确度已经十分好了,而一般的均值模拟
      达到此精确度则需要10^7个点。

    2. 对偶变量法。
      这是一种特殊的针对于定义域具有一定对称性的函数性质的技巧。基于这样的命题:如果f(x)是单调的,则Cov(f(X),f(1-X))\leq0,这里X\sim U[0,1]。则有

      I_N=1/(2N)\sum_{i=1}^{N}(f(X_i)+f(1-X_i))

      ,
      EI_N=I(f),
      经过证明知 Var(I_N)\leq \frac{1}{2N}Var(f),方差显然减少。 

      但是这种方法只是对于定义域有一定对称性的函数,如果对一种函数先验了解不多,而且对称性不知,则运用会出现失误,所以具有一定局限性。
      模拟\int_{0}^{1}exp(-x^2/2)/\sqrt{2pi}dx的代码如下:

      1
      2
      3
      
      a=runif(10^6)
      y=exp(-a^2/2)/sqrt(2*pi)+exp(-(1-a)^2/2)/sqrt(2*pi)
      b=sum(y)/(2*10^6)

      结果为 0.3413487,精确度也得到了提高。

    3. 重要性抽样:核心就是构造一个新的积分

      I(f)=\int_{0}^{1}f(x)dx=\int_{0}^{1}\frac{f(x)}{p(x)}p(x)dx

      其中p(x)[0,1]的概率密度,\int_{0}^{1}p(x)dx=1,p(x)\geq 0。这样得到类似于最初最简单的那种均值法的新的形式积分近似 I(f)\approx 1/N\sum_{i=1}^{1}f(Y_i)/p(Y_i)Y_i为遵循p(y)的独立同分布随机变量(i.i.d),而之所以这样构造,可以根据这幅图来理解,a importance sampling chart\int_{0}^{a}+\int_{b}^{1}f(x)dxI(f) 很少比例(即图像的边缘部分),尽可能的不让随机点去计算这些值,这样的“拟合”积分较好(面积),否则很多的随机点在 [a,b] 外,为计算小比例的东西做贡献,这样必定精度较差,所以我们要期望尽可能的值出现在峰值处。我们称其为“重要性”抽样。理论分析,

      Var_X(f)=\int_{0}^{1}(f-i(f))^{2}dx=\int_{0}^{1}f^{2}dx-I^{2}(f)

      ,

      Var_Y(f/p)=\int_{0}^{1}(f/p)^{2}pdy-I^{2}(f)=\int_{0}^{1}f^2/pdy-I^{2}(f)

      .
      只要选取适当的 p(y) 就可以使方差减少。
      当然不可能是Var(f/p)=0,这在图的模拟便可知,p(y)不可能与f(x)相合的同时,还有\int_{0}^{1}p(x)dx=1。但是p难以选择,需要对问题先验有了一定的了解程度才能适当选取。比如对于正态分布求\int_{0}^{1}exp(-x^2/2)/\sqrt(2pi)dx,我不太好选p,就没有做相关的模拟了。
    4. 控制变量法。
      核心也是构造,数学式子I_N(f)=i/N\sum_{i=1}^{N}(f-g)(X_i)+I(g).
      需要I(g)已知,但是这是同重要性抽样一样,都不知道,也需要对问题先验有了一定的了解程度,才能好的选择,使I(g)已知,并且有Var(f-g)\leq Var(f)。相关模拟我也没有再做

总之,正如书中所言,减小误差都只是具有指导思想的意义,具体问题,往往需要我们对被积函数有先验的了解,才能够减小误差,提高精度。确实如此,我的例子都是比较特殊或者简单的,如果对于复杂的问题,而且对该函数比较怪异少见,那么在以上的方法中寻找什么p(x),都是一项非常艰巨或者不可能完成的任务。因此,我下一结论,奇技淫巧之所以巧无外乎我们有了个千里眼和顺风耳,瞎碰瞎撞弄出一个巧妙的解法那是很难的,大部分都是前人解法突然激发灵感得来(此时此刻竟想起了贝叶斯来...)。所以对于专门解决高维积分强有力算法——Metropolis算法我是敬佩有加,这一算法2000年惊被评为20世纪十大算法之一。后面产生的什么模拟退火,拟Monte Carlo都与他有莫大的渊源。而伟大的东西终究充满了高深而深刻的东西,我目前楞死没有看太明白,已经被符号给堆死了。但是我想过段时间在查阅相关资料应该可以弄明白吧。

最后再扯句从中学到了的随机向量的构造——正态分布随机变量——Box-Muller法
以前学的是用n个均匀分布随机数由中心极限定律构造,但是那个要得到较好的精度需要N很大,故需要精确的还是用上面的那个方法更比较有保证,咱看公式也心安嘛。