R

计算模拟历史

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

既然有了频率,自然也就有了概率和条件概率。根据条件概率,当给出一个序列的字词后,预测下一个字词是什么,就变成了一个简单的最大似然估计问题。如果觉得这个序列太长,计算起来太麻烦,可以假设简化的马尔科夫结构,譬如假设下一个词的概率取决与之前的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  
Blog分类: 

(-1)^(1/3):从 C++ 到傅立叶变换

image 一直好奇一个关于 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 还是 logexp,后来就想到用

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

image

有了 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) 这个简单的问题,俗谚云:“看人挑担不吃力,事非经过不知难。”诚哉:)

Blog分类: 

How Blue Taiwan Is: 台湾的选举地图

Rplot

台湾2012年选举结束,国民党胜出,在 维基百科上看到了投票的结果,很好奇蓝绿投票的地理分布,于是用维基百科上给的投票数据,又从网上找来了免费的台湾 GIS 数据,在 R 里自己画了一下。这个 GIS 数据比较古老,和现在台湾的行政区划略有不同,譬如嘉义县和嘉义市、新竹县和新竹市还没有分开,而高雄市和高雄县、台南市和台南县、台东市和台东县却是分开的,还有没有新北市,只有台北县(所以上图的 Taipei 图例其实是新北市,忘记改了,呵呵)。好在边界的坐标没有大变,所以手动的作了调整,毕竟是免费的 GIS。arcGIS要更精确些(包括边界的 polygon ),但是那个要钱的:)

每个投票区域的“蓝”度由 蓝票/(绿票+蓝票) 的比例来决定,按照这个简单的算法,台湾版图上最蓝的是花莲,其次是台东和苗栗。最不“蓝”的是云林,嘉义和台南。图例的括号里的数字就是这个地区的“蓝”度。

从这个得票的地理分布和最终国民党和民进党最终 51.6% 比 45.6% 的得票率来看,民进党的竞选纲领有问题,明显的没有获得中间选民的支持,而是过于政策导向 (policy-oriented) 了。蔡英文的竞选过于侧重 base 选民,忽略了 swing 选民,失败也是在情理之中。不过有时候这确实也是一个很艰难的问题,如果不能得到 base 选民的支持,就无法从党内初选中胜出,参加最终的角逐。党内的 median voter 的 preference 如果和全部选民中的 median voter 的 preference 相差太远,不是一件好事。 同样的,就即将到来的美国大选而言,如果 Romney 能赢得初选,他对于 Obama 的威胁时最大的,但是正是因为他能和 Obama 争夺中间选民,反而使共和党内铁杆选民对他的支持减弱了,因为他会被认为不是真正的 conservative 。

下面是生成这副地图的 R  codes:

   1:  library(maptools)
   2:  library(RColorBrewer)
   3:  tw.map <- readShapePoly('taiwan/TWN_adm2.shp')
   4:  mcol=brewer.pal(9, 'Blues')
   5:   
   6:  votes = read.table('votes.csv',sep=",")
   7:  votes = votes[ ,2:5]
   8:  votes=votes[votes[,1]!="",]
   9:  votes
  10:  names(votes)<-c("loc","total","green","blue")
  11:  votes$bratio <- votes$blue/(votes$green+votes$blue)
  12:  votes$col <-floor(((votes$bratio-min(votes$bratio))/(max(votes$bratio)-min(votes$bratio)))*8)
  13:   
  14:  provname=votes[,1]
  15:  f=function(x,y) ifelse(x %in% y,votes[votes$loc==x,6],-1);
  16:  colIndex=sapply(as.character(tw.map$NAME_2),f,as.character(provname));
  17:  colIndex[c(1,15,17)] = colIndex[c(8,14,16)]
  18:  fg=mcol[colIndex+1]
  19:   
  20:  votes = votes[order(votes$bratio,decreasing = T),]
  21:  bkp=paste(as.character(votes$loc)," (",as.character(round(votes$bratio,2)),")",sep="")
  22:  plot(tw.map,col=fg)
  23:  legend("topleft", bkp,fill=mcol[votes$col+1], bty="o", cex= 0.75, pt.lwd=3)
Blog分类: 

Exception Handling & OOP in R

image大概不同的学科会造就不同的思考方式,记得若干年前和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:

image

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 出错,导致后半篇的字母全部被转化为小写 =_= 。

Free Tags: 
Blog分类: 

R 的 vectorization, again

思考的方式,大概会受一个人正在说的语言影响,写 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++ 的影响吧。看来习惯也是可以改的:)

Blog分类: 

为人性僻耽佳色

Monokai Textmate color theme天天和各种程序的编辑器打交道,有一个恶习,就是苛求编辑器的颜色搭配和字体。譬如用 Lyx 编辑 Tex 文档之前,会不断的去修改它的 preferences 里的各项颜色设定(路径在 "Documents and Settings\Username\Application Data\Lyx"),以便使用习惯的配色。

也是因为配色的缘故,我舍弃了 WinEdt,舍弃了 Tinn-R,舍弃了 Visual Studio 去使用 Sublime Text 编辑 Tex,R,或者是 VB,VC++;也是因为不支持 TT 字体(并且界面太丑),舍弃了 Texnicenter —— 当然,从功能上讲,Sublime Text 并不算强大,譬如对 LaTex 的支持远不如 WinEdt 或者 Texnicenter,以至于编译复杂的 LaTex 文档时 (譬如使用 pstricks package 画图,必须使用 LaTex –> DVIPS->PS2PDF 的时候),就不得不先在 Sublime Text 里面编辑好,然后用 WinEdt 或者 Texnicenter 来编译,因为 Sublime Text 只支持简单的 PDFLatex (顺便感慨一下世界如此之小, Sublime Text 的 LaTex 插件的开发者,竟然是我的微观经济学老师 。)

其实与其说喜欢 Sublime Text,倒不如说喜欢 Monokai 这款颜色搭配(如上图),不仅纯粹从视觉效果上来说,这样的灰底彩色搭配更容易减缓视觉疲劳(现在绝对不用白底黑字,看一个小时以后,眼睛肯定要花),并且这个色彩的名字跟我很有缘分(迷信!),所以在无可替代的时候,譬如 Lyx,譬如 Maple,就只好改这些程序的 style 设定了,好在现在一般的程序都支持客制化。为了方便,特提取了 Monokai 所包含 11 种颜色的 Hex 码如下,方便使用。

要搁在战国,我肯定是那个买椟还珠的郑人:“楚人有卖其珠于郑者,为木兰之柜,薰以桂、椒,缀以珠玉,饰以玫瑰,辑以翠羽。郑人买其椟而还其珠。”:)

#272822
#49483E
#75715E
#E6DB74
#AE81FF
#F92672
#66D9EF
#A6E22E
#FD971F
#F8F8F2
#F8F8F0

Blog分类: 

R 的 Vectorization

最先习惯的用的计量软件是 stata,然后是 vba for excel, 然后是 R,大概因为是这个学习路径,导致了在 R 中的一些不良习惯。 stata 的指令很多都是直接针对 vector 的,譬如下面的 panel data:

ID Year X Y
CHN 1990 45 36
CHN 1991 23 62
CHN 1992 55 90
JPN 1990 79 56
JPN 1991 60 72
JPN 1992 78 100

如果需要只需 CHN 的数据,在 stata 里只要:

keep if ID==”CHN”

就可以了。而 VBA 中因为 excel 表格的直观,一般是用 for  或者 do loop:

Sub chn()

n = Range("A50000").End(xlUp).Row
for i = n to 2 Step -1
if Cells(i,"A")!="CHN" then Cells(i,"A").entiredRow.delete
end if
next i

End Sub

因为 R 中的 if statement 不能直接针对 vector, 所以就习惯像 VBA 一样用 loop 来整理数据,但是因为 excel 中最大的数据量不过 6万五千多行 256 列,数据量并不大,所以用 loop 虽然慢,但是因为数据量的限制,还算能忍受,但是在 R 中,如果有上十万的 observations,再用 loop 就极其的痛苦了,虽然R中的loop 的效率要比 VBA 高(尽管不如 S),虽然刚换了 Duo Quad 的 CPU,但是仍然能算到死机 (虽然可以 remote access 学校的 cluster 上去,但是浪费学校的资源也不好),所以迫使自己重新学习 R 中的 vectorization,其实比如上面的例子,在 R 中用 vector 作也简单(假设上面 data frame 的名称是dat):

chn<-dat[grep("CHN",dat$ID),]

也挺简单的。所以程序还是一个摸索熟练的过程 ( 虽然我还没有发现如何像 stata 一样便捷的在 R 中的 panel data 里创建 lag variable,并且还在等待大牛人们完善 R 中的 panel data package,譬如 fixed-effect negative binomial model,但是 R 的灵活性,特别是在编程和模拟上,还是超出其他软件很多,特别是 R 是 免费的!)

Free Tags: 
Blog分类: 

开源软件的苦与乐

开源软件的一点感想,不过不是说 Drupal (当然 Drupal 作为开源软件的一种也有同样的问题),而是 R。原本只做数学分析模型,不做统计计量模型,但是无数的 reviewer 都说: 证据呢?好吧,作计量,但是简单的计量回归已经无法满足“人民群众日益增长的复杂要求”,所以做数值模拟,做蒙特卡洛,用到 R,虽然以前没怎么用过,但是编程语言都有相通之处,譬如和 Excel 里的 VBA 相比, loop 或者 flow control 的结构基本相同(除了不用写 end if 或者 next i,略省些力气),再就是把 Cells (i, j) 换成 Mat[i, j] (总之只要把一个 spread sheet 看成一个 matrix,基本的思路都是一样,当然 R 可以直接进行 vector 运算,譬如 elseif,这个要方便很多),做模拟和矩阵运算,R 的效率都很不错,code 也简单,譬如重复 probit 估测一千次,R 里面不过几行代码,不到一分钟的运算,这点比以前常用的 Stata 要好很多 —— 这是开源软件灵活性的好处。

但是作一些简单的操作,却又凸现了开源软件的协调问题 (coordination problem)。举个最简单的例子,在 panel data 里生成 lagged variable。这样的操作在 Stata 里只有一行指令:

gen lag_var=l.var

(当然你要先用 xt 设置好 panel data)。但是在 R 里面,尽管 time-series data 相关 package 里有定义好的 lagged variable 函数,但是在 panel data 相关的 package 里面却没有,因此如果想完成这样一个简单的操作,必须写一个复杂的函数,把需要 lag 的 variable 从panel data 的 data.frame  里按照unit抽出来,设置为 time-series data,然后生成 lagged variable, 再用 cbind 合并回去,然后用 by 按照不同的 unit,再用 rbind 把各个 unit 的生成的 row 合并回去,然后 do.call 。

如果作 dynlm 和作 plm package的人能合作一下,统一做一个根据数据类型生成 lagged variable 的函数岂不是很方便( Stata 里的 laggged variable 指令是同时适用于 time-series data 和 panel-data 的)——也许这就是商业软件的好处吧,至少相互协调的很好。

当然这还不是最麻烦的地方。最麻烦的问题在于开源软件给了你灵活性的同时,也要求你很多东西要 DIY,虽然有很多现成的 package ,但是这些 package 都是大家自愿贡献的,所以不可能要求他非常的全面,当你需要一些 package 里没有的东西的时候,就不得不自己动手了,譬如想对 panel data 作一下 fixed effects poisson/negative binomial regression, plm package 里没有,其他任何的 package 中也没有找到,自己写? 虽然这在 Stata 里只是 xtpoisson / xtnbreg 这么简单,虽然这在技术上也不算太高深(Wooldridge 那本关于 panel data的教科书里有算法),但是从头开始写起恐怕没有一个星期是弄不完,并且问题在于我只是一个打酱油的 –,- … …

结论:开源软件固然好,商业软件离不了:)

Blog分类: