顯示包含「不學無術」標籤的文章。顯示所有文章
顯示包含「不學無術」標籤的文章。顯示所有文章

2014年12月23日星期二

The tetrahedron puzzle

大家都知道甚麼叫正四面體吧?


將一個正四面體切成兩半,交由別人重新拼合成正四面體,有幾難?


研究人工神經網絡的著名學者 Prof. Geoffrey Hinton 說,原來這是很難的,許多麻省理工的教授都搞唔!以下是其演講 "What's wrong with convolutional nets?"連結)的一部份 transcript(從 13:50 開始),相當惹笑。
Here's something you wouldn't believe. I take a simple object, like a tetrahedron. A tetrahedron is a pyramid with a triangular base. And I slice it with a plane. That's a flat thing ... I slice it with a plane. Uhm, so I get two pieces. And then I take an intelligent person, and I give him the two pieces, and I say "OK, make a tetrahedron." And I make sure he knows what a tetrahedron is. And he can't do it. Now I clearly ... now you don't believe that, presumably. It's just two pieces. Surely you can put it together to make a tetrahedron.

I present the experiment today to MIT professors. I got one sample thirty years ago, that's a professor called Carl Hewitt. He is very smart. I gave him the two pieces and ... he looked at them for a long time ... he didn't really play with them but looked at them for a long time, and determined whether he could write down a proof whether it's impossible. [laughter] So, his time to solve the puzzle is infinite. [laughter] OK?

Today I've been doing the experiment to MIT professors, and the number of minutes they take to solve the problem is about the number of years they've been in MIT. [laughter] Roughly speaking, it's definitely very ... the length of time is very positively correlated to how long they have been in MIT. [laughter]

So now I'm going to show you this puzzle, because it's extraodinary. It can't be so hard. It's completely trivial, and there's an obvious way to solve it that people don't figure out. They figure out in a few minutes, but, OK, here's the two pieces ...

(之後 Prof. Hinton 用了幾分鐘示範一些人解這個謎題時是多麼困難)

Why is this puzzle almost impossible? Why is it so hard? 'Cos it's a two-piece jigsaw puzzle. [laughter]

And the MIT professors, and ... oh, incidentally, I tried this on a Google vice president, just to reassure the MIT professors. I gave these two pieces to the Google vice president and said, "This is really a hard task. Can you make a tetrahedron?"
欲知 Google 的 VP 同 MIT professors 相比,結果如何,請看演講片段。

2014年3月21日星期五

網文偶讀之百年一遇

  1. R.J. Oosterbann, Frequency and regression analysis of hydrologic data (pdf) 
  2. Floods: Recurrence intervals and 100-year floods; USGS 
  3. 看香港如何防“水浸” 排水干渠系统200年一遇;北京晚報,2012年07月22日 
  4. 香港渠務署網頁
兩年前七月廿一日的北京暴雨,據說令全中國非常感動,各人為了賑災,有錢出錢,有妹捐妹。內地大部份的報章都稱此暴雨為「百年一遇」,局部地方的單日降雨量甚至是「五百年一遇」云云,惹來網民嘲弄,謂「在我短短的一生裡,百年一遇的洪水見過10次,千年一遇的地震見過2次,唯獨四年一遇的全民大選還沒遇見過」。這句譏諷話獲內地網民瘋傳之餘,報章亦有人指摘「官方」以「百年一遇」之說來推卸責任,例如《金羊網》就有篇專欄文章說:
筆者認為,說白了,無論是推出“40年一遇”、“60年一遇”、“61年一遇” 之說,還是拋出“百年一遇”、“接近五百年一遇”之論,無非就是想告訴大家︰這是老天爺的錯,要怪就怪老天爺去吧。

「百年一遇」的定義
「百年一遇」指的其實是自然現象的重現期 (return period)。「重現期」是大學專科(譬如水利工程學)術語,但它牽涉的統計概念只屬高中程度。簡單來說,若某自然現象 A,於某單位時間(譬如一日)內發生的機會率為 p,那麼其重現期(譬如以日數計),就是 1/p 個時間單位。

詳細一點來說,當我們談及重現期的時候,背後假設了 A 於每個單位時間內發生與否,皆為獨立事件,亦即是假設 A 是 i.i.d. Bernoulli(p),而所謂 A 的重現期,就是距離下一次觀察到 A 所需的平均時間。譬如我們以日為單位,現在是第一日的開首,而下一次要第 $T$ 日才觀察到 A,那麼 A 的重現期,以日數計就是 $E(T)$。很明顯:
  • 第一日就觀察到 A(亦即 T=1)的機會率為 $p$;
  • 第二日才觀察到 A(亦即 T=2)的機會率為 $qp$($q=1-p$);
  • 第三日才觀察到 A(亦即 T=3)的機會率為 $q^2p$;
  • 第 n 日才觀察到 A(亦即 T=n)的機會率為 $q^{n-1}p$;
  • 故此 $E(T) =  p + 2qp + 3q^2p + ... + nq^{n-1}p + ...$(練習)。若讀者念過大一統計學的話,當然知道 T 依隨的實乃參數為 p 的 geometric distribution。
「重現期」既是科學術語,有關現象的定義自然要準確。假若我說「雨是 7.5 日一遇」,是沒意思的。換作「天文台沙田氣象站的雨,是 7.5 日一遇」好一點,因為指明了事件發生的地方。換成「天文台沙田氣象站錄得 40 毫米以上的單日降雨量,乃 7.5 日一遇的事件」就更好,因為說明了談話者所關心的是多大的雨。

網民嘲諷中共「官方」稱這次事件乃百年一遇,然而我找過一些內地大報,如《人民日報》、《光明日報》與《中國日報》,若不計轉載內容,它們均無提過該暴雨為百年一遇,例如《中國日報》一篇報道就只是說這次是 the most devastating downpour in the Chinese capital for 61 years。中國氣象局、北京市氣象局及其轄下機構好像也沒有以「百年一遇」來形容該次暴雨。所以,公道一點地說,當時的「百年一遇」之說,並未得到中共中央或氣象部門認可。

一些北京報紙,如《北京青年報》、《新京報》或《北京晨報》,倒報道過「百年一遇」一說,而且它們的消息來源皆為「北京市人民政府防汛抗旱指揮部副指揮潘安君」。據《北京晨報》報道:
市气候中心昨天公布的数据显示,比起21日历史罕见的大暴雨,自1951年有完整气象记录以来的京城历史上,单日降雨量排名“亚军”和“季军”的降雨日分别为1952年7月21日和1954年8月9日,但老天爷降下的雨水比起“冠军”却远远不及。
   昨天,市防汛抗旱指挥部副指挥、市水务局副局长潘安君通报了此次特大暴雨的四个“历史罕见”:降雨总量之多“历史罕见”,全市平均降雨量170毫米,城区平均降雨量215毫米,为新中国成立以来最大一次降雨过程,房山、城近郊区、平谷和顺义平均雨量均在200毫米以上,降雨量在100毫米以上的面积占本市总面积的86%以上;强降雨历时之长“历史罕见”,一直持续近16小时;局部雨强之大“历史罕见”,全市最大点房山区河北镇为460毫米,接近五百年一遇,城区最大点石景山模式口328毫米,达到百年一遇,……
從上述報道,可見:
  1. 潘安君只說過而北京當日的全市單日平均降雨量屬「新中國」成立以來之冠,並沒說此降雨量是「百年一遇」。
  2. 然而他確有聲稱「石景山模式口」當日雨量是「百年一遇」,更說「房山區河北鎮」的降雨量為「五百年一遇」。

如何估計重現期
前面提到一個現象的重現期,是距離它下次發生的平均時間。要估計這個期望值,從文獻所見,常用方法似有三種。

第一種是以簡單的相對頻數來計算事件發生的概率 p,從而推算重現期 $E(T)$。舉例說,若根據過往記錄,某 B 區於過去十年(n = 3652 日)當中,有 m = 13 日的單日降雨量超過 100 毫米,那麼,對於「B 區的單日降雨量達100 毫米以上」這個事件,我們所估計的 p 就是 $\frac{m}{n} = \frac{13}{3652} = 0.0036$,而此事件的重現期為 $E(T) = \frac1p = 281$ 日,或者 0.77 年。因此,若今天 B 區剛巧錄得 100 毫米以上的單日降雨量,我們可以稱這場雨為「0.77 年一遇」。

第二種方法是用有序統計 (order statistics)。假設我們將前例中 3652 個單日降雨量由大至小排列,發現當中排第 r = 13 的單日降雨量為 102 毫米,那麼,對於「B 區的單日降雨量達 102 毫米以上」這個事件,我們估計 $p = \frac{r}{n+1} = \frac{13}{3653}$(留意此處與第一種方法不同,分母為 n+1),而重現期為 $E(T) = \frac1p$,同樣大概是「0.77 年一遇」。

第三種方法是 distribution fitting,亦即是將過去的觀測記錄模配到理論上的統計分佈之中。若模配的是長期平均值的分佈的話(例如每年平均單日降雨量),由於大數定律,用 normal distribution 就可以了;若模配的是最大值的分佈的話(例如每年最高單日降雨量), 較常用的則有 Gumbel distribution 與 exponential distribution 兩種。如何將觀察資料模配到理論分佈之中,或者如何檢定模配的合適度,是大學統計學的內容,此處不贅。

假設我們手頭上只有某地過往 N 年的單日降雨量記錄,可以證明(練習),無論今日我們觀察到該地的單日降雨量為何,若用上述頭兩種方法,而且 Y>N+1,我們是絕無可能得出今日這場雨是「Y 年一遇」這個結論的。《北京晨報》說北京自 1951 年起有完整氣象記錄。其實中共建政後,北京於 1951 年建立第一所氣象站,當時是否已有「完整氣象記錄」,不得而知,但就算有,從 1951 年至 2012 年才不過 61 年,因此,若說北京雨災當日,「石景山模式口」的降雨量是「百年一遇」,甚至「房山區河北鎮」的降雨量為「五百年一遇」,是絕無根據的。

若用的是第三種方法,由於所用的分佈是伸延至無限大的連續曲線,無論如何大的降雨量,我們也的確可以估計出它的重現期。問題是,若降雨量像前述的北京暴雨那樣,完全高於過往最高記錄,則如此推算出來的重現期,實際上等於是用外插法 (extrapolation) 取得。用 normal values 來估計 extreme value distribution,或者借用近年新興的說法,以白天鵝來估計黑天鵝 (black swan),是任何受過良好基本統計訓練的大學生,都不會視為可靠的做法。

這不是說水利工程師不會用只得幾十或一百年的資料來估計「二百年一遇」事件的規模。事實上,以我們香港的渠務署為例,它對於市區排水幹渠系統(urban drainage trunk system,「幹渠」即是排駁大型集水區、防洪標準最高的渠道)與鄉郊的防洪堤堰 (flood protection bund) 的設計標準,是必須能應付「二百年一遇」的水浸。只是,工程師除了要對特定的重現期估計洪水或降雨規模之外,還要估計有關的信賴區間 (confidence interval)。更重要的,是渠務署會不時監察排水系統的表現,而不是紙上談兵。「重現期」是用來幫助工程師估計排水系統所需的設計參數,並非供官員於洪水雨災發生後用作遁辭。

2014年3月19日星期三

勿因蟲廢言(三)

基於莫名其妙的原因,昨日忽然有較多讀者閱讀本網誌上週的文章《勿因蟲廢言(二)》。即使此處並非人氣網站,文章刊登與瀏覽量高峰之間的時差,也從未達四日那麼長。無論原因為何,重讀舊文一遍之後,為免讀者誤會,且容我澄清一下。

是次 HKUPOP 「特首民望調查」,梁振英的平均得分低於 50。梁粉認為這個結果有誤導之嫌,並提出幾項論據,其中較注目的一項,指受訪者所打的 998 個分數中,有 91 個屬零分這個「極端分數」,拉低了平均分數。

我不認為這個論據成立,原因如下:

  1. 有時計算平均數,我們的確希望排除「極端」數值(離群值,outliers)的影響。然而所謂「極端」,指的不單是樣本數字的大小,還有該樣本數字出現的概率。換句話說,我們要排除的,是極端罕見而且又有極端數值的樣本的影響。舉例說,設想某個介乎 0 至 100 的數字 X,若出現 X=0 的概率為 0.49505,出現 X=100 的概率也是 0.49505,而出現 X=1, X=2, ... 以至 X=99 的概率,各為 0.0001。那麼,儘管 0 與 100 是 X 的可能數值中極小與極大的兩個,但它們才是正常的數字,絕不應排除。
  2. 我們並無先驗 (a priori) 理由認為「特首民望調查」中的 0 與 100 分份屬罕見。故此它們應視為正常樣本數字,而非離群值。
  3. 若論語意的話,由於問卷只解釋了 0, 50 與 100 三個分數的意義,所以受訪者揀這三個分數,比揀其他分數更正常。因此,0 和 100 更不應視為離群值。
  4. 何況是次 998 個分數中,接近 9%(91 個)是零分。有如此多零分,我們更有理由相信它們並非離群值,而且受訪者所打分數其實屬多峰分布 (multimodal distribution),而非類似鐘形曲線 (bell curve) 的單峰分布 (unimodal distribution)。
  5. 退一步說,就假設分數呈單峰分布,而 0 與 100 是離群值好了。由於有效分數只可介乎 0 與 100 之間,故此屬有界分布 (bounded domain distribution)。一般而言,相比無界分布,在有界分布中,離群值對平均數不會有嚴重影響。
  6. 實際上,今次樣本,在不加權的情況下,就算我們略去高低各一成樣本數字,而計算截尾平均 (trimmed mean),結果仍與簡單的樣本平均數無大分別,兼且仍低於 50 分。

由此可見,梁粉圍繞「極端數值」的論據並不成立。

梁粉的批評,有另一個毛病,就是無視了平均分應該加權計算這個事實。

由於社會因素(例如從日間到黃昏,留在家中的成年人,應以家庭主婦和長者佔多),用固網電話做民意調查,必然有某類人會較少受訪 (under-represented),而另一類人則訪問過多 (over-represented)。為了修正這個誤差,我們就要以受訪者的統計特徵(例如年齡、性別等等)來分組,將各組人佔受訪者總數比例,與人口普查所得的比例比較。若發現某組人於訪問中所佔比例太小,這一組樣本就要加重權數。這就是做民意調查時,受訪者經常會被問及年齡、性別等等個人資料的原因。

且不說今次調查中,打低於 50 分的受訪者,其實比打高於 50 分的還要多。就算給超過 50 分的受訪者真過半,若他們是 over-represented,一經加權修正,出來的結果也可能比 50 分低。

梁振英所得的平均分,加權後為 47.5(也就是 HKUPOP 公布的數字),未過 50,但此分數只是單點估計 (point estimate)。既是抽樣調查,就必然有統計誤差。如何分析這個統計誤差?有一篇由署名「陳電鋸」的大學博士生所寫網誌,就用了正統的 bootstrap method 來分析。結果顯示,即使考慮了統計誤差之後,我們仍有信心說,梁振英所得的加權平均分,明顯低於 50。

「輔仁媒體」刊登了另一篇反駁梁粉的文章。文章作者處理權數的方法並不妥當(概念上,權數應施加在樣本數字的出現頻率上,而非樣本數字本身;若果只是計算加權平均數這個單點估計,兩者效果並無分別,但若要做其他分析,就必須搞清楚概念),不過該文章也有有道理和有趣的地方,不妨一讀。

梁粉的另一項批評,是 HKUPOP 使用平均分不妥。這方面我是同意的,但我認為更根本的問題,並非在於「平均分」,而是曖昧不明的「分數」本身。

HKUPOP 的問卷只指明了 0-100 分這把量尺中其中三點(0、50 及 100)的意義,其他分數刻度的意思曖昧不明。儘管根據陳電鋸另一篇文章分析,多年來從 HKUPOP 調查計算所得的「特首民望」平均分數,與中大亞太研究所獨立調查所得大致吻合(不過他並無說清楚,其文章所指有很強相關的,究竟是兩個研究所的民望指數,還是兩組指數的升跌),然而,這頂多表示市民於大概同一時段內,心中所用的量尺大體一致,卻不能為升跌的幅度賦予任何意義。我們無從判斷,心中量尺刻度的意義,是否隨時而變。同是從 30 跌至 20,若於不同時間發生,意義是否一樣?即使現刻,我們亦無從判斷,能令民望從 47.5 升至 57.5 的政績,是否與能令民望從 47.5 跌至 37.5 的過失等量齊觀。

若然我們不能說明民望指數變幅的實際含義的話,「民望指數大跌」與「民望大跌」就難言是同一回事。若有不明就裏的受眾搞錯了,就算兩所研究機構並非存心誤導,伯仁也是為他們而死。

2014年3月15日星期六

勿因蟲廢言(二)

上週香港發生了兩單非常嚴重的事件,其一是警方竟然闖入區議會,抬走正在正常開會的區議員;其二是政府以一個極其荒謬的理由,要求明明應該只受《電訊條例》監管的「香港電視網絡」流動電視業務,亦要受《廣播條例》監管,令港視於正式開台之前,必須取得本地免費或收費電視節目服務牌照。

儘管香港淪為中國殖民地以後,情況一直變壞,市民也一直說「低處未算低」,然而今次這兩宗事件,卻是香港正式告別英治年代政治格局的里程碑。舊日講求公務員政治中立,和講究規章制度的精神,今日已完完全全淪為「親疏有別」、「法律因人而異」的人治局面。兩天前於電視上聽見中國總理李克強說「中國是法治國家,不論是誰,不論職位高低,法律面前人人平等」,我聽了只得冷笑。黨大於法的國家一邊自詡法治,一邊摧殘其殖民地原有的法治制度;自稱「人民民主專政」,卻擁有永遠執政黨的野蠻人,大言不慚,指點有半熟代議政制的文明人何謂民主,何謂法治,何謂普選……甚麼叫「匪夷所思」,香港歷史就是最佳註腳。

比起上述兩宗嚴重事件,今日我於「主場新聞」看到的鬧劇就顯得非常次要,然而此事上我總算有一些實質的東西可以說,故此不妨一談。鬧劇的起因,又是「梁粉」批鬥鍾庭耀。詳情請閱以下兩篇立場迥異的網絡文章:
公說公有道,婆說婆有理?
「梁粉」批評如下:
依據港大最新的民調,以100分為滿分,特首僅獲47.5平均分,當然就被評為不合格了。然而,只要打開原始資料,就會發現998個評分者中,原來有多達615人、即逾6成人均給予特首50或以上的合格分數,其中更有29人給予100分;僅有383人給予50以下的評分。那麼,為何特首的評分又會不合格呢?最大的問題在於有91人個受訪者給予0分,就是這些極端評分,令特首的平均分大幅度拉低。
「主場」及香港大學民意研究計劃研究經理李偉健則反駁:
 評論指有91個0分樣本「拉低」平均分,沒有提到29個100分樣本同時會「拉高」平均分。港大民意計劃研究經理李偉健向《主場新聞》解釋,民望調查詢問受訪者給予官員0分至100分的評分,相信受訪者誠實回答,無論樣本是0分或是100分,都應納入計算,除非是101分,在數值範圍之外才會剔走。

李偉健強調,歷來民望調查同樣沿用這方法,公布按評分計算算術平均值(Arithmetic Mean),「沒有篩走特別低、特別高的評分。」
開門見山。我認為「梁粉」的批評有其道理,但其為己方所作辯解,一樣有問題。另一邊廂,「港大民研」的統計方法也有毛病。

Lies, damned lies, and 梁粉's statistics
統計數字不會說謊,它有的只是統計偏差。說謊的,是運用它的人。"Lies, damned lies, and statistics" 這句名言,就是用來諷刺那些蓄意運用統計數字來製造假像的人。前述「梁粉」的批評,正好拿來作「統計語言偽術」的最佳範例。

從「特首民望調查」所得到的 998 個有效評分,平均分為 47.4(「港大民研」公布數字為 47.5,略有不同,這是因為他們按受訪者的統計特徵作加權平均),低於 50,但實際上 998 個分數當中,有 615 個為 50 分以上……至此,梁粉都沒有說錯。然而,他們沒說的是:
998 個分數當中,也有 663 個為 50 分以下。
感覺混淆嗎?或者這樣說吧,998 個分數當中,有 383 個低於 50 分,280 個等於 50 分,335 個高於 50 分。分數的分布如下:

從 0 到 100,共有一百零一個整數,而 50 正好居中。梁粉試圖以「50 分以上」這個標準來描繪一個梁振英有超過六成人支持的景象,可是據他們的邏輯,我們同樣可以說,以「50 分以下」這個標準來判斷的話,有超過六成人(而且這個「超過六成」的人數比起梁粉的「超過六成」更多)反對梁振英!

我不明白一眾梁粉何以如此介懷 47.5 這個只略低於 50 的數字。若是選舉的話,兩三個百分點也許是勝負關鍵,可是像印象分這種雖非玄學,卻也「不算精密科學」的東西,47.5 和 50,實在沒有分別。換了我是梁振英,看到如此數字,高興還來不及呢。

離群值與平均數
梁粉指出,998 個分數當中,有 91 個是 0 分,這些極端評分拉低了整體的平均數。這是正確的。「主場」卻反駁梁粉,說他們沒提及樣本當中亦有 29 個 100 分,會有拉高平均分的相反效果,也同樣正確,亦再一次顯示梁粉玩弄輸打贏要的統計語言偽術。

然而,撇除梁粉的拙劣技倆不談,若樣本中可能有不少「離群值」(outliers) 的話,到底我們應該如何估計統計母體的平均數?

港大民研的李偉健指「無論樣本是0分或是100分,都應納入計算」。就一般統計調查來說,這是過時的做法(但此處有一個 catch,要押後談)。現代統計學認為「穩陣」(robust) 的做法,本網誌之前的書評其實已經提過,就是利用截尾平均 (trimmed mean),也就是先截去最高和最低的 5-10% 數據,然後才計算平均數。

可是我們幾乎可以斷言,在「特首民望調查」中,無論用普通的算術平均,抑或用截尾平均,都不會有大分別。原因是一般來說,離群值最有殺傷力的情況,是母體數字本身為「無界」(unbounded) 的時候。是項調查當中,有效的評分本身有界(只可介乎零至一百),離群值的影響通常不會太壞,故此梁粉的批評,抓不到統計學的重點。

實際上,若截去今次樣本當中,高低各一成的數據的話,得出來(未經加權)的截尾平均為 48.1,與樣本平均數 47.4 相去不遠。

尺度不同,分數如何換算?
這倒不是說「特首民望調查」無問題。印象中,港大民研所做的民意調查,大部份(例如立法會選舉的選前調查和 exit polls)都很紮實。然而此項「特首民望調查」,卻非常礙眼。我很想問鍾庭耀一句:
How on earth is this rating meaningful? 
單單叫受訪者為梁振英打個分數,已經很有問題。問卷只提過零分(「絕對唔支持」)、五十分(「一半半」)與一百分(「絕對支持」)的意義,中間的尺度 (scale),人人卻不同細分。你我各給六十分,意思未必相同。你的分數如何換算成我的,完全木宰羊。現時港大民研的做法,實際上假設了所有人的評分尺度均一。由此引起的模型風險 (model risk),無法評估。舉個例說,若你看到梁振英的「民望指數」比上月高,你可能以為他真的愈來愈受市民歡迎,但實情可能是他的民望無變,只是今個月的受訪者的評分尺度較寬鬆,對無甚特別感覺的官員,也傾向打一個高分而已。

就算是奧運體操項目,評分有較多稍為客觀的細項憑依(動作要求、難度、時限等等),仍不時惹人爭議,各人對特首表現的評分尺度,又怎可能大致一樣?

不知尺度,何論變化?
好了,就假設香港有一個平均的評分尺度吧。套用經濟語言來說,就當人人都用一個一致「市場評分尺度」好了,但為何我們可以計算平均分?平均數並不一定是有意義的。一半人給零分,另一半給一百分,借用時下流行語來說,是社會撕裂的狀況;所有人都打五十分,卻更似人人認命。兩種情況截然不同,平均分都是五十分,那麼五十分究竟是甚麼意思?

以上例子當然太極端,極端到與雷鼎鳴對堅尼系數的批評如出一轍。假若港大民研只是拿這個平均分來判斷粗略民情的話,上一段的批評是不適用的。問題是,港大民研對待這個平均數時,彷彿其精密數值或它幾個百分點的變化,有甚麼微言大義似的。然而,即使香港有一個「市場評分尺度」,我們仍不知道這個尺度是甚麼樣子。同樣是跌十分,從一百跌至九十分,是否跟六十跌至五十,或十跌至零同樣大鑊?木宰羊。五十分所代表的「一半半」,和「及格」是同樣意思嗎?木宰羊。不及格的話,甚麼分數才算民怨沸騰,很想梁振英辭職?木宰羊。

不知背後的評分尺度的話,再精密的數字都是沒用的。弄得好像很精密,反而令人誤以為該數字很科學,其細微變化很有意義。

離群值真是離群值嗎?
前面說過,以普通的算術平均來估計母體平均數,乃過時做法。諷刺的是:
  • 對「特首民望調查」來說,由於整把由零至一百分的量尺中,只有零、五十及一百有清晰意義,所以這三個分數,比其他分數可靠。
  • 故此,吊詭地,0 和 100 兩個離群值,反而不應剔除。
  • 結果梁粉針對離群值的批評,意外地不適用。
  • 若硬要計算平均數,普通的算術平均,此處亦反而比截尾平均更恰當。
然而這不表示港大民研的做法正確。正正因為他們採用了語意不明的尺度,才造成這許多奇怪狀況。

結語一:less is more
如前述,港大民研的民意調查,一般都很紮實,但這項「特首民望調查」,用粵語來說的話,真係「畀位人插」。"Less is more" 這句說話聽來陳套,但此處適用。奉勸 Robert Chung,還是乾脆將問卷問題改成簡簡單單的「你想唔想梁振英繼續執政」之類好了,不要再搞那些懶細緻的評分吧。

結語二:廢話去死,自由萬歲
最後且談文字,不談統計。梁粉謂:
港大民意研究計劃的民調早陣子引起連串質疑,未知是否有見及此,今次港大再度公布特首評分時,民意網站已出現所謂的「原始資料」,雖然相關檔案的格式要以特定軟件才能打開,但內裡所刊載的正正是評分分布數字。
這不是廢話嗎?有甚麼檔案是任何軟件都可以打開的呢?何況所謂「特定軟件」和檔案格式,也不過是統計佬慣用的 SPSS 與它的 sav 格式吧。不想付鈔的朋友,可用免費的自由軟件 R 打開有關檔案。

相關網頁
伸延閱讀
  • 電鋸,你玩統計,統計玩你:「問題根本不在於 0 和 100 等等 outliers ,而是佔人口比重較多的組群對梁振英評分較低。

2014年3月5日星期三

勿因蟲廢言

張 has a point
恒基地產副主席李家傑,昨日大力批評港大民調的「特首民望調查」,共產黨報章及土共聞聲而動,紛紛幫忙叫罵,霎時將今年「驚蟄」提早了兩日。

這些「小爬蟲」(小爬蟲黃定光語)的評論,大部份都不值一哂,然而「梁粉」張志剛的批評,卻有一點道理。按AM730報道,他說:
 「如果問聽日投票,投梁振英定係唔投梁振英嘅話,你就只可以將個結果解讀為『如果聽日投票,梁振英可以攞到幾多票』,唔可以當作係施政支持度,而且如果講投票選舉嘅話,一定要列明對手。」
香港人根本就不能(或未能) 選舉行政長官,因此以這種「假想得票率」(鍾庭耀語)來估計任何其他數字,都必然有所歪曲。況且,如張志剛所說,就算要估計選舉得票率,也應該指明對手。若將「梁振英對余若薇」換成「梁振英對曾德成」,結果就恐怕有天淵之別。

鍾 knows what he is talking about
鍾庭耀如何擬出現在的問卷?按香港大學民意研究計劃網站的說明:
「港督或特首民望」調查一直沿用的提問方式為「而家想請你用0-100分評價你對港督/特首某某某既支持程度,0分代表絕對唔支持,100分代表絕對支持,50分代表一半半,你會俾幾多分港督/特首某某某呢?」及「假設明天選舉特首,而你又有權投票,你會唔會選某某某做特首?」。
而按鍾庭耀於2003-6-12的解釋:
筆者在12年前引入0至100的評分標準,是經過一番考究和深思。該等基準,在民主社會時有用之,但不算普遍。西方社會慣用百分率顯示政治領袖的「認許率」或「支持率」,提問方法大抵可以分成兩大類: 

(1) 你是否認許某某處理其崗位(如總統、首相等)的工作? 

(2) 如果明天進行選舉,你會否投票支持某某? 

一般而論,答案只分正反兩面,再加「不知道」。所謂「正反兩面」,可能是採用二分或四等法。二分法的答案可以是「認許」、「不認許」;「滿意」、「不滿意」;「是」、「否」;「會」、「不會」等;而四等法的答案可以有「很滿意」、「頗滿意」、「頗不滿」、「很不滿」;「肯定會」、「可能會」;「可能不會」、「肯定不會」等。 

西方社會涉及0至100評分調查的對象,大多數不是選舉出來的政治人物,或者只是醞釀參選的人物。提問者通常會要求被訪者想像一個攝氏溫度計,0度表示冰冷,100度表示沸騰,50度表示不冷不熱,然後要求被訪者打分形容其對某人的感覺。顯而易見,這種方法毋須被訪者想像任何選舉處境,或與其他人物比較。 

不過,這種提問方式得出來的數字,往往不能直接轉化成為選民基礎的強弱,對預測選舉結果的作用不大。民主社會習慣以選票較高下,選舉期間如是,民意角力也如是。因此,西方社會的民意指標主要採用二分法來量度領導人的「假想得票率」,或「工作表現認可率」,類溫度計的數字則屬次要。西方社會稱民意調查為「opinion poll」,就蘊藏了民意公決的概念。 

香港的民主進程只屬初階,除了從地區直選出來的立法議員外,用「假設得票」的問題來量度民望,對被訪者而言,似乎沒有太大意義。因此,筆者在12年前便採納了0至100分的標準,作為主流測試方法。不過,時至今日,「問責」觀念高唱入雲,巿民的民主訴求逐漸熾熱,民望指標亦要相應調整。更何況,部份官員已學會偷換概念,把40多分的支持度評分說成是獲得四成多的支持比率,還說這樣的民望可以媲美世界各國。這些言論,或許能淡化巿民的不滿,但是本身有違科學精神,對民意不公,對社會無益。
由此可見,鍾庭耀本人是很清楚其問卷是有甚麼潛在流弊的,他也是經過深思熟慮後,才將問卷設計成現在這個樣子。你可以不同意他的見解(我就不同意),但他足足十年前已公開解釋其問卷設計,而這是可供公開學術辯論的。將他的調查形容為對民意的「操控」(李家傑語),絕不公允。

Approval ratings
話雖如此,我也一向不大滿意鍾的「特首民望調查」(他將支持率與反率對的差額,用一個自創兼懶學術的詞語「民望淨值」來包裝,就更叫人眼冤了)。要採用他口中所謂「二分法」抑或「類溫度計的數字」(0至100分),只屬旁支末節。最根本的問題,是他究竟想量度甚麼?

既然這是一項「特首民望調查」,想量度的,當然是市民有多想行政長官下台了。港大問卷的兩條問題當中,張志剛已指出第二條問題的弊端。至於第一條問題,它用的是甚麼字眼?
「而家想請你用0-100分評價你對港督/特首某某某既支持程度,0分代表絕對唔支持,100分代表絕對支持,50分代表一半半,你會俾幾多分港督/特首某某某呢?」
這個「支持」 是甚麼意思?是滿意的意思?是我會以實際行動幫助他(例如加入「愛字頭」的鋤奸隊,或盡可能在共幹面前替他美言幾句)?還是不管我有多(不)滿意,我都贊成他繼續坐在他的位子上?用語言來表達思想,很難完全精確,但是「支持」這個詞語,意思卻未免太過模糊。

如鍾庭耀所說,有些外國的民意調查,會用上'approve'(鍾庭耀譯作「認許」)這個字眼。例如美國的市場調查一哥蓋洛普,就有每週追蹤調查,問:
Do you approve or disapprove of the way Barack Obama is handling his job as president?
亦有一些調查,用的字眼為「滿意」。例如英國的老字號市場研究公司 Ipsos MORI 的問卷,就會問:
How satisfied or dissatisfied are you with the way David Cameron is running the country /doing his job as Prime Minister?
留意問卷問的並非國民有多滿意首相的政績,而是問他們有幾滿意首相的施政。這兩者差之毫釐,謬以千里。問政績的話,遇著時勢艱難,官員難有作為,政績很差,但國民可能諒解,滿意其施政,認為官員已盡力扭轉乾坤,只是形勢比人強。就好像第二次世界大戰後期的北非戰場,德國的隆美爾元帥頻頻吃敗仗,但不少軍事史家反而對這個階段的他,有很高評價一樣。

我們不想問卷問題太長,但受訪者又未必了解「政績」與「施政」的分別,故此我還是覺得蓋洛普的問法較好。然而直接挪用又有點不妥,皆因 'approve' 一詞,隱含說話者有權干預的意思,但現今體制下,特首並非民選,立法會又有非民選的功能組別,兼投票時有分組點票,要說港人 'approve' 特首施政,真是無從談起。其實要量度官員的民望,我們可以單純了解市民有幾想換人。要避開張志剛所指出的困難(假設明日有選舉,亦要指明對手),我們其實不必扯甚麼「明天有(假)選舉」,只需簡簡單單地將問題更改如下:
你贊唔贊成梁振英繼續施政?
當然,文字還可以斟酌(後記:「你想唔想梁振英繼續執政」或者更佳),但我想說的是,鍾庭耀問卷的毛病,更深層反映的其實是語文問題。

2014年3月1日星期六

冰滑橋搖 ── 中學物理課本錯誤兩則

我念中學那個年代,各校物理課多採用外國教科書(實際上本地出版的,當時好像只有中大楊綱凱教授所編那本)。這些課本對現實世界的現象,時有錯誤解釋,我也是後來看其他書籍,才發現原來這些解釋是錯的。針對這些錯誤的批評,有一些流傳得比較快和廣,令相關錯誤很快從教科書消失。例如我那個年代的課本,許多都用 Bernoulli's principle 來解釋「用風筒打斜吹向乒乓球,但乒乓球在空中固定不倒」這個實驗 (Youtube),但後來人們知道這是誤解,正確的解釋應該是 Coandă effect 加上 Magnus effect,而這兩種效應的成因,其實超出了中學物理學的範圍。故此,儘管這個乒乓球實驗很有趣,但根本就不應拿來做課堂例子。(芸芸物理學原則中,Bernoulli's principle 也可能是最常被中學課本濫用的一個,但這是另話。)


然而也有一些錯誤,是流傳至今的。日前上網,發現兩個我念中學時所讀到的物理解釋原來又是錯的時候,心裏依然很驚訝。

為何人們可以溜冰?

Robert Rosenberg, Why is Ice Slippery? Physics Today, December 2005, pp.50-55.

水(及一些其他物質)受壓時,熔點會降低,這個現象稱為 pressure melting。以往課本多以此來解釋何以人們可以溜冰:溜冰的冰刀壓著冰面,令冰塊溶化,潤滑冰面,減低摩擦力。也有一些課本以摩擦生熱 (frictional heating) 來解釋冰面水分。

可惜,根據不同物理學家所做實驗,人們發現 pressure melting 這個解釋完全是錯的 ── 冰刀的壓力的確會降低冰塊的熔點,但憑常人的重量,施加在冰面的壓力,最多只可降低冰塊的熔點三、四度。當冰面接近零度時,以 pressure melting 來解釋冰面滑溜,乍聽還有理,但問題是,實際上冰面容許人們溜冰的溫度,可以低至零下三十五度,此時冰刀的壓力,卻不足以令冰面溶解。

其實物理學界對於冰面為何滑溜到可以溜冰,至今仍未有定論。暫時最令人信服的解釋,並非摩擦生熱,而是一種稱為 premelting 的原理(勿與前述的 pressure melting 混淆),其大意是,結晶體即使低於熔點,結晶面仍可存在準液態薄膜。連學者都仍未了解的事情,其真正的物理解釋當然也超出中學範疇。實際上,單是要證實冰面上有薄薄的水膜(而不論此水膜是否容許溜冰的關鍵),竟然也要出動到核磁共振儀器!

Tacoma Narrows Bridge 塌橋事故

K. Yusuf Billah and Robert H. Scanian, Resonance, Tacoma Narrows bridge failure, and undergraduate physics textbooks, American Journal of Physics, 59(2):118-124, 1991.

當一個可以振動的系統遇到週期性的外力,而施力頻率接近系統的自然振動頻率時,就會發生大幅振動現象,稱為共振 (resonance)。我那個年代的中學教科書(甚至大學教科書),不少都以共振來說明一些塌橋事故,或者以共振來解釋士兵走過吊橋時,為何不能步操,必須打亂步伐。事實上,過往確曾有大量士兵過橋時,發生塌橋事故,其中傷亡最嚴重的,當數 1850 年法國 Angers 懸索橋斷裂意外。當時有一營士兵正在過橋,儘管他們事先已知道要打亂步伐,但士兵為了在大橋搖幌時平衡身體,反而不自覺地令步伐一致。當大橋從懸索開始斷裂時,有 483 名士兵在橋上,最後有 226 人死亡。

1850 年連攝影都未普及,更遑論攝錄機。中學教科書引述得最多的,並非 Angers 大橋斷裂事故,而是有聲有畫,1940 年末在美國發生的 Tacoma Narrows Bridge 塌橋事故。




這宗極富戲劇性的經典意外中,無人傷亡,唯一死者,是被主人遺留在車上,名為 Tubby 的犬隻。(Tubby 的狗種為 Cocker Spaniel。香港將此狗種譯作「曲架犬」,巧合地與此塌橋事故語帶相關。)好些教科書都將這次意外形容為風與橋的共振,卻無解釋當時風力何以帶週期性。後來有人提出,並不是風力帶週期性,而是風吹過橋體時,造成稱為「Kármán 渦流街」(Kármán vortex street) 的週期現象。
然而有趣的是,且不說這個由渦流引起的震動 (vortex-induced vibration) 技術上算不算共振,原來早已有人驗證過這種震動並不會危害大橋。話說當 Tacoma 峽灣橋通車後不久,人們已發現大橋會隨風上下起伏。由於搖動幅度頗大,有人甚至報稱橋面隆起至看不見前車,而時人將此橋戲稱為 "Galloping Gertie"。親歷其境的過客當中,有一位華盛頓大學教授 Burt Farquharson。他事後利用風洞實驗,證明這種上下搖幌是安全的。而且,出事當日,大橋乃左右扭動至斷裂,搖幌模式與平日不同。

那究竟是甚麼機制導至大橋倒塌?這要到 1991 年,前述普林斯頓大學與約翰霍金斯大學兩位工程學者的論文出版後,才有公認的解答。事件的元凶,是一種稱為 aeroelastic flutter 的現象。大體來說,當風吹過橋體時,會造成一些小旋渦,令橋面扭動。許多時,這種擺動就像鐘擺般,會慢慢停下來(即 "damping",「阻尼」)。然而,若風速、橋面闊度與扭動的頻率皆符合特定條件的話,橋面扭動時形成的攻角,又會造成更多旋渦,而這些旋渦又令橋面扭動得更厲害,直至大橋斷裂。留意這種現象並非共振(事件中,施加在橋體的外力並無週期性),而是渦流與橋身擺動形成的正反饋 (positive feedback)。



究竟我們在中學(或大學時),還學了幾多錯誤的物理解釋?

2012年11月23日星期五

2012 Peking University grad school entrance exam (Higher Algebra)

From math.SE:

1) Let $\xi_1, \xi_2, \ldots, \xi_n$ be all the roots of a polynomial $g(x)$ (defined over the complex field) with rational coefficients. Suppose $f(x)$ is an arbitrary polynomial with rational coefficients. Is $\prod_{i=1}^n f(\xi_i)$ necessarily a rational number? Prove your assertion.

2) Show that the following determinant is nonzero:

$$
\left|\begin{matrix}
1 & 2 & 3 & \ldots & \ldots & \ldots & 2010 & 2011\\
2^2 & 3^2 & 4^2 & \ldots & \ldots & \ldots & 2011^2 & \color{red}{2012^2}\\
3^3 & 4^3 & 5^3 & \ldots & \ldots & \ldots & \color{red}{2012^3} & 2012^3\\
\vdots\\(k-1)^{k-1} & k^{k-1} & (k+1)^{k-1} & \ldots & 2011^{k-1} & \color{red}{2012^{k-1}} & \ldots & 2012^{k-1}\\
k^k & (k+1)^k & (k+2)^k & \ldots & \color{red}{2012^k} & 2012^k & \ldots & 2012^k\\
\vdots\\
2010^{2010} & 2011^{2010} & \color{red}{2012^{2010}} & \ldots & \ldots & \ldots & \ldots & 2012^{2010}\\
2011^{2011} & \color{red}{2012^{2011}} & \ldots & \ldots & \ldots & \ldots & \ldots & 2012^{2011}\\
\end{matrix} \right|.
$$

3) An order $n$ matrix $A$ has exactly one nonzero entry on each row and each column, whose value is either $1$ or $-1$. Show that $A^k=I$ for some positive integer $k$.

4) Define the "product" of two $n\times n$ matrices $A=(a_{ij}), B=(b_{ij})$ as

$$
A\circ B = \left(\begin{matrix}
a_{11}b_{11} & a_{12}b_{12} & \ldots & a_{1n}b_{1n}\\
a_{21}b_{21} & a_{22}b_{22} & \ldots & a_{2n}b_{2n}\\
\vdots\\a_{n1}b_{n1} & a_{n2}b_{n2} & \ldots & a_{nn}b_{nn}
\end{matrix}\right).
$$
If $A$ and $B$ are positive definite, show that $\mathrm{rank}\,A\circ B\le(\mathrm{rank}\,A)(\mathrm{rank}\,B)$.

5) Suppose $f_1,f_2,\ldots,f_{2012}$ are 2012 different linear transformations on a vector space $V$. Does there exist a vector $\alpha\in V$ such that $f_1(\alpha),f_2(\alpha),\ldots,f_{2012}(\alpha)$ are mutually different? Prove your assertion.

6) Suppose $A$ and $B$ are positive definite matrices of order $n$. Prove that they can be simultaneously congruence-diagonalized by some invertible matrix $T$.

7) Let $f(\alpha,\beta)=g_1(\alpha)g_2(\beta)$ be a symmetric bilinear form on a Euclidean vector space $V$ over a field $P$. Show that there exist some linear function $h(x)$ and some $k\in P$ such that $f(\alpha,\beta)=kh(\alpha)h(\beta)$.

8) For any $n$-dimensional Euclidean vector space $V$, show that there are at most $n+1$ vectors such that the angle between any two of them is obtuse.

11) Given that the following linear transormation is a rotation in $\mathbb{R}^3$. Find the rotation axis and the angle of rotation.
$$
\left(\begin{matrix}x'\\ y'\\ z'\end{matrix}\right)
=\left(\begin{matrix}
\frac{11}{15}&\frac{ 4}{15}&\frac{ 2}{ 3}\\
\frac{ 4}{15}&\frac{13}{15}&-\frac{1}{ 3}\\
-\frac{2}{3}&\frac{1}{3}&\frac{2}{ 3}.
\end{matrix}\right)
\left(\begin{matrix}x\\ y\\ z\end{matrix}\right)
$$

2012年9月20日星期四

Julia 初體驗の立法會選舉勝算 DIY

我慣用 C++ 與 Matlab/Octave,偶爾也用 Python 及 R。近年眼見不少有趣語言出現,我尤為喜歡 Ruby, D(兩者其實都不算新), ScalaChapel,可是除了 Ruby 我依然計劃會抽時間學之外,其餘都只得三十秒熱度。直至最近偶然碰到 Julia,覺得應該先與她(實在無法說「它」呀)把臂同遊。

C++ 的發明人 Bjarne Stroustrup 將 C++ 形容為 "a general purpose programming language with a bias towards systems programming"。依此說法,Julia 大概就是 "a general purpose programming language with a bias towards scientific computing" 了。若不計 Julia 語法上對矩陣的直接支援,我想一般 programmers 應該不會將她聯想成 domain-specific language 吧。

我昨天才下載 Julia,約會了一日,她給我的第一印象,是她絕對有潛力成為 Matlab 殺手或 R 殺手。無論是語法的簡潔程度、data structures 的數量、語法上對 functional programming, generic programming 及 parallel computing 的支援,抑或程式的執行速度,Julia 都明顯超越對手,網上不少 Matlab 與 R 用家亦對她頗為讚賞。她也借用了 Ruby, Python, Matlab 與 C/C++ 語法當中一些優良部份,我學習時倍感親切。

Julia 今年一月才出 1.0 版本,算係有女初長成,距離亭亭玉立還有一段日子,現在仍有不少未成熟的地方,不過已經夠足我做練習用。前文提過,計算立法會選舉各競選名單的勝算,可用多項式分佈的常態逼近當成投票的分佈。利用蒙地卡羅模擬實驗,就可以計算出各名單的勝算。文友電鋸於選舉前已經做過類似的計算,此處只是當成我第一次的 Julia 編程練習。

先說明計算細節。設 $\mathbf{p} = (p_1,\ldots,p_n)^\top = $ 民調所得各名單的支持度 ($\sum_i p_i = 1$),而 $n$ 是樣本數,譬如 NOW 新聞台於九月七日報道(調查時窗為九月二至六日)的結果為

$$\mathbf{p}=\frac{1}{101} (4, 7, 5, 7, 1, 2, 16, 9, 3, 3, 8, 4, 1, 8, 1, 12)^\top,\ n=503.$$
(由於 NOW 新聞台四捨五入,以上向量內各數字的總和為 101 而非 100。)我們的做法,是模擬多次投票實驗。每次實驗,均由 n=503 位選民,每人隨機投一張名單一票。投票的概率由 $\mathbf{p}$ 決定。換句話說,每一名選民都會有 $p_1=\frac4{101}$ 的機會投票予第一張名單、$p_2=\frac7{101}$ 的機會予第二張名單,餘此類推。當 503 人都投完票,就可以按比例代表制查出名單上各人是否當選,查核完畢,就完成了一次實驗。重覆同樣實驗許多次 ──  譬如 1,000,000 次,就完成了整個模擬過程。整個過程當中,若排在七號名單第二順位的余若薇當選了300,000 次,她當選的概率就估計為 300,000/1,000,000 = 0.3,其他人的當選概率也用同一方式估計。

n = 503 位選民每人按 $\mathbf{p}$ 的概率來投票,即是說各名單得票 $(X_1,\ldots,X_m)$ (今屆新界西有 m=16 張競選名單)的分佈為 $\textrm{Multinomial}(n, \mathbf{p})$。多項式分佈的常態逼近公式,可參考本網誌前文,當中用到的正交矩陣 Q 的構作方法,則見我另一篇網誌。總括來說,設

$$
\begin{eqnarray}
v &=& \frac12 \left(
\begin{bmatrix}0\\ \vdots\\0\\1\end{bmatrix} -
\begin{bmatrix}\sqrt{p_1}\\ \vdots\\\sqrt{p_m}\end{bmatrix}
\right),\\
M &=&
\begin{bmatrix}\sqrt{\frac{p_1}n}\\ &\ddots\\&&\sqrt{\frac{p_m}n}\end{bmatrix}
\left(I_m - 2\frac{vv^\top}{\|v\|^2}\right).
\end{eqnarray}
$$若每次投票實驗,我們皆能夠生成 $m-1$ 個服從標準常態分佈的隨機數字 $Z_1, \ldots, Z_{m-1}$(重申,m 是競選名單數目,n 為投票人數),則每次實驗各名單的得票率可模擬為:

$$
\begin{bmatrix}\frac{X_1}n\\ \vdots\\ \frac{X_m}n\end{bmatrix}\approx \mathbf{p} + M\begin{bmatrix} Z_1\\ \vdots\\ Z_{m-1}\\ 0\end{bmatrix}.
$$
生成得票率之後,將它正規化為 $m$ 倍:
$$\mathbf{f} = (f_1, \ldots, f_m)^\top = m\left(\frac{X_1}n, \ldots, \frac{X_m}n\right)^\top,$$
之後就可以點票。正規化後的黑爾數額為 1(原本黑爾數額為 1/m,乘以 m 倍就變成 1)。每個 $f_k$ 都是一個實數,其整數部份 $\lfloor f_k\rfloor$ 代表第 k 張名單因超過黑爾數額而取得的議席數目,小數部份 $r_k = f_k - \lfloor f_k\rfloor$ 代表餘額。比較各餘額的大小,就知道餘下 $m - \sum_k \lfloor f_k\rfloor$ 個議席落入誰家。最後程式如下:
 

我和 Julia 還不是很熟。上面有些迴圈,相信可以用比較 functional programming 的方式寫得精簡一些。

按上述 NOW 新聞台的調查結果,由 Julia 所計算,十六張名單中各候選人(是候選人,不是候選名單)的當選概率依次為:
  1. 郭家麒 (1.0)
  2. 譚耀宗 (1.0)
  3. 李卓人 (0.999993)
  4. 田北辰 (0.99991)
  5. 李永達 (0.99954)
  6. 梁耀忠 (0.999528)
  7. 陳偉業 (0.996816)
  8. 麥美娟 (0.996756)
  9. 陳樹英 (0.918242)
  10. 余若薇 (0.821678)
  11. 陳恒鑌 (0.716208)
  12. 梁志祥 (0.716037)
  13. 陳一華 (0.338804)
  14. 何君堯 (0.338146)
  15. 龍瑞卿 (0.066209)
  16. 曾健成 (0.065039)
  17. 譚駿賢 (0.004338)
  18. 麥業成 (0.002589)
  19. 陳強 (0.002561)
  20. 張慧晶 (0.000565)
... 等等。一百萬次實驗,在我的電腦上需時共廿九秒。至於何謂有「極高機會」當選,何謂「較高機會」、「機會均等」、「機會較低」,就見仁見智了。

若只想知道那九位候選人有最高機會當選,那其實毋須搞甚麼模擬實驗,因為以上當選概率的高低名次,基本上與民調得出的支持度的高低排列相同。因此,模擬實驗結果中勝算最高的九位候選人,就是民調結果中支持度最高那九位。這類勝算計算的目的,其實並非要找出最有機會當選的是誰,而是要反映這些勝算較高的候選人,與其他候選人的差距。

這個方法也有它的毛病,當中最嚴重的,是它沒有考慮政黨配票的情形。以上例來說,民建聯梁志祥與陳恒鑌的支持度一直低企,約四、五個巴仙左右,可是到了選舉日,他們的得票率比起民調結果大幅上升,相反,譚耀宗的得票率就比民調結果低許多(變成約 8%)。要考慮配票,就要考慮各政黨互相鬥法,結果可能要計算隨機博奕下的 Nash equilibrium。除了較複雜之外,均衡點是否存在,是否唯一,亦造成很大的技術困難。

2012年9月17日星期一

Fate/Legco: 機器學習與選舉工程

有兩夫婦與一位鄰居遇到海難,三人各乘一艘只可載一人的充氣快艇逃生。大海茫茫,暴風雨迫近。他們估計所餘燃料只夠二人逃走,究竟誰要犧牲?最後三人相持不下,全數沒頂,這是誰的責任?
前言
剛過去的立法會選舉,民主黨於新界西全軍覆沒。另一邊廂,公民黨的兩人競選名單成為新界西票王,獲 72185 票之多,扣除首席所需的基本票額之後,餘額卻不夠令名單中排次位的余若薇連任。有人指摘公民黨策略錯誤,既分薄了民主黨票源,又平白浪費選票。

本文將運用機器學習 (machine learning) 技巧,說明根據過往經驗,今屆各次滾動民意調查當中,即使民主黨的支持度處於最高峯 (14%) 的時候,其形勢也有很大隱憂。多數時候,民調結果更顯示民主黨有全滅危險。民主黨要保住一席的話,只要棄車保帥就可以,因此它遭全滅,完全是它自身策略錯誤所至。

過往最常見的選舉分析方法,是從滾動調查所得的政黨支持度,按多項式分佈的常態逼近 (normal approximation to multinomial distribution) 來計算勝機。一般期刊文章,若要以票站調查結果來模擬比例代表制之下的選舉結果,就是用多項式分佈。這樣做的最大好處,在於可以隨時更新評估結果。原則上,此方法更可以納入整體形勢,而不是對每張名單都只按它自己的支持度來計算勝機。然而,遠在選舉提名期尚未結束,政黨連有甚麼對手都未清楚之時,這個方法就不適用。

本文將提出一種靜態的研究方法,它只依賴過往的選舉結果,並不需要(但也可以用)最新的民意調查結果。它無法動態地考慮通盤選舉形勢,但它可以於選舉提名期間,就用來評估分拆或分併競選名單的風險。作者將估計勝算的問題聯繫到數據科學中的「分類問題」(classification problem)。解決問題的方法稱為 logistic regression,以今日的標準來看,只是統計學與數據科學 (data science) 的初等技巧,但足以應付眼前的新界西個案。

基本須知
本文假設讀者知道何謂比例代表制,亦明白本地立法會選舉比例代表制所用的最大餘額法當中,黑爾數額 (Hare quota) 的意思。從某個角度而言,黑爾數額是一張擁有多名候選人的競選名單當中,每名候選人所能消耗的票數上限。與此相關的是特羅普數額 (Droop quota),它是保證一名候選人當選所需的票數下限。近年網上多了談及 Droop quota 的文章(例一例二),可是為求完整,以下亦會稍作說明。

以今年新界西選區為例,各黨派一起競爭 9 個議席。問題:無論其他名單得票若干,民主黨的李永達名單,至少要有幾多得票率,方可保證任何情況下當選?

答案是 $\frac1{10}$。理由:設 9 位當選者的得票率為 $p_1, p_2, ..., p_9$。若李永達取得 $\frac1{10}$ 的票仍不夠當選,那即是說每名當選者的得票率 $p_i$,必然較 $\frac1{10}$ 為高,故此

$$
\begin{align*}
&\phantom{=}李永達與\ 9\ 位當選者的總得票率\\
&= \frac1{10} + p_1 + \ldots + p_9\\
&> \frac1{10} + \underbrace{\frac1{10} + \ldots + \frac1{10}}_{9 個}
\quad(因為每個\ p_i\ 都大過 \frac1{10})\\
&= 1.
\end{align*}
$$
亦即是說,李永達與 9 位當選者的得票率總和,竟然大於 1,這顯然是不可能的。同一道理,若選區有 n 個議席,只要排於名單首位的某君取得
$$d_1(n)=\frac1{n+1}$$ 的票,不管其他對手取得幾多票,此君亦鐵定當選,否則就會出現「部份候選人的得票率總和竟然大於 1」這件邏輯上不可能發生的事。上面這個 $d_1$,就是一般文章所講的 Droop quota,也就是保證排於名單首位的候選人取得一席所需的安全線

機器學習:得票率與選舉結果的關係
Droop quota 只是一條安全線,它是當選的充份條件 (sufficient condition) 而非必要條件 (necessary condition)。候選人越過它,則鐵定當選,但不越過也有可能當選。以今屆立法會選舉新界西的結果為例(下表,綠色當選,紅色落選),十六張名單爭奪 n=9 個議席,安全線為 $d_1(9) =  \frac1{9+1} = 0.1$。從表中可見,九名當選者之中,只有郭家麒一人越過安全線(即是得票率 $\ge 0.1$)。



現實中,經常有候選人未達安全線但仍當選。當然,得票率愈接近安全線,候選人就愈篤定,反之風險愈高,如何量度這個風險?

我們可以憑歷屆選舉結果,評估得票率與當選機會之間的關係,但首先要就選區議席數目調整得票率。舉例說,對一個有 9 個議席的選區(安全線 $d_1(9)=0.1$),若一張名單的得票率為 v = 0.08,應該有不錯的勝算。然而,換了是一個只得 4 個議席的選區(安全線 $d_1(4) = 0.2$),0.08 的得票率就未免離安全線太遠。因此,若一張名單的得票率為 v 而選區議席數目為 n,我們會考慮以下這個「規範得票率」:

$$
\textrm{normalized proportion of votes}\ v' = \frac{v}{d_1(n)} =(n+1)v.
$$v' 愈接近 1,即代表該選舉名單愈接近安全線。搜集從 1998 年起四屆立法會的選舉數據,可將歷屆各參選名單的規範得票率與選舉結果圖列如下。(由於 $v'\ge1$ 就必然當選,故下圖並不包括 $v'\ge 1$ 的例子。)
圖一:1998-2008 年立法會選舉各候選名單的規範得票率
從圖中可見,過往當 v' 介乎 0.5 至 0.8 的時候,當選與落選的例子混雜。這引出以下問題:若某候選名單估計會有某規範得票率 v'(譬如 v'=0.75),那麼,根據過去經驗,我們認為該名單的首席候選人會否當選?

這就是機器學習理論所謂的「分類問題」。換句話來說,我們希望憑規範得票率 v' 的數值,就可以將候選名單歸入「當選」或「落選」其中一類。當然,從上圖可知,當選與落選名單的規範得票率之間,並無一條清晰界線。因此,機器學習理論的做法,是將過往的選舉結果稱合 (fit) 到一條概率曲線之上,如下圖。譬如圖中顯示,若某名單的 v'=0.75,我們即相信該名單有 0.83 的機會當選。顧名思義,一般分類問題的目的,是為了替目標對像分類,不過我們這裏的着眼點,是評估一張競選名單的當選機會,而不是賭它會不會當選,因此我們不會說「v'=0.75 時,我們相信該名單會當選」,而只會說它有 0.83 的勝算
圖二:用概率曲線稱合 1998-2008 年的規範得票率
如何稱合概率曲線?答案是用 logistic regression。Logistic regression 是統計學與機器學習理論之中的初等方法。詳情相信大部份讀者沒有興趣知道,故此處不贅。上圖是作者利用歷屆立法會選舉結果,從統計軟件 R 計算所得,稱合的勝算曲線為:
$$f(v') = \frac{1}{1+\exp(7.899974-12.664143v')}. $$
方法利弊
從規範得票率的歷史數據來學習當選機會,有利有弊。如文首所說,它最大的好處,是只依靠歷史數據,完全不需理會現狀,甚至遠在上屆選舉剛剛結束,今屆選舉提名期尚未開始之時,即可找出稱合的勝算曲線。

然而這也是它最大的毛病。我們固然可以於民意調查開始之後,將政黨的最新支持度當成得票率,再從勝算曲線讀出當選機會,這樣做卻忽略了對手的選情。這好比打麻雀只盯着自己摸回來的牌,而不顧其他三家一樣,並非善用情報的方法。

此外,由於這個方法只依賴歷史數據,但各政黨的排陣、支持度與選區議席數目皆因時而異,故此歷史經驗未必適用。事實上,今屆有選區(新界西與新界東)議席增加為九個,就前所未見,不過以前已經試過一區有八個議席,所以這還不算是大問題。最嚴重的,是現況完全超出過去經驗的時候,稱合所得的曲線就會變成廢物。例如根據 1998-2008 年的數據,若我們想稱合一張名單會取得第二席的勝算曲線,就會遇着這樣的情況,不過本文重點是取得第一席的概率,故詳情此處不贅。

換句話說,儘管 logistic regression 背後有堅實的統計學基礎,但歷屆選舉並非處於同樣條件的統計事件,因此我們實際上並非用 logistic regression 來解決一個統計問題,而僅僅是把它當成一種插值法,用來將曲線稱合到數據當中。這也算是沒辦法之中的辦法,卻也是讀者必須留意的一點。

圖表詮釋
今屆民主黨於新界西用李永達、陳樹英兩張名單參選。去屆該黨於新界西取得 23.25% 的選票。今屆安全線為 $d_1(9)=0.1$,若該黨得票率與上屆相若,則已超過兩張名單越過安全線所需,所以,毋須任何麻煩的分析,只要民主黨能夠令支持者比較平均地投票給兩張名單,就幾乎可以保證全取兩席。

只不過,「得票率與上屆相若」 這個假設,於民意調查甫一開始時,就已經站不住腳。



以上七次滾動調查結果,以最尾一次最利民主黨。暫且假設這個結果可信。根據此結果,民主黨的支持度有 14%。原則上同黨各名單的支持度是不可相加的,原因是選民可能「投人」而不是「投黨」,不投李永達,也未必會改投陳樹英,更何況有些人本來投的就是遊離票。只不過以舊鑑新總要有個標準,因此本文仍假設選民乃按黨派投票,而同一黨的支持票是可以於各名單之間轉移。這樣,我們應如何看待民主黨這 14% 的支持度?

記民主黨的得票率為 $v_\max$,並設李永達名單的得票率為 v,因此它的規範得票率為 $\frac{v}{d_1(9)} = 10v$,而陳樹英名單的規範得票率為 $10(v_\max-v)$。故此,根據歷屆數據,李、陳兩張名單的勝算分別為 $f(10v)$ 與 $f\left(10(v_\max-v)\right)$。若民主黨新界西的總得票率真的如調查所得,即 $v_\max=0.14$,那麼李、陳兩者的勝算曲線有如下圖,其中藍線代表李永達,黑色虛曲線代表陳樹英,橫軸為李永達的得票率(不論是李永達的曲線還是陳樹英的曲線)。例如九月七日的民調顯示李永達的支持度為 0.08,圖中顯示了這個情況之下,陳樹英的勝算為黑色曲虛線的高度,亦即 0.43,而李永達的勝算為藍色曲線的高度,即 0.9。

從這幅圖也可推出兩人皆勝或兩人皆輸的機會率。當李永達的支持度為 0.08,陳樹英的支持度就是 0.14 - 0.08 = 0.06,比李永達低。若陳樹英當選,由於李永達得票比她多,因此亦必然當選。換句話說,當陳樹英的得票比李永達低的時候,她的勝算實際上就是兩者皆勝的機會率,而李永達的落選機會率就是兩人皆輸的機會率。同一道理,可知圖中紅色部份的高度,代表兩人全滅的機會率,綠區高度代表兩人全勝的機會率,而黃區高度代表兩人中只得一人當選的機會率。以李永達的得票率為 0.08 為例,兩人全勝的機會率為 0.43(圖中深綠色直線長度),兩人全滅的風險為 0.1(赤色直線長度),而李勝陳敗的機會率為 0.47(暗黃直線長度)。
圖三:從九月七日民調所得的民主黨勝算曲線
分析
政黨派兩張(實際上等於一人的)名單參選,目標不外以下兩者:
  1. (全攻型)務求全取兩席。
  2. (防守型)確保一席,並伺機提高另一張名單的勝算,情況轉壞則棄車保帥。
九月七日的民調指民主黨所享的支持度為 14%。若此支持度不變,而民主黨務求兩張名單皆勝,就要靠配票來提高全勝的機會率。下圖顯示,若要將全勝的機會保持於五成以上,李永達的得票率就必須維持於 (0.062, 0.078) 這個範圍之內,得票率的誤差不能大於 0.016(圖中綠色橫線長度),亦即所有支持票的 (0.016/0.14) × 100% = 11.4%;以去屆該區的總投票人數 398292 人來估計,即誤差不能大於 6373 票。靠配票令兩張名單得票相差不多於六千幾票,換了是民建聯應該能輕易辦到。民主黨的支持者並非鐵票部隊,但是上述這個準確度照計尚算寬鬆。

吊詭的是,要提高全勝機會,亦自然會提高全滅的風險,原因是綠區的頂峯即是紅區的尖端。最有全勝把握的時候,也是全滅風險最大之時。圖四之中,若民主黨要將全勝機會保持在五成以上,全滅的風險就會介乎 0.12 至 0.27 之間。0.27 的全滅概率,或 1 - 0.27 = 0.73 的全勝機會算不算高,我不知政黨怎麼想,我自己就覺得「唔湯唔水」,有點微妙,不過也不是去到要令人非常警惕的地步。故此,若民主黨要採取全攻型策略而不棄保,也是合理的。

只是,這背後隱藏了一個假設,就是民主黨的支持率真的如民調所說,有 14%。
圖四:總得票率 0.14,全勝機會 0.5 或以上
民意調查當然並非毫無誤差。九月七日的調查,樣本數為 503,而民主黨的支持度為 0.14。若以常用的 95% 信賴區間來計算,民主黨支持度的誤差可達 $1.96\sqrt{\frac{0.14(1-0.14)}{503}} = 0.03$ 之多。換句話說,即使民主黨實際的支持度為 0.14 - 0.03 = 0.11,也不是甚麼奇怪的事。


事實上,前述 NOW 新聞台的多次民意調查,都顯示民主黨的支持度低迷,只得 11% 左右,直到最後兩次民調,才忽然上升至 14%。如此情況下,究竟我們應該相信之前的結果,認為後來兩次有統計誤差,抑或相信民意真的變強,實在見仁見智。然而,問題是,即使民主黨的支持度真的變強了,民調的誤差也沒有 3% 那麼大,我們仍不能排除 14% 這個數字有誤差,而民主黨的真正問題,是即使它的支持度只下降少許,勝算也會大幅改變

前面提過,若李、陳兩張名單的總得票率為 $v_\max$,李永達的得票率為 $v$,則陳樹英的得票為 $v_\max-v$。因此,若李永達的得票率不變,但民主黨的總得票率 $v_\max$ 下降一個單位,陳樹英的勝算曲線就會向左移動一個單位。下圖顯示了民主黨的總得票率為 0.14, 0.13, 0.12 及 0.11 時的情況。圖中可見,若李永達的得票率為 0.08,而民主黨的總得票率由 0.14 跌至 0.13,則全勝機會會由 0.4 以上大降至 0.2 以下;若總得票率降至 0.11,而李永達的得票率並無改變,則民主黨的全勝機會更是近乎零。
圖五:民主黨的全勝機會,對總得票率十分敏感
就算民主黨的真正得票率只下降 0.01,配票也會變得異常困難。見圖六。當 $v_\max = 0.13$,全勝機會最多也不到六成。若要將全勝機會維持於五成或以上,配票的準確度就要在 0.006 之內,或支持票的 (0.006/0.13) × 100% = 4.6%;以去屆投票人數來計算,即是 2400 票左右。換了是民建聯,恐怕也未必做得到。何況選舉當日的票站調查結果以及配票都會有誤差,若對配票準確度的要求太高,配票策略很容易泡湯。困難之餘,全滅的風險亦變得相當大,介乎 0.33 與 0.41 之間。為追求不到六成的全勝機會而冒上超過三分之一的全滅風險,實在不是明智的做法。
圖六:民主黨的勝算曲線 ($v_\max=0.13$)
要減低全滅的風險又如何?見圖七。若要將風險局限於兩成以下,全勝的機會就頂多得 0.34,並不值得期待。

由此可見,即使民主黨的支持度只下跌 0.01,守勢已比攻勢來得現實和明智。
圖七:民主黨的勝算曲線 ($v_\max=0.13$)
若支持度下降到 95% 信賴區間的下限 0.11 的話,情況就更嚴重了。見圖八。兩條曲線圍成一個紅色的死亡大三角,民主黨全勝的機會頂多只有三成,全滅風險卻可以高達七成。若如上例般,想將全滅風險保持於兩成以下,則全勝機會最多只剩 4%。若這樣的情況下仍希冀全勝,簡直異想天開,棄車保帥才是唯一合理策略。
圖八:民主黨的勝算曲線 ($v_\max=0.11$)
討論
以上各種情形顯示,若民主黨的支持度真的達到如選前最後兩次民意調所指的 14%,它不棄車保帥,還算合理。今次選舉,結果民主黨的得票率只有 11.77%。如果我們據此指摘民主黨之前過份樂觀的話,只是事後孔明。然而,按九月七日的民調結果,只要民主黨的總得票率下降 0.01,全攻型策略已經變成不夠現實。鑑於民主黨的全勝機會對總體得票率太過敏感,而民意調查結果及配票也會出現誤差,因此,除非主事人懷着 wishful thinking,或者於九月七日民調結果公布後,認為自己有可靠方法提高總得票率,否則,為免全滅,民主黨當時是應該棄車保帥的。

(至於誰是車,誰是帥,就很難說。我自己就比較希望棄李保陳,但這是另話。)

既然選舉結果已經揭盅,我們難免要事後孔明一番。民主黨於新界西最後得票率為 0.1177,李永達為 0.0658。憑過往經驗,這個情況的勝算有幾高?見下圖。
圖九:事後孔明圖 ($v_\max=0.1177, v=0.0658$)
Fate/Legco
有人指摘余若薇參選分薄了民主黨票源,連累民主黨全滅。余若薇分薄了民主黨票源固然是事實,然而,根據上述分析結果,按當時形勢,民主黨是應該棄車保帥的。它最後得票 11.77%,九月七日的民調結果則為 14%,兩者都夠一張名單越過安全線。要避免全滅,保持原有的一席,技術上十分容易,但它沒有這樣做。因此,今次全滅完全是民主黨自己的責任,與人無尤。

然而,亦有人說,若余若薇不參選(或者不落力競選第二席,而只是保郭家麒當選),則民主黨可全取兩席。故此,就算李永達落選的責任不在余若薇,她仍是連累了陳樹英丟了另一席。

我不能同意這種論調。誠然,若余若薇不參選,就算支持余若薇的票並不完全流向民主黨,添給陳樹英的部份應該也足夠令她當選。事後看來,要將余若薇的票過給李永達與陳樹英,令民主黨取得兩席,亦比將陳樹英的票過給李永達與余若薇,令他們連任容易。然而,說余若薇連累陳樹英,講到這個議席彷彿本來就屬於陳樹英一樣,是十分奇怪的。若問余若薇為何不退選,為何不反過來問,何以陳樹英不退選?候選人參選,是為了提供一個選擇給選民;選民投票,是為了選出能夠代表他們的議員。議席並非議員的私有物,民主派的候選人,對選民來說,更非隨便一個也可以代表他們,個個一樣的替代品。即使議席是議員的私有物,余若薇本來就有一個議席,她有何義務要將議席拱手相讓給本來並無議席的陳樹英?

有兩夫婦與一位鄰居遇到海難,三人各乘一艘只載一人的充氣快艇逃生。大海茫茫,暴風雨迫近。他們估計所餘燃料只夠二人逃走,究竟誰要犧牲?最後三人相持不下,全數沒頂,這是誰的責任?

有人指摘該名鄰居不肯自我犧牲,成全兩夫婦。你話呢?

參考網頁
伸延閱讀

2012年9月15日星期六

Normal approximation of multinomial distribution

How to simulate relative frequency outcomes of a multinomial experiment using normally distributed random numbers? The answer is surprisingly simple, but for some curious reason, it is seldom mentioned on the internet.

Let $\mathbf{X} = (X_1, \ldots, X_m) \stackrel{\textrm{i.i.d.}}{\sim} \textrm{Multinomial}(n, \mathbf{p})$, where $\mathbf{p}=(p_1,\ldots,p_m)$ is a probability vector whose entries sum to 1. In other words, we are talking about $n$ independent trials of a multinomial experiement, in which the probability of getting outcome $i\in\{1, 2, \ldots, m\}$ in each trial is $p_i$, and $X_i$ is the frequency count for outcome $i$ after all $n$ trials are completed. We have the following:

Theorem. Let $u = (\sqrt{p_1},\ldots,\sqrt{p_m})$ and $Q$ be a real orthogonal matrix whose last column is $u$. Suppose $Z_1, \ldots, Z_{m-1}$ are $m-1$ i.i.d. standard normal random variables. Then
$$
\mathbf{X}^\ast = \left(
\frac{X_i-np_i}{\sqrt{np_i}}
\right)_{1\le i\le m}
\stackrel{d}{\longrightarrow}\quad Q \begin{bmatrix} Z_1\\ \vdots\\ Z_{m-1}\\ 0\end{bmatrix}
$$as $n\rightarrow\infty$.

Proof. See sec. 11.1 of Hans-Otto Georgii (2007), Stochastics: Introduction to Probability and Statistics, Walter de Gruyter, Berlin.

*****************************
Note that the denominator of the $i$-th entry of $\mathbf{X}^\ast$ in the above theorem is $\sqrt{np_i}$, not the usual $\sqrt{np_i(1-p_i)}$ we see in the normal approximation formula to binomial distribution. We will immediately see why the binomial case is just a special case of the above general formula. Let $m=2$ and write $\mathbf{p}=(p,q)^\top$. Take $Q$ as $\begin{pmatrix}\sqrt{q} & \sqrt{p}\\ -\sqrt{p} & \sqrt{q}\end{pmatrix}$. So, the above theorem says that
$$
\begin{bmatrix} \frac{X_1-np}{\sqrt{np}}\\ \frac{X_2-nq}{\sqrt{nq}}\end{bmatrix}\stackrel{d}{\longrightarrow}
\begin{pmatrix}\sqrt{q} & \sqrt{p}\\ -\sqrt{p} & \sqrt{q}\end{pmatrix}
\begin{bmatrix} z\\ 0\end{bmatrix}= \begin{bmatrix} \sqrt{q}\,Z\\ -\sqrt{p}\,Z\end{bmatrix}.
$$Hence, by Continuous Mapping Theorem, we can divide both sides entrywise by $(\sqrt{q}, -\sqrt{p})^\top$ and get
$$
\begin{bmatrix} \frac{X_1-np}{\sqrt{npq}}\\ \frac{nq-X_2}{\sqrt{npq}}\end{bmatrix}

\stackrel{d}{\longrightarrow}
\begin{bmatrix} Z\\ Z\end{bmatrix}.$$
Since $X_1+X_2= n$, we have $nq-X_2=X_1-np$. Hence the above convergence reduces to the familiar normal approximation formula to binomial distribution.

The theorem requires the use of a real orthogonal matrix $Q$ whose last column is $u$. How to construct such a matrix? See my previous blog entry.

By the theorem, the covariance matrix of the approximating distribution is given by $Q\ \mathrm{diag}(1,\ldots,1,0)\ Q^\top = I-uu^\top$. This matrix is, of course, degenerate because the $X_i$'s are not independent of each other (as $\sum_{i=1}^m X_i=n$).

2012年9月9日星期日

A handy completion of orthogonal matrix from a column vector

Given a unit vector $u=(u_1,u_2,\ldots,u_n)^\top$, we seek to construct an orthogonal matrix $Q$ whose first column is $u$.

If $u=e_1=(1,0,\ldots,0)^\top$, clearly we can pick $Q=I$. Suppose $u\ne e_1$. Then the Householder reflection $Q = I - 2vv^\top$ will do, where
$$
v = \frac{u-e_1}{\|u-e_1\|}.
$$
In particular, if $u=\frac{1}{\sqrt{n}}e=(\frac{1}{\sqrt{n}},\frac{1}{\sqrt{n}},\ldots,\frac{1}{\sqrt{n}})^\top$, we may choose
$$
Q=\pmatrix{a&a&\cdots&\cdots&a\\
a&b+1&b&\cdots&b\\
\vdots&b&b+1&\ddots&\vdots\\
\vdots&\vdots&\ddots&\ddots&b\\
a&b&\cdots&b&b+1}
$$
where $a=\frac{1}{\sqrt{n}}$ and $b=\frac{-1}{n-\sqrt{n}}$.

One question remains. When $u\approx e_1$, how stable is this construction?

2012年4月14日星期六

淺談 Google 的 PageRank 算則

之前報讀了 Udacity 的 CS101: Building a Search Engine,開課後才發現它比想象中淺得多,大多數時候,與其說該課程教人建構搜尋器,不如說它是 Introduction to the Python Programming Language,不過課程終段的確牽涉到 Google 的搜尋結果排名算法,只是也許為了遷就學生的程度,避談了本來可以用數學來解釋得清清楚楚的部份。這裏我就不妨補充一下。

Google 的搜尋結果排名法,稱為 PageRank。經過多年發展,現時 PageRank 算法的內容已經不為人知,不過 Larry Page 與 Sergey Brin (Google 兩位創辦人)起初是有將描述此算則的論文擺上網的,有興趣的讀者不妨看一看。

簡單來說,早期的 PageRank 是一個投票機制,每引用一個網頁,就等於投了該網頁一票。舉例說,假設下圖十一個方格都是提及動畫 Steins;Gate 的網頁,其中每個箭咀代表一項連結,亦即是 a 網頁有一道連向 A 網頁的連結,而 c 網頁內有三道連向三個不同網頁的連結等等。由於 b, c 各得三票,a 網頁卻無票(無人引用),驟眼看來,b, c 兩網頁理應比 a 網頁更權威。


然而,若純粹以票數多寡來作排名指標,會引起明顯問題。譬如:
  • 儘管圖中 A, B 兩個網頁都各得一票,但 B 網頁應該比 A 網頁重要,原因是投票給 B 的 b 網頁有人引用,投票給 A 的 a 網頁卻無。投票者較權威,得票者亦應該更權威。
  • 儘管 B, C 均各得一票,但 B 的排名應該比 C 高,原因是投票給它們的 b, c 兩個網頁同樣權威,但是 b 網頁將票全投給 B,而 c 網頁卻投了三個網頁,所以 C 網頁未必如 B 網頁那麼有用。
以上兩點,第一點其實是說我們不應只考慮得票多寡,還要考慮投票者本身的重要性(即排名),而第二點就是說,除了要考慮得票多寡(即有幾多條 incoming links)之外,亦要考慮投票者的投票若干(即投票者有幾多條 outgoing links)。要解決這兩點,可以用以下做法。假設整個整個互聯網中,有 $n$ 個網頁提及 Steins;Gate 動畫,而我們以符號 $p_i$ 來代表等 $j$ 個網頁的排名分數,分數愈高則名次愈高。那麼,我們可將 $p_j$ 定義為下列方程式的解:
$$p_i = \sum_{\textrm{page } j \textrm{ has link to page } i} \frac{p_j}{\textrm{no. of out-links of page } j}.$$
根據以上方程,若 $j$ 網頁投了票給 $i$ 網頁,那麼 $j$ 網頁本身的分數 $p_j$,亦反映在 $p_i$ 的計算之中,因此顧及了上面提到的第一點。此外,等式右邊的分母是投票者所投的票數。若 $j$ 網頁引用了許多網頁,而 $i$ 網頁只是其中之一,那麼這個大分母會將分數拉低,因而解決了上面提到的第二點。

Google 最初期的 PageRank,其實還有另一步 ── 它實際上將 $p_i$ 定義為
$$p_i = \frac{1-d}{n} + d\left(\sum_{\textrm{page } j \textrm{ has link to page } i} \frac{p_j}{\textrm{no. of out-links of page } j}\right).$$
換句話來說,除了前述的排名方法之外,PageRank 還考慮了另一個排名方法,就是各網頁皆有同等名次。將一票分給 $n$ 個網頁,每個網頁就有 $\frac{1}{n}$ 票。將兩種票數以權數 $d$ 來作加權平均,才是真正得票。Page and Brin 將這個權數 $d$ 稱為 damping factor。他們起初所用的數值為 $d=0.85$。

至於何以要作加權平均,他們含糊其詞,說 PageRank 可以視為搜尋用戶的行為的數學模型,而上式 $\frac{1}{n}$ 的部份,可視作當一名用戶厭倦了點擊網頁的連結以後,隨機探訪網上一個任意網頁。他們於另一篇文章又提到這樣可以給每個網頁一些內在價值,有助避免某單一網頁有過大影響云云,不過我感覺他們當時其實並無一套清楚的理論,解釋加權的原因,只是單純覺得這樣做有助改善搜尋器的表現而已。讀者切勿誤會,以為我看扁他們,實情剛剛相反。儘管我以為他們並無一套理論解釋何以要加權,但就正正因為如此,才顯出他們是上等的工程師,有靈敏的直覺。

PageRank 的加權平均方法,其實是有理論基礎的。假設 $C_j$ 代表網頁 $j$ 有幾多條 out-links,而 $g_{ij}$ 代表「網頁 $j$ 引用了網頁 $i$」的布林值,亦即當網頁 $j$ 引用了網頁 $i$ 的時候,$g_{ij}$ 為 $1$,否則為 $0$。於是我們可以將 $p_i$ 的定義以矩陣形式寫成
$$\begin{eqnarray}
\mathbf{p} &=& \left[\frac{1-d}{n}
\begin{pmatrix}1&\ldots&1\\\vdots&&\vdots\\1&\ldots&1\end{pmatrix}
+ d\mathbf{G}\begin{pmatrix}\frac{1}{C_1}\\&\ddots&\\&&\frac{1}{C_n}\end{pmatrix}\right] \mathbf{p}\\
&=& \mathbf{A}\mathbf{p} \quad\textrm{ (say)}.
\end{eqnarray}$$
我們要問的問題有三個。第一,$\mathbf{p}=\mathbf{A}\mathbf{p}$ 此方程是否有解?第二,若它有解,又是否只有獨一的解?若有多個解,豈非說一個網頁有多個排名?當然,若 $\mathbf{p}$ 是方程式的解,將 $\mathbf{p}$ 乘上任意一個實數 $\lambda$,亦將是一個解,因此我們必須將方程式的解正規化 (normalized)。由於 $p_i$ 可以理解為「$i$ 網頁乃有用的搜尋結果」的機會率,所以我們要求 $\sum_i p_i=1$。第三,亦由於我們將 $p_i$ 視為機會率,故此亦要求 $p_i\ge0$。所以,綜合起來,我們要問的,就是:
方程式 $\mathbf{p}=\mathbf{A}\mathbf{p}\quad (p_i\ge0,\ \sum_ip_i=1)$ 是否有獨一的解?
不難驗證,$\mathbf{A}$ 是一個非負 (entrywise nonnegative) 矩陣,且每列元素總和為 1 (column sum = 1)。數學上,這樣的矩陣稱為隨機矩陣 (stochastic matrix)。矩陣理論中,有一條與隨機矩陣相關的定理,稱為 Perron-Frobenius Theorem,內容是說,若 $\mathbf{A}$ 是一個元素為正的隨機矩陣,那麼它會剛好有一個特徵值 (eigenvalue) 為 1,而且所有其他特徵值的絕對值皆小於 1;此外,撇除放大與縮小 $\mathbf{p}$,方程 $\mathbf{p}=\mathbf{A}\mathbf{p}$ 有獨一的解,兼且此解的所有元素為正數。

有了這條定理,PageRank 算則之中權數 $d$ 的功用就呼之欲出了 ── 若不取加權平均(或等價地將權數設為 $1$),$\mathbf{A}$ 只是個非負矩陣,而非正矩陣,因此不能保證 $\mathbf{p}=\mathbf{A}\mathbf{p}$ 有獨一的解,亦即可能存在多種排名次序,結果要用其他方法解決排名問題。然而,一旦有 $0\le d<1$,我們即可確保 Google 大神分得出排名先後。

回說 Udacity 的 CS101。主講的 Prof. Evans 避談數學,將 PageRank 以迭代方式定義為 $\mathbf{p}_t =\mathbf{A}\mathbf{p}_{t-1}$。這其實於數學上稱為 power method,是一種計算特徵向量的方法。Brin and Page 的文章說他們用了某種 simple iterative algorithm 來計算 PageRank,所指似乎就是 power method。此方法的名稱,源自 $\mathbf{p}_t =\mathbf{A}\mathbf{p}_{t-1} = \mathbf{A}^2\mathbf{p}_{t-2} = \ldots = \mathbf{A}^t \mathbf{p}_0$,亦即利用 $\mathbf{A}$ 的 power 去計算特徵向量 $\mathbf{p}$。這個方法的好處是容易理解、實作簡單,缺點是當 $\mathbf{A}$ 最大和第二大的特徵值相當接近的時候,power method 會收斂得非常慢。一般而言,power method 並非尋找特徵向量的好方法,不過做網頁排名的話,$\mathbf{A}$ 是一個超巨大的稀疏矩陣 (sparse matrix),兼且由於 web crawler 會不時更新矩陣內容,因此簡單的方法可能更有利。實際情況如何,就在我知識範圍以外了。

上文是我將多年前的舊文略為變更而成,目的是以最少的線性代數解釋早期 PageRank 的數學原理。美國數學學會 (American Mathematical Society) 另有一篇 feature column 做類似的事,但遠較本文詳盡,有興趣的讀者不妨一看。本文並未參考 AMS 的文章。

2012年3月18日星期日

Caltech's online ML course

繼 Stanford U. 與 UC Berkeley,今次輪到 Caltech 開辦網上課程,主題是 Learning From Data,是一個 Machine Learning course,四月開課,從三月廿六日起接受報名。課程網頁強調此乃真正的加州理工課程,會直播課堂,亦絕無將內容淺化:
Real Caltech course, not watered-down version
Broadcast live from the lecture hall at Caltech
整個課程為時兩月,共十八課,看來非常認真。據稱功課有深有淺,有理論性的題目,也有編程習作,學生成績會登在公開的分數版上,太介意分數者要三思。

2012年3月6日星期二

Ruby 初體驗

史丹福大學的網上課程出現問題,未能如期開辦,於是我唔嫁又嫁,上了兩星期 SaaS 的課。課堂乃 UC Berkeley 本科生課堂內容的錄影,完全是正宗名校課堂內容,只不知功課程度與本科生的差幾遠。沒有同學、助教或教授支援,實在有點吃力,結果,第一份功課已經覺得太難,搞到要上網抄功課。抄襲本是學生大忌,但這並非正規課程,而且根據課程簡介所述:
Those who submit homework 1 and receive a passing grade will receive a coupon good for 100 hours of small instances of EC2 for use on the remaining homework assignments plus a coupon to upgrade their free GitHub accounts to a Micro account (both good through the end of course).
有此物質誘因之下,我實在很想學曉如何在 Amazon 設置雲端應用,那管該應用是如何白癡。第三週的功課已經要求學生做 application deployment,真懷疑自己能否完成。

SaaS 課堂用的是 Ruby。Ruby 感覺上有些似 Smalltalk,是很徹底的 OOP(連 class 都是 Class 的一個 instance!)加 dynamic typing。猶如 Java 推出初期,C/C++ 與 Java 的用家之間有 language war 一樣,Ruby 與 Python 的粉絲之間,好像也有互相攻訐。有些粉絲將某些 features 吹噓到天下無敵,以求將對手比下去。我不熟悉 Ruby,對 Python 也只是知道一點點,不敢妄言,但感覺上,Ruby 某些所謂 killer features (例如 block)好像也不是那麼突出,反而是一些小處更有趣,例如 Ruby 的 method names 可以包含 ?(問號)或 !(嘆號)。前者用來表示該 method 乃是非題,會回傳布林值。譬如我想知道某變數 x 是否屬於整數類別,就可以說 x.kind_of? Integer,換了是其他語言,類似的表述就可能寫成 is_integer(x),儘管差不多同樣清晰,但是 method name 有一個問號,讀者就絕不會搞錯這個 method 的意圖。至於感嘆號,乃用來表示該 method 具破壞性,舉例如 y 是一個 array,那麼,呼叫 y.sort! 會作 in-place sort,亦即 y 本身經過排序之後,內容會改變。相比之下 y.sort 只會回傳一個經過排序的 y 的 copy,但 y 本身保持不變。比起 Python 用 y.sort()sorted(y) 來表示兩個版本,Ruby 明顯高出一籌。

初學一種新語言,難免有不習慣的地方。儘管 Ruby 與 Python 都有不少 functional programming 成份,但不知何故,我總覺得 Ruby 的 funcional programming 味道較濃。我唸大學時讀的是數理科目,基本上,除了覺得 lambda functions 相當有用之外(特別是做 numerical optimization 的時候),很少接觸其他 functional progamming 技巧。學了 Python 之後,也只是覺得 list comprehension 很厲害,但其他 functional programming 部份有乜用,我是很疑惑的。事實上,就連 Python 之父 Guido van Rossum 也說一些 functional programming 的程式很難明白,因此他設計 Python 3.x 時,一度想將 Python 2.x 之中的 lambda, map(), reduce(), filter() 除掉三個,只保留 reduce(),但他又說
So now reduce(). This is actually the one I've always hated most, because, apart from a few examples involving + or *, almost every time I see a reduce() call with a non-trivial function argument, I need to grab pen and paper to diagram what's actually being fed into that function before I understand what the reduce() is supposed to do.
當然,如著名開發者 Joel Spolsky 所講,MapReduce 可以是超有用的,不過若電腦的 interpreter 或 compiler 並不利用 map 或 reduce 作平行處理,又或者我們用到的迴圈有 side effect,那麼我覺得 map 或者 reduce 並不很有用場。利用它們,你也許可以將程式由幾行變成一行,但代價是對一般人來說,程式反而變得較難閱讀,而這也是我對 Ruby 感到最不習慣的地方 ── Ruby programmers 似乎有一種意識形態,就是愈少行數,即代表程式愈 elegant。老實說,我覺得這種想法不但搞混了 compactness 等同 conciseness 兩種概念,還太過精英主義,有些變態。

2012年2月17日星期五

Udacity 的網上課程

之前史丹福大學的網上課程大受歡迎,大概需求大增,準備需時,來季課程要比原訂時間押後數星期才開始。我原本報讀了 SaaS, Anatomy 及 NLP,但 SaaS 竟然要學生買課本!身為 cheap 精,決定打退堂鼓。況且 SaaS 課程會用 Ruby on Rails 的 framework,而我本身的 programming skills 已經夠屎,暫時實在無 mood 學多一種 programming language。

去季教授 AI 的 Prof. Sebastian Thrun 另起爐灶,開了一個 Udacity 網站。本月二十日起,將開設CS 101: Building a Search Engine 及 CS 373: Programming a Robotic Car 兩個課程,相當吸引。Prof. Thrun 領導的史丹福大學隊是第二屆 DARPA Grand Challenge(由美國國防部資助)的冠軍,也是兩屆比賽中首個能夠讓越野車自動走完整個賽程的隊伍(首屆直情無車走完全程),他毫無疑問是有關領域首屈一指的專家。當然,現在這個課程應該不會碰到硬件層面,大概仍是教授 particle filter 等等的基本 AI/ML 原理,但已經足夠吸引我報讀。課程簡介說會以 Python 作程式語言,這又省卻了學習新語言的功夫。

至於 CS101: Building a Search Engine,由於是入門課程,可想而知,大抵又是教授初期的 Google PageRank 那套。稍有數學根底的人都知道,初期的 Google PageRank 只不過是矩陣理論入面 Perron-Frobenius Theorem 的一個簡單應用,其原理,十五分鐘已經可以講完。現在化為一個長達六週的課程,我希望是牽涉一些其他方面的知識,例如怎樣寫 web crawler,如果 tokenize 一個網頁等等。這些非數學類的編程知識,是我最希望學到的。

報讀了這兩個新課程後,就要考慮是否放棄其他已經報讀的課程。Anatomy 當睇 Discovery Channel,課業應該很輕,照去。NLP 感覺最困難,未知能否兼顧,但這個 area 依然方興未艾,值得學習。看來仍是見步行步,睇送食飯。

2011年12月21日星期三

The joy of learning

呢個標題,是抄自電鋸的。史丹福今季的 AI 與 ML 課程剛剛完結,距離下一季開課還有一個月,於是我昨晚急急腳上網報讀同樣剛剛完結DB 課程。雖然任何人事後都可以上網睇返 lectures 的影片,但係無報讀的話,就無得觀看功課內容,更遑論做功課或參加考試。

據說 DB 乃三個課程之中,難度最接近正式的本科生課程者。我一口氣完成了第一週的課程,覺得 DB 在三個課程當中,真係「有 D 嘢」,不但老師教得最好,功課亦設計得最佳,令我每條題目都學到嘢。我無諗過自己學 DB 呢類睇落好 Q 悶的課題,結果反而最留心聽。例如以前好多次試過想學 XML,但係拿起書讀幾頁就已經見周公,但係睇史丹福的 lectures,做了幾題寫 DTD 的 hands-on exercises,就即刻入腦。第四份功課講 Relational Algebra,最後一題好難,但係又唔會令人氣餒。最後解決到的時候,真係好有滿足感,玩得好過癮。下一課係 SQL,期待更有趣的挑戰。

2011年11月17日星期四

史丹福 2012 年春季公開課程

Stanford U. 的工程學院於來季會開設多項免費、公開但無學分的網上課程,而該校其他學院與另外一些大學院校亦紛紛響應,例如 UC Berkeley 就開設下述的 SaaS 課程。比起 Stanford 今季只開三科,來季起碼合共有七個九個超過十個不同科目,執筆時數目仍在增加,不能盡錄。有興趣的朋友不妨報名。

身為過來人,奉勸各位,無論課程有多吸引,最好還是只讀一科,否則時間實在太吃緊(我今季報讀了兩科,而且課程中近半都是我原本就懂得的東西,但仍是感到「困身」和吃力),但各位可以先報多個課程,然後才決定主攻那一科,而放棄其餘。
後記:我報讀了 SaaS 及 NLP,開學後看情形才決定是否放棄其中一科。NLP(是 natural language processing,不是那種騙錢的 NLP)是我一直想知的範疇。SAAS 則不然 ── 事實上我覺得自己「未夠班」去學 ── 但它的課程看來十分實用,又似有 hands on 的 programming exercises,所以一試。教授們的英語也是我考慮報讀的原因。看兩個課程的簡介短片,四位教授的英語都非常清晰,100% 聽得懂,比起今季 AI class 兩位教授,實在好太多了。

2011年10月19日星期三

SEE

前文提到史丹福大學設立了三個免費、公開但無學分的網上課程,其實它還有其他課程,雖不公開報讀,但課程的材料是公開的,例如有一科 iPhone Application Development,有志開發 iPhone Apps 的朋友不妨一看。校方稱它的公開課程計劃為 Stanford Engineering Everywhere (SEE),目的是為了與 MIT OpenCourseware 競爭。

AI and ML

無心寫文章,吹吹其他水。

月前於《小城科學》blog 及電鋸處看到史丹福大學開了 database, artificial intelligencemachine learning 三科免費、公開但無學分的網上課程,儘管我全部都有興趣,但是 database 方面,總覺得若無實際問題在手,齋聽書還是讀不通;況且現今資料庫的應用,十居八九都與網絡有關,可是自己對網絡一竅不通,所以還是作罷。其餘兩科,本來只報其一,時間上比較鬆動,但心癢之下,還是兩科都報讀了。

根據校方數字,最後每科都有幾萬人報讀。面向如此大的群體,又除了聲譽之外沒有實利,課程自然較本科生所念的淺,許多材料會被 heavily dumbed down。例如看 AI 課的學生論壇,AI 的本科課程 CS221 於頭兩週過後,就要學生用課堂中所教的 A* search 做一個 project,寫一個類似用於 Pac-Man 遊戲的算則。參加 AI 網上網程的學生,不但毋須做 project,功課也大多只是網上選擇題,難度低很多。(後記:出乎意料,原來 AI 課的功課與考試與給予正規學生的相同,但 ML 的網上功課就和正規生的有別。)

AI 第一週所教的,大多與其他大學課程有重疊,例如用 BFS, DFS 搜索樹形圖等等,這些也是運籌學 (Operational Research) 的標準內容,不過學習一下 CS 佬的觀點也不是壞事。然而,不知是教授本身還是 CS 佬的慣例,課堂中好些語彙的用法,似乎都大大偏離學界常規。例如科學上,當我們說 discrete problem 與 continuous problem 的時候,discrete 的可以是 finite,也可以是 countably infinite,總之就是 countable。只有有 uncountably infinitely many states 的問題,才稱為 continuous problem。可是根據教授的說法,discrete problem 就是有 finitely many states 的問題,其餘的一律稱為 continuous problems。路徑的長度是另一個例子。教授稱 BFS 為 shortest first search,他又說 BFS 與 uniform cost search 都能夠找出 optimal path。然而 uniform cost search 尋求的,是最低代價 (cost) 的路徑,而 BFS 所得到的,只是一條最少節點的路徑,而完全不理會路徑的代價。要說 BFS 保證找到 optimal path 也可以,只不過這個 "optimal",指的是節點或層級的數目,而非路徑的代價。在其他學科中,例如圖論或運籌學,arc cost 歸 arc cost,no. of arcs 歸 no. of arcs,兩者決不輕易混為一談。

教授於 lectures, quizzes 跟 homework 的遣詞用字,亦常常過於含糊。這並非我獨有的印象,也是學生論壇裏的主流意見。甚至有些 quizzes,連教授到底想問甚麼,我也搞不清楚。感覺上,AI 兩位教授不是很 well prepared,有點急就章,不過他們的 lectures 很有啟發性,例如有一處談到 A* search,教授問,A* search 行得通,當中的 intelligence 究竟從何而來?要留意,並非人人也將 search method 當是 AI 的,例如早年深藍擊敗國際棋王卡斯巴洛夫,後來負責設計算則的許峰雄來港,就提及他的算則不過是 brute-force search,算不上是 AI。然而教授的問題,為何 A* search 行得通,就真的令我不禁要停下來,想一想,而他的答案,也令我有恍然大悟之感。

相比之下,ML 的教授較注重包裝,無論是 presentation 抑或 website 的設計都很講究。ML 的教授 Andrew Ng 風格有點「執手教」,即是連很顯淺的東西也唯恐你不明白,所以解釋得很仔細。要解釋詳盡但不冗贅,並非人人都做得到,華人尤其傾向太注重技術細節,令人失去 big picture,沒有 motivation,而 Andrew Ng 是罕見的例外。不過和 AI 課相比,我還是覺得 AI 課較能刺激思考,不知是否西人與華人的治學方式始終有別。只是現在始終開課不久,日後也許會有所不同。

儘管 Andrew Ng 的 presentation 很好,但也有些我不喜歡的地方,尤其是他有時舉一些很不設實際的例,很容易「教壞人」。譬如他解說 linear regression,以樓價與樓面面積的關係作例子。樓價的研究確實有用得上 linear regression 的地方,但是實際做法是有成例的,例如樓價幾乎一定要 take logarithm,而且,由於樓面面積幾乎一定不是決定樓價的唯一主要因素,若不考慮其他因素(例如地區、座向、層數、交通等等),regression 得出來的結果差不多肯定沒用,但是引入其他因素的話,又幾乎必定牽涉 hedonic regression 的概念。實例可以淺化,但不能偏離正軌,「老作」一個例子,然後硬套入現實場景,很容易誤導學生,令他們以為隨隨便便放幾個變數,就可以做 linear regression。其實隨便抽一本計量經濟學 (Econometrics) 的書,也可以找到許多實例,若無實例在手,還是只講抽象例子為妙。

此外,ML 課的內容也不時有錯。ML 要用到其他學科的技巧,而教授不是那些科目的專家,所以犯錯也情有可原,但是向學生胡亂解釋,就會造成真正問題。例如課堂中有處(大意)指,若要 minimize $\|X\theta-y\|^2$,其中 $X$ 是 $m\times n$,而 $n$ 遠大於 $m$(亦即 $X$ 是闊闊的矩陣),就應該用 gradient descent 而非 normal equation 來解決,原因是在 normal equation 之中,要計算 $(X^\top X)^{-1}X^\top$ 的話,由於 $n$ 很大,會很花時間云云。問題是,由於 $n>m$,$X^\top X$ 並不滿秩 (rank-deficient),所以它根本就不能反逆!最奇怪的,是 Andrew Ng 於 lecture 中有提及 pseudoinverse,可見他應該知道,當 $n>m$ 的時候,我們要計算的,應該是 $X^+$ 而不是 $(X^\top X)^{-1}X^\top$。為何仍有上述錯誤,真是木宰羊。(又後記:先前我跳過了許多我本身懂得的內容,現在打開來看,才發現教授並不真正熟悉 multiple regression ── 在 "normal equation" 的 video 約 10:23,他說用 Octave 解 normal equation 的時候,應該用 $\theta=\textrm{pinv}(X'\ast X)\ast X'\ast y$。技術上這沒有錯,但一般通用、等價而且較簡單的答案,其實是 $\theta=\textrm{pinv}(X)\ast y$,實際計算上,我們更罕會先算 $\textrm{pinv}(X)$,再算 $\theta=\textrm{pinv}(X)\ast y$,而是用諸如 QR factorisation 等等的數值方法尋求方程式 $X\theta=y$ 的解。Andrew Ng 取 $(X^\top X)^+X^\top$ 而捨 $X^+$,令人愕然。)

課堂中有關 feature normalization 的討論,更是完全錯誤。教授說 feature normalization 的目的,是為了令 gradient descent 加快收斂,但這兩件事,其實風馬牛不相及。試想像,若 features variables 未 normalized 之前,objective function 的 contour plot 本身已是同心圓狀,那麼,經過 feature normalization,contour plot 變成橢圓形,gradient descent method 豈非收斂得更慢,而不是更快?真正要改良 gradient descent,化橢圓為正圓,應該用 conjugate gradient method。Feature normalization 其實只是單純地從按每個 feature 的數值範圍 ── 而非 objective function landscape ── 去改變該 feature 的 learning rate,與加快/減慢 gradient descent method 的收斂,關係不大。

ML 的功課也設計得很奇怪。第一週的功課有兩部份,首部份是 ordinary linear regression,必答;第二部份是 multiple regression,是 bonus part。教授大概想弄得愈淺愈好,結果所謂功課,不過是在每個教授預先寫好的 script file 中加入一行指令。然而 multiple regression 部份要求加入的程式指令,其實與 ordinary linear regression 部份的完全相同,因此只要 OLS 部份答對,就等於懂得 bonus part,根本沒有額外挑戰。

另外,教授聲稱可以用 Matlab 或 Octave 來做功課,但是 submit 功課的 script file 其實呼叫了 Octave 的 urlread() function,所以 Octave 其實是不裝不行。然而我電腦 (Windows XP) 上的 Octave 又好像很 buggy,只要呼喚任何 plotting functions 就會 crash(後記,問題已解決;詳見此),令我不得不先用 Matlab 做好功課,再用 Octave 提交,但校方提供的 Octave script file,又時不時與 Matlab 不相容。例如 Matlab 的 function 應該用 "return" 來結束,但 ML 的 Octave script file 就用 "end";向量的長度,在 Matlab 是 length(),但是 ML 課的 script 就用 Octave 的 numel()。結果我要先修改那些聲稱與 Matlab 相容的 Octave script files,才可以做功課,十分麻煩。