2012年1月26日星期四

2012年1月24日星期二

車公籤錯別字

每逢農曆新年,香港都有權貴往沙田車公廟,為香港求籤。求籤本質上只是機率遊戲,求得上籤,討得好意頭,高興亦無妨,若得下籤,亦不必太過在意。只是近年不知何故,權貴一而再,再而三地求得政治意味極濃的籤文,令大眾對籤文的興趣愈來愈大。今年由新界土皇帝劉皇發所求之籤,亦不例外。籤文曰:
何為邪鬼何為神
神鬼如何兩不分
但管信邪修正外
何愁天地不知聞

解曰:凡事皆吉。
車公廟的籤文,別字甚多,它實在有另一版本,籤文與現時車公廟所印有異。網上早已有人將兩個版本的車公籤文都列出,例如「福山堂」的網頁,右上方就有 (1)、(2) 兩項,按 (1) 即可見車公廟現時籤文,按 (2) 則列出另一版本。兩年前劉皇發求得第二十七籤謂:
須防人不肖
眼前卒皆為妖
秦王徒把長城築
去禍來因自招
另版的末句則為「福去禍來因自招」。以文理而論,「福去禍來」應是正寫,但車公廟的人不知是不小心、不關心還是力有未逮,沒有發現或更正籤文錯誤。

網上有些人說此有別於車公廟的版本乃疏正所得,但在我印象中,此版本乃印於舊日籤書之中,兩個版本孰先孰後,難以斷定。況且,若此另版為疏校版本,它應該只有很少錯誤,但實情它依然謬誤處處。例如《新春秋》網誌就有人指出,以上第二十七籤之中,兩個版本皆有的「君須防人不肖」一句,「不」字應為「子」字之訛,而第二句中「鬼卒」一詞,也可能是傳鈔「走卒」之誤。今年沙田鄉事委員會求得第七十六籤,車公廟籤文曰:
人情反覆似波濤
勿謂金蘭結義高
上古守明金重利
囊空不見舊同袍
另版則將第三句改為「上古守明不重利」,同樣不通。其實八年前劉皇發亦求得此籤,當時精研術數的專欄作家王亭之已經指出,籤文第三句應為「上古守朋今重利」,這樣方可呼應籤文頭兩句。(後記:倉海君留言,指此句為「上古守名今重利」,解得更通。)

有些錯誤,更是連我等只有普通中文水準的粗人都看得出。譬如第二十六籤,另版斷曰「不信報應,惡趨善走,事急抱佛,神仙難救」,善惡二字,位置明顯掉亂了。(後記:回心一想,此處籤書可能沒斷錯,「惡趨善走」可能是「惡來善去」的意思。)我懷疑此所謂疏正了的版本,與車公廟的籤文一樣,同是從舊日手寫的籤書抄來,只是另版印於籤書之後,就再沒有改動,反而車公廟的籤文就愈抄愈錯,因此造成另版是疏正版的印象。

今日劉皇發求得的第二十九籤,另版中第三句為「但管邪修正」,從字面看,應該無錯,反而車公廟將「拒邪」二字誤植為「信邪」,就鬧了個大笑話,而將「處」(处)字誤植為「外」字,也令人覺得車公廟一方漫不經心。籤文解曰「凡事吉」,更是荒謬。對照另版
斷曰:家宅不和,占病不妙,自身不吉,出入不安,婚姻不合,求財不得。
忠奸不分,顛善倒惡,愛出風頭,非改不可。
顯然本來應是「凡事不吉」(但非凡事皆凶,原因此為中籤)。現在將不吉改為皆吉,已不止是手民之誤,而係廟方為了討好善信,而不誠不實了。

去年第十一籤「威人威威不是威」,車公廟的解籤佬竟然解為「本港經濟平穩,好壞參半,但有祖國支持。商家和工人要合力創新機」,已可見解籤佬極盡扭曲籤文勁擦鞋之能事。今年解籤佬解說謂「凡事皆吉,香港今年有些問題要想想如何解決,如果能修正就全港都好」,總算沒那麼離譜,但是將籤文中的「修正」當成更正或改正的意思,就令人有點啼笑皆非。籤文第三句,
但管拒邪修正處」(只要是抗拒邪道,修持正道的地方),
根本不是甚麼艱深的中文吧,不明白解籤佬何以解錯。記者今日訪問梁振英對籤文的看法,他謂大家對鬼神有興趣,但他稱對蒼生更有興趣,回答得相當大體,也顯示他的確比唐英年有才(最少較有口才),不過他也一樣將「修正」理解為更正或改正。

話說回來,雖然籤文開首「何為邪鬼何為神,神鬼如何兩不分」 兩句,很爆很啜核,亦雖然劉皇發於特首小圈子選舉上並非持中立立場,又雖然他很乞人憎,不過今日記者訪問劉皇發時,主動將鬼、神之說比喻特首選舉,硬要套劉皇發的話,並以此說來問唐、梁二人的看法,就未免卑鄙,兼有造新聞之嫌。姑勿論原本此比喻是由那家報館提出,但各媒體皆一樣作故仔,我對此是深感失望的。

2012年1月16日星期一

Large data sets in R

近來正在學習如何處理 large data sets,非常頭痕。

何謂 large data set,不同人有不同定義。在某些應用範疇,要起碼以 TB 來計算的才稱得上是真正的大筆資料;在另外一些範疇,幾個 GB 已算,總之 RAM 不夠裝載的,就是大筆。以我的情況來說,資料只有數百 MB,完全可以讀入記憶體,就其他人的標準而言,只是很細筆資料。然而我用的只是單機 PC,幾百 MB 到兩三 GB 的資料量,不但讀入資料須時,還佔去大量記憶體。完全讀入記憶體,並不划算。

要提高效率,首先要將資料轉為二進位檔,令讀取速度較讀取文字檔為快。然而此舉措不能省回多少讀入整筆資料的時間。由於我們每次使用資料,往往只須整筆資料的其中一小部份,因此更重要的,是要找方法 on-demand 地直接從硬碟讀取所需部份。

長遠來說,也許將資料匯入某些 SQL database 最徹底,不過暫時我只希望找到一些 quick and dirty 的方法,可以使 R/Python/Octave/C++ 都可以快速使用同一筆資料,其中,R 最令我頭痕。

舉例說,我有一個 temp.csv 檔,大約 500MB,當中每橫行均有四欄資料
n, t1, t2, s
其中 n 是整數、t1 與 t2 是一個 2008 年或以後的 UTC time strings(即 2012-01-16 00:03:59 之類),而 s 是一列長短不一,但介乎 0 與 280 的整數(例如 "1 5 236 247 248" 或 "16 135 277" 等等)。若單純將此 csv file 存為 R 的 rda format,一來每次仍須讀入整筆資料,二來,出乎意料地,讀入整筆資料的時間依然很長(在我的電腦上要大約兩分鐘),原因除了資料量較大以外,好像還因為 R 要花時間解讀 s 這個字串。因此,我首先要做的,是將 s 轉化為一個 fixed width 的欄位。280 bits 大約有 9 bytes,所以我們可用 9 個 32-bit integers 作儲存。每當 s 包含某數字 m ── 例如 m=236,我們可以將第 $\lfloor m/32\rfloor$ 個 byte 的第 $m \textrm{mod} 32$ 個 bit 設為 1。這可用 Python 輕鬆搞掂:
import struct
import time
import calendar

t2008 = calendar.timegm(time.strptime('1/1/2008', '%d/%m/%Y'))

def seconds_since_2008(time_str, fmt_str):
    if time_str == 'NULL':
        return -1
    else:    
        return calendar.timegm(time.strptime(time_str, fmt_str)) - t2008

# The last column (s) of the csv file is a list of ordered distinct
# integers between 0 and 280.  We represent it a 32-bit unsigned
# integer array of length 9.  If the list contains an integer t, we
# will flag the (t%32)-th LEAST significant bit (i.e. the (t%32)-th
# bit from the RIGHT) of the (t//32)-th byte as 1.
s_array_size = 9

s_array = [0]*s_array_size  # an array that will contain the numbers in s

infile = open('c:/temp/temp.csv', 'r')
outfile = open('c:/temp/temp.bin', 'wb')

iteration=0
for line in infile:
    if iteration==0:
        headers = line.split(",")
    else:
        data = line.split(",")
        s_array = [0]*s_array_size

        n = int(data[0])

        # The following three entries are represented as number
        # of seconds since 2008-01-01 00:00:00. NULL values are
        # represented by -1.
        t1 = seconds_since_2008(data[1], '%Y-%m-%d %H:%M:%S')
        t2 = seconds_since_2008(data[2], '%Y-%m-%d %H:%M:%S')

        # s
        s = data[3].split()
        for item in s:
            m = int(item)
            s_array[(m//32)] |= (1<<(m%32))

        outfile.write( struct.pack('<iii', n, t1, t2) )
        for m in s_array:
            outfile.write( struct.pack('<I', m) )
    iteration +=1
infile.close()
outfile.close()
問題是,如何用 R 讀回數據。

R 有三個處理大筆資料的 packages:ff, bigmemory, mmap(其實還有一個 lazy.frame,不過看來並不很成熟)。ff 功能最廣,但也最複雜;bigmemory 較側重如何用 shared memory 讀寫大量數據,亦支援 memory-mapped files ,資料格式亦似乎最有彈性。ff 和 bigmemory 均可與 biglm package 配合,用來做 linear 或 generalized linear regression。兩個套件皆曾經獲獎(亦見 bigmemory 作者的 presentation)。

mmap 的作者乃現職 quant,作者自稱動輒要應付以 TB 計的數據。以匯入資料的介面而言,三個 packages 之中,以 mmap 最直觀,但不知 mmap 能否配合 biglm 或其他套件使用。由於暫時我仍在非常前期的清理數據與基本分析階段,所以暫且先嘗試 mmap。前面我用 Python 將一個叫 temp.csv 的檔案轉為一個稱作 temp.bin 的二進位檔。以下將用 mmap 來讀它。首先定義 temp.bin 每行的格式。先前我用一個 32-bit integer 來表達 n, t1 與 t2,而 s 就用 9 個 32-bit integers 來表達,故此每行共有 12 個 32-bit integers:
> library(mmap)
> record.type <- struct(int32(), int32(), int32(), int32(), int32(), int32(), int32(), int32(), int32(), int32(), int32(), int32())
這種寫法看起來有點蠢,不過它富彈性的地方,在於每個欄位的資料種類也可不同,此處我每一欄都是一個 32-bit integer,但若第四戙是一個雙精密度的浮點數字,你可以將第四個 int32() 改成 double(),餘此類推。定義格式之後,就可以將 temp.bin 映射到一個 mmap object:
> m <- mmap(file="c:/temp/temp.bin", mode=record.type)
這樣,就可以讀寫 temp.bin 這筆數據。mmap 的原理,是有需要時,才將硬碟中所需的數據映射到記憶體之中,因此毋須讀入整筆資料。你可以用 m[i,j] 的方式提取數據:
> m[1,2]
[[1]]
[1] 83017033
此處即可見到 mmap 與 data frame 的分別。若 m 是一個 data frame,當鍵入 m[1,2] 之後,並不會有以上 "[[1]]" 那一行。事實上,若要顯示 m 的第一行的頭三個數字,mmap 與 data frame 會有資料結構上的分別。
> m[1,1:3]
[[1]]
[1] 85818

[[2]]
[1] 83017033

[[3]]
[1] 83017098
若換了是 data frame,應會有如下結果:
> m[1,1:3]
     V1       V2       V3
1 85818 83017033 83017098
出現此分別的原因,在於 mmap 將一個 table 視為 a list of columns,而每個 column 又是一個 list,因此,要讀取一個 column,所用的語法就與 data.frame 有點不同。例如要讀取第二個 column,你就應該用以下其中一種方法:
> m[][[2]]
> colnames(m) <- ('a', 'b', 'c')
> m[]$b
> m[][['b']]
若純粹要顯示第二個 column,你仍可以像顯示 data frame 般鍵入 m[,2],但是若要做頭兩戙的 linear regression,你就不能說
> library(lm)
> lm(m[,1] ~ 1 + m[,2])
Error in model.frame.default(formula = m[, 1] ~ 1 + m[, 2], drop.unused.levels = TRUE) : 
  invalid type (list) for variable 'm[, 1]'
而應該說
> library(lm)
> lm(m[][[1]] ~ 1 + m[][[2]])
用家必須小心留神。又例如要找出 m 第二戙有幾多個 0,我們不能用以下方式查詢
> length(m[,2][m[,2]==0])
> length(m$b[m$b==0])
> length(m['b'][m['b']==0])
而必須使用下列其中一種方法:
> length(m[][[2]][m[][[2]][]==0])
> length(m[]$b[m[]$b[]==0])
> length(m[][['b']][m[][['b']]==0])
這些過多的方括號很煩人,但起碼能夠 get things done。

最後,當用完一筆資料時,記得釋放數據,否則不知會否鎖死資料,令其他程式不能讀取。
> munmap(m)
mmap 套件還可配合另一個稱為 indexing 的套件(只在 R-forge,CRAN 好像還未有,而且暫時只能在 *nix 上使用),令功能更強,兼可配合資料庫使用。見套件作者夫子自道

2012年1月14日星期六

一毫子都唔值的預測

我說,蔡英文會贏。馬英九有些論點都很有道理,宇昌案表面上亦有值得調查之處,但最終決定勝負的,並非誰人的理據或能力較強(其實論施政能力,我猜也是蔡高一籌),而是誰能掌握民心。選民覺得誰是自己友,就會投誰一票。單是拋出「台灣共識」來對比「九二共識」,高下立判。

話說回來,以前我對蔡英文的印象是蠻不錯的,但近來卻愈發覺得她像陳水扁,是個巧言令色之徒,只不過更技巧更高罷了。

2012年1月11日星期三