趣物

以前闲着无聊的时候曾经做过一个《资治通鉴》的字频统计，单以频率计，中国历史不过是“王”与“人”，“义”与“忠”，“将军”与“刺史”，“长安”与“洛阳”。

既然有了频率，自然也就有了概率和条件概率。根据条件概率，当给出一个序列的字词后，预测下一个字词是什么，就变成了一个简单的最大似然估计问题。如果觉得这个序列太长，计算起来太麻烦，可以假设简化的马尔科夫结构，譬如假设下一个词的概率取决与之前的n个词而不是整个序列，这基本上就是计算语言学里的 n-gram 算法了。

所以我们可以用《资治通鉴》作为语料得出经验条件概率，然后来随机模拟出历史文本，产生原汁原味（至少是统计意义上的）史书 （技术细节见附录）。 虽然这只是文字游戏，但是仍然能从概率上看出《资治通鉴》记述的历史中，最容易重现怎样的事件。

譬如下面这则 （random seed = 2000）：

撰 刘 崇 俊 以 惟 岳 又 从 入 关 ， 宣 等 从 太 子 也 ， 惧 履 危 亡 之 事 ， 发 步 骑 二 十 骑 自 北 至 北 寺 狱 ， 竟 不 使 宗 庙 社 稷 。 宗 元 为 柳 州 司 马 。 坚 大 怒 ， 士 气 彫 沮 ， 无 事 ， 更 为 后 拒 ， 倍 急 于 亡 命 聚 众 二 万 会 麻 秋 、 姚 、 宋 赤 眉 将 逢 安 为 新 都 ， 剽 掠 。

我们可以这样翻译：

刘崇俊因为惟岳（人名，可能是李惟岳。刘崇俊是五代人，李惟岳是晚唐人，相差不算太远）一起入关（姑且认为是潼关，但是李惟岳在和河北，真实历史是不会入关的），而宣（人名）等人却追随太子而去，害怕这是危亡社稷的事情，于是发步骑二十骑（区区二十骑！估计是武林高手），从北面到北寺狱（这是东汉时候黄门署属下的监狱，鞫禁将相大臣的，好吧，晚唐也有宦官之祸，这里东汉的宦官乱入了，不过二十骑到北寺狱，难道是要劫狱？），最终也没有拜谒宗庙（不把皇帝放在眼中啊）。

柳宗元被任命为柳州司马（二王八司马嘛，承接上文，还是和阉祸有关）。

（某）坚大怒（不知道是孙坚？苻坚？杨坚？），士气不振，好在也没有什么大事，继续抗拒（王师？）。 因急于亡命，聚众两万人与麻秋汇合（麻秋登场，那么坚应该是苻坚了，但是麻秋很早就被苻坚的伯伯苻健杀死了）。

（与此同时）姚、宋（姚崇，宋璟？）率领赤眉军把逢安（这个地名是自动产生的，历史上似乎没有，权且当作是四川 蓬安 吧）作为新的都城，四处剽掠。

我们来梳理一下这段模拟历史的脉络：

这大概是一个王朝末年的乱象。 地方农民起义（赤眉），建立政权（所以有新都），负责讨伐的将领反而形成军阀割据，一些军阀随权臣（刘崇俊）入京，干预朝政（惟岳），一些军阀在地方反叛（坚），勾结外敌（麻秋是羯族人）。这一切的原因可能是因为朝廷宦官弄权已久（北寺狱），忠良被贬（柳宗元）。 军阀入京大概是打着清除宦官的名义（所以要发兵北寺狱），但是同时他们也不把社稷放在眼里。 京城在军阀到来前似乎已经被反贼攻克，所以皇帝和太子分道逃亡。如今皇帝似乎已经回到京城，而太子却还在外面招兵买马（宣等人追随），似乎有不臣之心。

简而言之： 中央朝政腐败，宦官专政，两宫不和。地方盗贼风起，军阀割据，外患不断。

难道这就是随机生成的中国历史最典型的场景？ :-)

附录

文中使用的通鉴文本是从维普网上下载的，做了一些简单的清理，上传到了百度云（链接)。 为了避免古文分词的麻烦，在作 n-gram 的时候以字为单位，所以用 gsub 在每个字的后面加入空格。 n-gram 选择 n=2.

 library(ngram)        
 file<-"C:/Users/Zeng/Downloads/zztj.txt"     
 str=scan(file,"character",encoding="GB2312")      
 s = concat(str)      
 s = gsub("(.)", "\\1 ", s)         
 ng = ngram(s, 2)  
 o = babble(ng, 100,seed=2000)  
 Encoding(o)<-"UTF-8"  
 o  

π 节（3月14日）已经过去了，所以现在写 π 就当是向魔幻圆周率致敬吧。其实数学的发展最需要牛人，就像祖冲之，冷不丁的从历史里给你扔出一个 355/113 (=3.1415929204) 来，等你回过神来，这个 355/113 是如何得到的已经湮没在历史的尘埃中了，所以牛人之后还需要不太牛的传承者。也是因为中国的数学研究一直缺乏传承（后来很大一部分因为科举），在元朝时数学大跃进之后，到了明朝，就大踏步的倒退到了简单算术（虽然珠算得到了很大的发展）。

相比较而言，欧洲的数学一直到十六世纪初都不过尔尔（古希腊的数学传统到那时已经断掉），但是在突然冒出了笛卡尔，莱布尼茨，牛顿，伯努利等人以后，数学的发展风驰电掣，而到了欧拉，简直是神一样的人物了。

莱布尼茨给出了一个用无穷级数求π的方法:

pi= 4*sum (-1)^n/(2n+1)

在 Q 里快速的实现一下：

 {sum raze 4%(1 -1)*/:(0N,2)#1+2*til x} 1000000
 3.1415916536

上例一共用了一百万项，精度略低于祖冲之的密率。

然后是欧拉用 p-级数的方法：

pi = sqrt( 6 * sum (1/n^2))

在Q里快速实现一下：

` {sqrt 6 * sum 1 % {x*x} 1_til x} 1000000 3.1415916987

同样一百万项，收敛速度和莱布尼茨的无穷级数差不多，但是作为大神，欧拉进一步证明了上述p-级数与素数乘积的关系:

sum 1/n^s = prd 1/(1- p^-s)   p 是所有素数

当 s = 2 时， 我们就可以得到

prd 1/(1-p^-2) = sum 1/n^2 = pi^2 /6 

这等同于

如果我们从 1 到 N 的N个自然数中随意抽取两个数，那么这两个数互质（最大公约数为一）的概率在当 N 趋近于正无穷的时候趋近于 6 / pi ^ 2。

以上推论的证明很简单，如果两个数互质，那么它们没有共同的素数因子。任意一个自然数被一个素数 p整除的概率是 1/p，譬如，能被5 整除的数字每间隔5个出现一次。任意两个自然数能被 p 整除的概率是 1/p^2. 至少一个不能被整除的概率是 1-1/p^2。把这个结论扩展到任意的 p，那么这个概率就是 prd (1- 1/p^2) for all p。 简单的变形，我们知道这个概率是 6/pi^2。

简单的蒙特卡洛一下：

{sqrt 6% avg 1={$[y=0;:x;:.z.s[y;x mod y]];}./:(0N;2)#(1+ x?x)} 1000000
3.1431107401

从一百万个自然数中随机选取50万对，然后算互质的概率。上述模拟中求最大公约数时用的是欧几里得算法:{$[y=0;:x;:.z.s[y;x mod y]];}。

转了一圈，终于回到了古希腊：）

开始继续写blog，用简书作为国内的镜像，因为我的blog还在被防火墙中。现在想起来，blog被封也该有近十年了，起因是清朝，后来清宫戏大热，也挺戏剧的：）

好吧，为了不让这个镜像再惹麻烦，决定不再评论嘉庆皇帝因为恪守着“明亡于万历”的祖训，恪守着“滋生人丁永不加赋”的旧制而无法应对由于大量白银涌入引起的通货膨胀而导致的财政枯竭，这样太让容易被误认为是在“打比方”（尽管现在最大的问题大概是通货紧缩）。

查询了一下blog的MySQL数据库，得到了上面的统计数据图。2006年写了大于300篇的blog，可见那时有多闲，那时候应该已经过了资格考试，正在准备写论文，然后07年开始准备写论文，08年开始准备写论文，09年开始准备写论文，10年论文居然神奇的写完了。写完了论文毕了业，就没有时间写 blog了，12年只写了9篇，13年1篇，14年0篇。于是上图呈现出对数正态分布：）

作为正在文盲化的资深失语症患者，觉得还是需要经常写写字，有助于防止成为痴呆的文盲：）

（这是一篇关于很枯燥的技术，很枯燥的历史文本，和不太枯燥的统计的 blog）

看过一篇关于《全宋词》词频统计文章，挺有趣的，想用类似的方法处理一下《资治通鉴》，所以就趁周末花了几个小时作了一下。

词是长短句，统计两个字组成的词频比较合适，《通鉴》是古文，文字结构不同，所以我统计了单字频，两字词词频，三字词词频，四字词词频，和五字词词频。同时也记录各个统计单位（字或词）出现的卷数。《通鉴》294卷，从三家分晋到五代结束共共1362年，所以卷数可以作为时间的度量。

《全宋词》的词频是用 R 作的。R 虽然是不错的统计软件，也是我的最爱之一，但是 R 并不适合作文本分析，更不适合来作数据库操作。所以就用了 C# 和 Kdb +3.0。 C# 用来分析文本，.Net 是懒人的福音，并且多线程运算非常简单，能够大大提升文本处理速度，Kdb+用来储存数据，它差不多是性能最好的 in-memory 数据库了，从它的网站上能下载到免费版本。这个分析里数据库是重头戏，因为需要查询数百万行的数据 row，如果用 MySQL，估计会龟速到死。另外 Kdb + 本身只有 300多K，不用安装，很方便。还有就是 Kdb+ 的 Q 语言也能满足编程需要。

Kdb+ 的网站提供了各种语言 API 的源码，C# 的 API 不支持多线程，所以需要在适当的地方加锁。Kdb 唯一的问题是不支持 UTF-8。它用的是 UTF-7，所以在注入中文文字数据的时候可能会出现乱码，为了省事，从 C# 里 publish 数据的时候，直接 publish 为三字节的 int[] 了。query kdb 时用了一个免费的 GUI QPad。QPad 似乎是用 Java 写的，它的编码默认是 UTF-8，所以在 query  Kdb 的时候直接把三字节的 int vector  cast 成 char，在 QPad 里显示的就是中文了，所以也很方便。

产生数据的 C# 代码非常简单，发布数据的时候自动生成 Kdb 的 schema。使用的《资治通鉴》的文本是网上广为流传的国学网简体版，在生成数据前，先用 C# 作了预处理，主要是用正则表达式替换掉了现代语言的“污染”（譬如：“后一页”，公元xxx年 等）

下表是各个字、词频的数据量：

从上表看，《资治通鉴》应该有近两百六十万字。

单字的字频统计如下：

“之”字当之无愧的排在了第一位。第一个非虚词是“王”，它包含了姓和爵位，第一个动词是“曰”。“人”的频率也很高，“将”，“军” 在双字词频中也会遇到。“帝”字排名32，“后”字排名33。但是因为是简体字，“后”并不专指皇\王后。

下面是价值观念的排名：

义、忠、孝排名在前，智排名最后，倒正印证了司马温公那句话：“凡取人之术，苟不得圣人、君子而与之，与其得小人，不若得愚人。”当然这个数据里噪音很多。

另外还有很多有趣的东西，就不一一叙述了。下面看一下两个字的词频：

“将军”出现的频率最高，在第一卷里就出现了，“度史”显然是“节度使”里出现的，虽然在203卷才出现，但是它居然出现了3202次，唉，唐朝啊！“节度”一次出现的要比“节度使”早。“皇帝”一次最早在第6卷出现，其实那时还是昭襄王元年，但是因为文本中出现了“秦始皇帝上”。

两字地名出现的最多的是“长安”，排名43，最早出现在第5卷，不过那里的“长安”并不是长安城，而是赵国的长安君。“洛阳”其次，排名81，最早出现在第2卷，三家分晋不久，洛阳附近就成了三晋与秦国的战场。

三个字的词频：

比较有趣的是“之子也”，老子英雄儿好汉。

四字字频：

鉴于“州刺史”在三字字频中频繁出现，所以频率出现比较高的各个州刺史的频率单列出来：

豫州刺史第58卷时才登场，但是雄踞第一，而豫、荆、兖、徐、扬也勾勒出了中国政治地理的热点。顺便提一句，最早登场的豫州刺史是王允，而最早登场的荆州刺史是杨震。

五字词已经没有太大的意义：

最后看看慕容家的英杰们谁的全名被提到的次数最多：

似乎是慕容彦超险胜慕容垂……慢着！慕容垂最初的名字是慕容霸，而慕容霸被提及了10次，所以慕容垂以 36 次远远胜出 （慕容缺这个全名并没有出现在《通鉴》中：））。

1999年，有了我的第一台笔记本电脑。在某一期《大众软件》或者类似刊物附送的光盘里，发现了 eDiary ，一个免费的写日记的软件。那个时候网络还不算发达，web 2.0 的概念还很遥远，所以开始把原来用纸笔写的日记，改在电脑上敲出。这是使用 eDiary 的开始。

因为那个时代自己建网站还不是很方便，所以也不能及时从网络上得到一款软件的更新，最多只能等待附送光盘的杂志或者去那时很有名气的“华军软件园”里看一看，所以最初的这个 1.0 的版本用了很久。

刚开始的时候，没有备份数据的习惯，那个时代没有网络云备份，一般电脑也不能刻录 CD ，更别说 DVD，流行的只有 2.5 寸软盘。所以早先的日记也零零散散丢过一些，还有一些存在软盘里，不知道搁在那里了。现在电脑上 eDiary 里可以追溯的最早一片日记是 2001年1月31日。这一日，日记的内容是空白。第一篇有文字的日记是 2001年2月7日，这一天除了琐事以外还记录了：

离考试还有36天。这36天要做的事情：完成国内题第二遍，北美题的第三遍（以复习总结为主），坚持背单词，做ETS的POWERPRE，还有LSAT的逻辑题，顺便读一读GMAT的阅读。记住要不断的总结。

考试说的是 GRE，从这里推断，原定的考试时间是 2001年3月15日。但是认真往后看去，真正的考试时间是 2001年4月20日 ——往后足足拖了一个月，看日记就像看镜子里的自己。

2003年11月29日, eDiary 发布了 v 2.53。然后这个作者就神秘的消失了，仿佛离开了地球。那时他已经有了一个在网易申请的免费空间，还有一个免费的域名自动转向服务。多亏了这个域名自动转向服务，从那里知道了作者最新的网站。作者有了自己的域名。不过还是沉寂。

意外的是，这个最初运行在 Windows 98 上的软件，在后来的岁月里，可以稳定的运行在 Windows 2000， Windows XP，Windows Vista 一直到 Windows 7。并且这个软件一开始就是“绿色”的，不需要安装。

所以尽管时间已经到了2011年底，却仍然使用着这个 2003 年发布的软件。好的软件是独立于时间的。

然后到了2011年12月底的某天，无意点开了软件作者的网页，居然发现他更新了！推出了 v3.0 beta 1，距离上次更新已经有了 8 年多！

一个 decade 差不多就这么过去了，外面的世界发生了翻天覆地的变化，萨达姆倒台了，拉登被杀了，苏丹分裂了，然后作者竟然更新软件了！：）

前不久在 Virtual Box 上装 Windows 8，做的第一件事就是确认 eDiary 是否能运行，总担心某一天它无法再适应新的操作系统，看到作者又重返地球，心里踏实了很多，至少不用再担心软件陈旧，无法使用。

又能安心的写日记了。“以铜为鉴，可正衣冠；以古为鉴，可知兴替；以人为鉴，可明得失。”日记里的自己就是“古”“人” 。

再次跳跃回1999年，新浪提供电子邮件服务，注册了 kzeng # sina.com 的邮件，但是没怎么用过，那时主要是是 263 ，然后是 hotmail，gmail。再跳回到2012年，终于开通了新浪的微博 —— 没写过一条微博，只有一个关注， eDiary 。

一直好奇一个关于 power function 算法的问题，直到自己动手写 power function。

在不同的编程语言里，遇到过同一个问题：(-1)^(1/3) 是多少。 很显然，在实数集里，这个表达式是有意义的，正如右图 Google Calculator 给的结果，在实数集里，它等于 -1。但是在很多编程语言的实数集运算中，这个表达式是无意义的。

譬如在 C++ 里，用 cmath library, 当你计算 (-1)^(1/3) 时，你得到输出结果是：

 -1.#IND 

在 C#里，用 System.Math, 作同样的计算，得到的输出结果：

NaN

在 R 里，不使用复数集，得到的结果也是一样的：

[1] NaN

同样的在 VBA 里，如果不 call Excel的 power function ( i.e. Application.WorksheetFunction.Power) ，而是直接使用 ^, 得到的结果仍然是 run-time error. 其它的编程语言也类似。

在复数集中，如果用 MATLAB, 得到的结果：

 0.5000 + 0.8660i

用 Maple，结果：

.5000000001+.8660254037*I

用 R 的复数集(i.e. as.complex(-1)^(1/3))，得到的结果：

[1] 0.5+0.8660254i

很显然，所有的语言的 power function 用的是同一种算法。这种算法无法得到实数解，而复数解得到是同一个数值。因为 power function 太基本，虽然有疑问，但是也没有过多的想这个问题，直到后来用 Q。

Q 的语法和 C++/C#/JAVA，或者 MATLAB/R/MAPLE 都不太一样，"^" 符号的定义和 C++/C# 相同，不是 power function。刚开始用的时候，不知道 Q  的 power function 是 xexp，觉得 power function 又不难，自己写一个吧，但是真正开始写，却又卡壳了：(int，int) 的函数好写，那 (double, double) 的呢？

拿出算法圣经《 Numerical Recipe 》（第三版），但是却发现它没有给出 power function 的算法，大概是太基础了吧，所以自己又想了一下，幸亏 Q  里的 log 和 exp 还是 log 和 exp，后来就想到用

exp(y*log(x)) = exp(log(x^y)) = x^y

来定义 power function，解决非整数的问题。这样的以来，一般的问题都解决了，但是因为用到了 log(x) , x 的值必须非负（0 的问题可以很简单的处理），所以一旦 x<0，这个算法就不适用了 —— 这时才突然的想到莫非上面的那些问题的症结正在此？实数集的问题是由于 log(x), x<0 在实数集里无定义，那复数集呢？在 MATLAB 里试了

>> exp(1/3*log(-1))
ans =
   0.5000 + 0.8660i

果然是这样的。所有的算法都依赖于 log 函数来获得 power 函数的值，这导致了上述问题在实数集无定义，而在复数集因为 log（-1)  =  3.1416i  这个默认值，导致了 0.5000 + 0.8660i  这个结果。

但是问题还没有结束。看到 log（-1），自然想到了 2*log(i), 然后自然而然的想到傅立叶变换里常用的 trick 可以解出 log（-1）的一般表达式（为了省事儿，下面用 LaTex写了）：

有了 log（-1）的通解，我们可以让 power function 获得任意 x^y, x<0 的所有解。譬如 (-1)^(1/3)简单测试一下，在 MATLAB 里，

x=2*i*(pi/2+2*(-3:3)*pi)
exp(1/3*x)

得到：

x =
        0 -34.5575i        0 -21.9911i        0 - 9.4248i        0 + 3.1416i        0 +15.7080i        0 +28.2743i        0 +40.8407i
ans =
   0.5000 + 0.8660i   0.5000 - 0.8660i  -1.0000 - 0.0000i   0.5000 + 0.8660i   0.5000 - 0.8660i  -1.0000 + 0.0000i   0.5000 + 0.8660i

结果里面包含了它的所有三个解。因为当 x<0, x^y = (-1)^y*(|x|)^y，所以只要有 (-1)^y ，就可以得到任何负数的 power function ( exp(a+bi）也可以用上面的方法转化成三角函数来解）。

上面的长篇累牍都起源于一开始的时候不知道 Q 的 power function 是 xexp，但是如果不是自己去写 power function，恐怕也没有机会搞明白 (-1)^(1/3) 这个简单的问题，俗谚云：“看人挑担不吃力，事非经过不知难。”诚哉：）

从有Samsung  Galaxy Note 的消息开始，就很期待这款 5.3 寸屏的手机（5.3 寸啊！5.3 寸啊！5.3 寸啊！），但是不出所料，这款手机要很迟才能在北美上市（甚至连会不会上市都还是一个未知数），所以当 Galaxy Note 在欧洲和亚洲发布之后，就只能关注 eBay，等待它从亚欧流入北美。

当欧洲刚刚发布这款手机的时候，eBay 也同步有人开始售卖，一开始的价格大约在 $1000 左右，并且数量很少。等了一段时间，香港的发布会过后，随着大量香港卖家的加入，价钱开始迅速下跌，很快跌倒了$800左右。作为以数字为生的人，当然会“萌”任何时间序列的漂移、扩散和跳跃（不负责的翻译 drift,  diffusion, & jump），但是不断的刷新去查看eBay页面是在令人厌恶，于是就写了一个简单的 C# 程序，定时去“挖掘”eBay 页面上的价格。

但凡提及买卖/价格，就不能不考虑风险（风险也是钱啊），特别是在 eBay 这样高风险的地方。eBay 正好在推出一个新的 beta 产品页面，在 这个页面上 eBay 已经利用自己的算法选出了风险和价格最优的产品，这样一来，挖数据就简单多了，就像附在后面的 C# codes 里显示的那样，只要 request 页面数据，把数据转化成 string，然后利用 C# 自己的 string search/index method，找到相关数据 CSS class， 读入数据即可（数据量很小，所以任何优化都不用做）。然后把数据不断存入 txt 文本。

开始测试时是每两个小时 quote 一次，后来改成一个小时 quote 一次。屏幕上跳出的数字，很像交易所的证券，所以就干脆把它用 chartSeries （R 的 quantmod package）画出来，然后就从数据里发现了很有趣的规律：

1. 首先总的价格趋势是下跌的，因为开始的高价格完全是因为 supply 的不足，再加上消费类电子产品本身的贬值以及与之相竞争的 HTC，Samsung Android 手机的推出（譬如在北美正式销售的 Galaxy SII skyRocket）

2. 其次从 Nov 11 到 Nov 14 之间，在微观结构上是两个卖家相互竞争导致价格下跌。为了拿到 eBay 页面上的产品推荐，风险相同的卖家不得不通过降价来相互竞争，但是降价的幅度一般都比较小，特别是从 Nov 13 开始，基本上就是几毛几毛的降；

3. Nov 15 左右出现了一段价格的稳定期，大约是一家放弃了价格的竞争；

4. Nov 16 其中一家的货卖完了，剩下的一家觉得自己暂时处在低风险卖家的垄断地位，所以遽然的开始升价，我们看到了一个大大的 jump；

5. 但是 jump 过后显然市场的反应冷淡，并且我们在上面说过价格的总体趋势是下跌的，如果不能及时卖出，卖家始终有一个 negative carry （也就说价格对于时间的一阶导数是负的），还有最重要的是差不多这个时候，Apple 开始 ship unlock 的 iPhone 4S ，价格大约在 $649 ，对于 Galaxy Note 的价格也是一个打击；

6. 所以过了 Nov 17 同一卖家又开始降回比原来更低的价格水平

7. 然后 Nov 18， 价格又开始大幅度的下调，因为有一个新的竞争者出现，而现在手里还攒有大量货的早期卖家基于自己的进货成本，不得不加大降价的力度， theta bleeds :)

当然数据有限，很多只是我的猜测，不过数据本身挺有趣的，不仅本身是一个很 behavior economics的测试，如果数据点足够的多，还能 fit 出一个 term structure 模型来。。。

。。。 所以到后来忘了，我只是来买手机的。。。

附 C# codes:

using System;

using System.IO;

using System.Net;

using System.Text;

using System.Text.RegularExpressions;

using System.Threading;

/// <summary>

/// Fetch eBay Price 

/// </summary>

class WebFetch

{

    static void Main(string[] args)

    {

        while (true)

        {

        // used to build entire input

        StringBuilder sb = new StringBuilder();

        // used on each read operation

        byte[] buf = new byte[8192];

        // prepare the web page we will be asking for

        HttpWebRequest request = (HttpWebRequest)

            WebRequest.Create("http://www.ebay.com/ctg/Samsung-Galaxy-Note-32GB-Black-Unlocked-Smartphone-/110509414");

        // execute the request

        HttpWebResponse response = (HttpWebResponse)request.GetResponse();

            // we will read data via the response stream

            Stream resStream = response.GetResponseStream();

            string tempString = null;

            int count = 0;

            do

            {

                // fill the buffer with data

                count = resStream.Read(buf, 0, buf.Length);

                // make sure we read some data

                if (count != 0)

                {

                    // translate from bytes to ASCII text

                    tempString = Encoding.ASCII.GetString(buf, 0, count);

                    // continue building the string

                    sb.Append(tempString);

                }

            }

            while (count > 0);

            string page = sb.ToString();

            //page = Regex.Replace(page, @"<(.|\n)*?>", "");

            // print out page source

            string targetString = "bbx2-pv";

            int first = page.IndexOf(targetString);

            string price = page.Substring(first + targetString.Length + 2, 7);

            targetString = "ship tfsp bbx2-s";

            first = page.IndexOf(targetString);

            string temp = page.Substring(first + targetString.Length + 2, 20);

            targetString = "</s";

            first = temp.IndexOf(targetString);

            string shipping = temp.Substring(0, first);

            if (shipping == "Free shipping") shipping = "$0";

            DateTime time = DateTime.Now;

            string path = @"C:\Users\Kai\Dropbox\Mis\Android App\GalaxyNote.txt";

            using (StreamWriter sw = File.AppendText(path))

            {

                sw.WriteLine("{0}\t{1}\t{2}", time.ToString(), price, shipping);

            }

            Console.WriteLine("{0}\t{1}\t{2}", time.ToString(), price, shipping);

           Thread.Sleep(1000 * 60 * 60);

        }

    }

}

大概不同的学科会造就不同的思考方式，记得若干年前和LD走在加州 Albany （伯克利附近）有些破败的马路上，突然在一片空地上看到两只狗，LD脱口而出“两只”，我也脱口而出“活的”。所以LD是学理的，我是学文的。

又好比为称赞人，学文的也许会顺口说道：“东陵虽五色，不忍值牙香。”，学理的会问：“这是说啥的？”；学文的会说：“这说的是瓜，用来称赞你的。”；学理的说：“这不是在骂人瓜娃子吗？”。。。学文的会暗自庆幸没有说：“浪乘画舸忆蟾蜍，月娥未必婵娟子”（“姮娥遂託身于月，是为蟾蜍”——张衡·《灵宪》）。

编程语言也是一样，不同之处是 epistemology 上的 （大概也只有学文科的人才会把 ontology, epistemology 这些东西挂在嘴边），就好比接触 C++ 以前从来没想到过 R 的程序可以如此这般如此这般的来写：例如以前从来没有想过 R 里的 (line by line ) debugging,  exception handling , 以及 OOP 之类的东西。一个简单的例子，譬如突然想大致画一下任意股票期权的 volatility skew。

纯粹（旧的） R 式的写法：

   1:  library('quantmod')

   2:  library('RQuantLib')

   3:  x <- getOptionChain("AAPL",Exp ="2011-07-16")

   4:  optionVal <- as.numeric(x$calls[,2])

   5:  strike <- as.numeric(x$calls[,1])

   6:  getSymbols("AAPL",src="google")

   7:  n <- length(strike)

   8:  mtrty<-businessDaysBetween(from=Sys.Date(),to=as.Date("2011-07-16"))/252

   9:  under<-Cl(AAPL)[nrow(AAPL)]

  10:  imvol <- array(NA,n)

  11:   

  12:  for ( i in 1:n){

  13:  imvol[i]=EuropeanOptionImpliedVolatility("call", optionVal[i],under,strike[i],0,0.02,mtrty,0.2)$impliedVol

  14:  }

quantmod library 提供了直接从 Google， Yahoo 等来源读取金融数据的接口，RQuantLib library 是 QuantLib 在 R 里的 port。 imvol 是一个vector， 用来储存数值算法得到 implied volatility，但是因为这只是一个粗略的计算，数据中存在着一个很大的问题：underlying 的数据和 derivative 的数据不是 synchronous 的（high frequency 数据常有的一个问题），某个 strike 的 derivative 交易的频率比较低，所以它的价格滞后，不能体现出它的真正市场价格。也是因为这个问题的存在，我们在运行这段 code 时候很快就出现错误： 

   1:  Error in EuropeanOptionImpliedVolatility.default("call", optionVal[i],  : 

   2:    root not bracketed: f[1e-007,4] -> [2.533890e+000,6.291277e+001]

因为 derivative 的价格和 underlying 的价格不是 syn 的，所以合理的 implied vol 的数值解不存在。程序因此中断，如果看 imvol 这个 vector：

[1] 1.331497 1.250151 1.481016       NA       NA       NA       NA       NA       NA
[10]       NA       NA       NA       NA       NA       NA       NA       NA       NA
[19]       NA       NA       NA       NA       NA       NA       NA       NA       NA
[28]       NA       NA       NA       NA       NA       NA       NA       NA       NA
[37]       NA       NA       NA       NA       NA       NA       NA       NA       NA
[46]       NA       NA       NA       NA       NA       NA       NA       NA       NA
[55]       NA       NA       NA       NA       NA       NA       NA       NA       NA
[64]       NA       NA       NA       NA       NA       NA       NA       NA       NA
[73]       NA       NA       NA       NA       NA

只有三个有效值，后面的值全是 initialize 的 missing values。按照旧的方法，debug 找出出现问题的地方，逐一处理，但是这样很麻烦，有时候也没有必要。所以就想到借鉴 C++ exception handling 的 try - catch方法，如果出错，我们 throw 一个 exception，然后继续其余部分的 code。R 里面也有类似的方法，所以我们可以重新写这个 loop：

   1:  for ( i in 1:n){

   2:    result<-try(EuropeanOptionImpliedVolatility("call", optionVal[i],under,strike[i],0,0.02,mtrty,0.2),silent =TRUE)

   3:    if(class(result) =="try-error") next

   4:    else {imvol[i]=EuropeanOptionImpliedVolatility("call", optionVal[i],under,strike[i],0,0.02,mtrty,0.2)$impliedVol}

   5:  }

在上面的 code 里我们作了一个简单 exception handling，先 try implied volatility 的数值解，把结果定义为一个叫做 result 的 object，如果出错，那么 result 应该属于 try-error 这个 class，所以我们作一个简单的判断，如果出错了，就直接 next 跳过。最后我们的 imvol vector 是：

[1] 1.3314972 1.2501506 1.4810157        NA        NA        NA 1.0397009        NA
 [9]        NA 1.0749333        NA 0.9936863        NA        NA 0.7165703        NA
[17]        NA        NA        NA 0.4677549        NA        NA        NA 0.5650445
[25]        NA        NA 0.4720218 0.2539875 0.3110036 0.2952723 0.2441016 0.2276886
[33] 0.2379125 0.2209396 0.2250593 0.2162720 0.2144386 0.2117865 0.2078974 0.2044520
[41] 0.2014596 0.1984621 0.1990857 0.1978392 0.1984587 0.1999078 0.2007216 0.2027609
[49] 0.2042993 0.2119021 0.2126273 0.2281725 0.2231475 0.2280039 0.1892973 0.2453698
[57] 0.2485359 0.3305246 0.3448491 0.2755123 0.2963304 0.3480423 0.3299750 0.3087564
[65] 0.3814543 0.3350278 0.3369351 0.4618993 0.3985197 0.3226662 0.4635653 0.3968701
[73] 0.3984441 0.4608958 0.3962065 0.4239336 0.4901043

尽管中途出错，但是仍然不会中断我们的 loop。

更进一步，我们可以像 C++ 一样在 R 里定义 class 和 相应的 method，从而进一步提高 code 效率（主要是 reusability ，这也是 OOP 的主要目的之一）。R 里有两种定义 class 的方式，一种旧的 S3 class，一种是新的 S4 class。 这里用 S3 作为例子，我们首先定义一个 impvol class 和它的 constructor:

   1:  impvol<-function(tick, exp="2011-07"){

   2:    getSymbols(tick,src="google",from=Sys.Date(),to=Sys.Date())

   3:    under <- as.numeric(Cl(get(tick)))

   4:    x <- getOptionChain(tick,Exp =exp)

   5:    optionVal <- as.numeric(x$calls[,2])

   6:    strike <- as.numeric(x$calls[,1])

   7:    option<-cbind(strike,optionVal)  

   8:    impvol<-list(tick=tick,

   9:              exp=exp,

  10:              under=under,

  11:              option=option)

  12:    class(impvol)<-"impvol"

  13:    return(impvol)

  14:    }

R 的 class 的定义和 C++ 还是有区别的，在 R 的 class 里面，我们只定义了它的 member variables ( tick, exp, under, option )，和它的 default constructor.  它的 method 并没有在 class definition 里 declare。定义了这样一个class 之后，我们就可以方便的获取任意stock 的数据了，譬如：

   1:  s1 <- impvol("F")

我们初始化一个关于福特的 impvol class 的object，把它定义为 s1 ，然后我们得到：

   1:  > s1

   2:  $tick

   3:  [1] "F"

   4:   

   5:  $exp

   6:  [1] "2011-07"

   7:   

   8:  $under

   9:  [1] 14.18

  10:   

  11:  $option

  12:        strike optionVal

  13:   [1,]     10      4.30

  14:   [2,]     11      3.30

  15:   [3,]     12      2.25

  16:   [4,]     13      1.35

  17:   [5,]     14      0.62

  18:   [6,]     15      0.21

  19:   [7,]     16      0.08

  20:   [8,]     17      0.03

  21:   [9,]     18      0.01

  22:  [10,]     19      0.02

  23:  [11,]     20      0.01

  24:   

  25:  attr(,"class")

  26:  [1] "impvol"

  27:  > 

上面列出了 s1 的 member variables，它包含 tick（是string obj），exp（是date obj)，under（是 numeric obj），option （是 matrix），它的属性是 impvol object。借助于 class，我们可以很轻松的获取任意的股票数据。注意我们这个 class 其实是从  list 这个 class 的 derived，所以作为 derived class，就像 c++ 的inheritance 一样，它也具有 list 的属性。

然后我们定义 method，这里主要是作图，画implied volatility，这里也和 c++ 类似，我们需要 overload 一个已存在的函数 plot():

   1:  plot.impvol<- function(x, ...){

   2:    n <- nrow(x$option)

   3:    mtrty<-businessdaysbetween(from=sys.date(),to=as.date("2011-07-16"))/252

   4:    imvol <- array(na,n)

   5:    for ( i in 1:n){

   6:    result<-try(europeanoptionimpliedvolatility("call", x$option[i,2],x$under,x$option[i,1],0,0.02,mtrty,0.2),silent =true)

   7:    if(class(result) =="try-error") next

   8:    else {imvol[i]=europeanoptionimpliedvolatility("call", x$option[i,2],x$under,x$option[i,1],0,0.02,mtrty,0.2)$impliedvol}

   9:  }

  10:  result<-cbind(x$option[,1],imvol)

  11:  plot(result,col="2",main=paste("implied vol of", x$tick),xlab="strike",ylab="implied vol",xaxt = "n")

  12:  axis(1,at=x$option[,1],labels=t)

  13:    }

因为 plot 是一个 generic function, 我们用 plot.impvol 来 overload 它，注意在我们的 plot 函数里，我们实际上先计算了它的 implied vol。检视我们的 class 的method：

   1:  > methods(class="impvol")

   2:  [1] plot.impvol

测试我们的 method:

   1:  s1 <- impvol("f")

   2:  s2 <- impvol("goog")

   3:  s3 <- impvol("aapl")

   4:  s4 <- impvol("c")

   5:  par(mfrow=c(2,2))

   6:  plot(s1)

   7:  plot(s2)

   8:  plot(s3)

   9:  plot(s4)

我们得到下面的 vol skew的 plot：

x 轴每一个 tick mark 都是一个 strike，如果那个 strike 对应的 implied vol 是 na，plot 自动就把那个数据点 drop 了。以 google 和 apple 为例，我们发现 deep ITM 的 call option 的 implied vol 有很多 missing values，并且还有异常的数值，其中的原因就是我们上面说的 liquidity 和 trading frequency的问题。最后一个 citi 的图，strike 在 $1 附近还有数据点，更能说明这个问题：这个数据源里还包括了citi并股前的数据。

OOP 大大简化了 R 的编程，所以学习 C++ 其实同时也很有助于 R 的使用，就像学点儿数学对于文科生也很有帮助一样：）

P.S. 因为编辑 CSS 出错，导致后半篇的字母全部被转化为小写 =_= 。

思考的方式，大概会受一个人正在说的语言影响，写 codes 似乎也是一样。如果是写 C++，思考的习惯大概会和 Matlab 不一样——不是说 syntax 或者是 concept，而是解决问题的抽象方式。

一个简单的题目，模拟一个 Markov Pure Jump Process (MPJP), transition matrix 是一个 8x8 的矩阵，模拟一百万次，先在 R 中写，习惯的先写一个1:8的loop 模拟一条路径，然后重复一百万次：

   1:  exp_sim=function(n){          

   2:    sim=array(0,c(n,64)) 

   3:    for (k in 1:n){ 

   4:      for(i in 1:8){ 

   5:        j=i 

   6:        t=rexp(1,rate=lam[j])           

   7:        while (t<5){ 

   8:          j=sample(1:8,1,prob=P[j,])    

   9:          t=t+rexp(1,rate=lam[j])       

  10:          }                             

  11:        sim[k,((i-1)*8+j)]=1 

  12:        } 

  13:      } 

  14:    result=array(0,c(8,8)) 

  15:    for (i in 1:8){ 

  16:      for (j in 1:8){ 

  17:      result[i,j]=mean(sim[,((i-1)*8+j)]) 

  18:      } 

  19:      } 

  20:      result 

  21:  }

结果消耗的时间可想而知，耗费了大约 340 秒才完成，实在太慢了，以为是 R 的问题，又在 matlab 里重新写过，运行大约耗费了 23 秒，正感叹 Matlab 速度快，优化作的好。。。可马上也意识到在 matlab 里写 codes 的思路竟然完全不同——根本没有去作那个 1:10^6 的 loop，而是按照习惯，一开始就生成了一个一百万行的 vector，避开了这个耗时的 loop。。。所以又参照 matlab 的 codes，重新写 R 的 codes：

   1:  exp_simv2=function(n){ 

   2:    mysample=function(vec){ 

   3:      m=length(vec) 

   4:      temp=numeric(m) 

   5:      for (i in 1:8){ 

   6:     temp=temp+(vec==i)*sample(1:8,m,prob=P[i,], replace=TRUE) 

   7:      } 

   8:      temp 

   9:      } 

  10:      

  11:      result=array(0,c(8,8)) 

  12:      for (i in 1:8) { 

  13:          t=array(0,n) 

  14:          state=i*array(1,n) 

  15:          t=t+rexp(n,rate=lam[state]) 

  16:          while (min(t)<5){ 

  17:              newstate=mysample(state) 

  18:              state=newstate*(t<=5)+state*(t>5) 

  19:              t=t+rexp(n,rate=lam[state]) 

  20:              } 

  21:              for (j in 1:8){ 

  22:                  result[i,j]=sum(state==j) 

  23:              } 

  24:          } 

  25:          result=result/n 

  26:          result 

  27:  }

像 matlab 一样思考，避开了那个一百万次的 loop，虽然写的时候有些别扭，但是运行以后发现，速度果然快了不少，效率大约提高了10倍，最后只耗时30秒左右，基本上和 matlab 的速度相仿。

习惯的力量真可怕 –_-#

P.S. 不过也发现了一个新的让人欣慰的习惯，不管是写 R 还是 matlab，已经不自觉习惯于 encapsulation，这肯定是受 C++ 的影响吧。看来习惯也是可以改的：）

喜欢 Windows Mobile 操作系统的一个重要原因就是可以使用的 GPS 软件很丰富。很长一段时间在 HD2 上使用的主要导航软件是 iGO 8。虽然 Garmin XT ，TomTom和 CoPilot 8 Live 等软件也不错并且各有长处（譬如Garmin XT可以使用 Google Search 已经下载实时的交通状况，CoPilot 有 TTS），但是从开车的易用性上讲还是和 iGO无法相比，所以从还在用 Touch Pro 的时候起，iGO 就是导航的首选，但是运行在 HD2 上的时候， iGO 8 却出现了一个不可思议的 bug：Text to Speech (TTS）无法使用。

以 HD2 的硬件配置， TTS 应该不是一个大问题，但是在 Touch Pro 上就可以流畅运行的 TTS 到了 HD2 上却非常不稳定，虽然在重启手机特别是停止 Sense UI 以后 TTS 可以正常使用一段时间，但是为了 TTS 而不断的重启手机，或者放弃 Sense ，抑或是使用一些功能比较弱的 rom 还是很不方便。

经过反复试验并结合别人的经验，总算找到一个比较完美的解决方法：使用低质量的 TTS 。因为不想覆盖我的高质量的讲 US English 的 Susan 和 Dave；我现在用的是低质量的讲 UK English 的 Kate 和 Simon，反复测试之后毫无问题，并且效果还算不错（当然如果你习惯英音的话），总算解决了这个问题。

另外， iGO 最近又推出了一款新的导航软件 iGO Primo，使用一段时间，越来越喜欢这款软件了，甚至打算用它替代 iGO 8 了。原因？看了下面的截图就明白为什么 Primo 更胜一筹了。

 
上面这幅截图是 iGO 8 用 Fast 模式算出的路线。我需要从 NJ 的 Newark 机场到 NJ 的 Jersey City。iGo 8 计算的路线硬是让我钻过 Holland Tunnel 到曼哈顿 downtown 打个 turn 再回新泽西！！！

第二幅图是 iGO Primo 给出的结果，同样是 Fast 的算路方式，同样的地图（iGO 8 的地图和 Primo 通用），它给出的路线是正确的，根本没有必要到纽约市中心掉头。

幸好这是一段熟路，所以不会被 GPS 误导；可是万一到了生路呢？相比之下，还是觉得 Primo 更省心一些，并且界面也简洁美观了很多（其实和 iGo Amigo 有些相似）。所以现在的 GPS 软件排序：1）Primo；2）iGo 8；3）Garmin XT

Updates:

忘记提 Primo 的 TTS 了，Primo 上可以用高质量的 TTS，不过用的不是 loq 的语音引擎而是 nau_v5，下载链接在这里。nau_v5 在 iGO 8 里可以显示选项，但是还是不能发生，因此 iGO 8 中最好还是用低质量的 loq 语言档。所以对于 Primo 来说，又是一个 plus。