Lisp-Stat 翻譯 —— 第二章 Lisp-Stat教程

第二章 一份Lisp-Stat教程

    本章打算將使用Lisp-Stat做爲統計計算器和繪圖器作一個介紹。在介紹完基本的數值和繪圖操做以後,介紹如何構建隨機的和系統的數據集合,如何修改數據集,如何使用一些內建的動態繪圖工具,如何構建線性迴歸模型。最後兩節給出了簡略的介紹,關於如何寫你本身的函數和使用功能數據進行高級建模的工具。javascript

2.1 Lisp解釋器

    你與Lisp-Stat系統的交互是由你和lisp解釋器之間的對話組成的。想象一下,你坐在計算機前打開了你的系統,或者更好點兒,你獲取了某個版本的Lisp-Stat,而且啓動了它。當準備好要開始對話時,Lisp-Stat在一個空白行的前邊給出了一個提示符,就像這樣:java

>

若是你鍵入一個表達式,解釋器就會打印這個表達式的計算結果來響應你。例如,若是你給解釋器一個數字,而後按回車鍵,那麼解釋器就會響應你——簡單地 打印 這個數字,而後 在下一行 退出並給你一個新的提示符。算法

 

> 1
1
>

對數字的操做能夠經過組合數字和一個符號並把它們複合成一個表達式來執行,就像這樣(+ 1 2):shell

 

> (+ 1 2)
3
>

就像你猜測的同樣,這個表達式就是將數字1和2加在一塊兒。這種操做符放在操做數前邊的表示方法叫前置表示法。這種表示方法起初可能讓你很迷惑,由於它背離了標準的數學實踐,可是卻會帶來一些頗有意義的優點。一個優點是能夠容納進行任意數量的參數操做。例如,編程

 

> (+ 1 2 3)
6
> (* 2 3.7 8.2)
60.68
>

    爲了理解解釋器是如何工做的,要記住的一條基本的原則是全部的東西都是能夠被求值的。數字求值獲得它自己。windows

 

> 416
416
> 3.14
3.14
> -1
-1
>

還有一些其它的基本數據形式求值後也會獲得它自己。它們包括邏輯值數組

 

> t
T
> nil
NIL
>

還有字符串,它們將用雙引號括起來。數據結構

 

> "This is a string 1 2 3 4"
"This is a string 1 2 3 4"
>

分號";"是Lisp註釋符,任何在分號後鍵入的字符,直到你下次鍵入回車,其間的全部內容都會被解釋器忽略。app

    若是解釋器接收到一個符號(symbol,Lisp術語),像"pi"或者"this-is-a-symbol",來計算的話,它將首先檢查是否有一個值與該符號關聯,若是有值關聯就返回那個值:less

 

> pi
3.141593
>

若是沒有值與該符號關聯,解釋器將發出一個錯誤信號:

 

> x
error: unbound variable - X
>

符號pi是系統預約義的,做爲數字π的一個近似值。從這點來看,符號x是沒有值的。下一節咱們將看到如何給一個符號賦值。

    當你敲寫符號名稱時,能夠用大寫形式也能夠用小寫形式。Lisp內部會將小寫字母轉換爲大寫形式。

    組合表達式的求值有一點點複雜。一個組合表達式是由一些放在小括號內部對的,用空格分開的元素的列表組成的。列表的元素能夠是字符串、數字或者其它的組合表達式。在對形如(+ 1 2 3)這樣的表達式求值時,解釋器把這個列表當成一個函數表達式。列表的第一個元素表明一個Lisp函數——能夠接受一些參數而且返回一個結果的代碼段。本例中的函數式加法函數,由符號"+"表示,列表的其他元素是它的參數,這裏是一、2和3。當要求對錶達式求值時,解釋器調用加法函數,把它應用到參數上,而後返回結果:

 

> (+ 1 2 3)
6

    函數的參數在函數使用以前已經被求值。在前邊的例子裏,參數都是數字,所以求值結果爲自身。可是在下邊這個例子裏:

 

> (+ (* 2 3) 4)
10

在加法函數使用以前,解釋器不得不先對第一個參數(* 2 3)求值,再傳遞給加法函數。

    數字、字符串和符號是一些在Lisp裏能夠直接使用的基本數據類型。我將會指出這些基本數據類型共同當作簡單數據。還有,Lisp提供了一些方法來造成更復雜的數據結構,那些叫複合數據。最基本的複合數據就是列表——list。列表可使用"list"函數來構造:

 

> (list 1 2 3 4)
(1 2 3 4)

解釋器打印的這條結果很像咱們已經用過的那個複合表達式:一個放在小括號裏的,由空格分開的元素序列。事實上,Lisp表達式就是簡單的列表。

    一些其它複合數據形式也是可用的,包括矢量(vector)和多維數組。這些在4.4節再討論。

    在開始第一個統計應用程序以前,讓咱們再看一個解釋器的特性。假設咱們想要解釋器打印一個列表(+ 1 2)給咱們,直接把這個列表鍵入解釋器是不行的,由於解釋器會堅持對該列表求值:

 

> (+ 1 2)
3

解決辦法就是告訴解釋器不要對列表求值。這個處理過程叫作引用(quoting)。就像這樣,在一個手寫的句子裏,將一個表達式周圍加上引號,而後說:把這個表達式當成字面意思理解。例如,使用引號標記下表這兩個句子:Say your name! 和 Say "your name"!這兩句話的意思是不一樣的。下表的表達方式告訴你如何在Lisp裏引用表達式:

 

> (quote (+ 1 2))
(+ 1 2)

    "quote"不是一個函數,它不遵照上邊描述的函數求值規則:它的參數不會被求值。"quote"能夠被稱爲一個special form——特殊是由於在處理參數上它有特殊的規則。須要的時候我還會介紹其餘的special forms。這些基本求值規則與這些special forms(特殊形式)一塊兒組成了Lisp語言的語法。

    "quote"使用頻率如此之高以致於專門開發了一個速記符,就是標記在表達式以前的單引號:

 

> '(+ 1 2)        ;單引號速記符
(+ 1 2)

也就是說,上邊的表達式與(quote (+ 1 2))是相同的。

    一旦你學會如何啓動一個計算機程序,確保指導如何退出程序是個好主意。在大多數Lisp系統裏你能夠經過使用exit函數退出。在提供%的shell提示符的UNIX操做系統裏,使用exit函數應該會返回shell:

 

> (exit)
%

其它操做系統可能會提供額外的退出方式。例如,在Macintosh操做系統裏,你能夠到「文件」菜單裏選擇「退出」選項退出。

練習 2.1

先預測,再上機實驗。

a) (+ 3 5 6)

b) (+ (- 1 2) 3)

c) '(+ 3 5 6)

d) (+ (- (* 2 3) (/ 6 2)) 7)

e) 'x

f) ''x

2.2 初級計算與繪圖

    本節介紹一些在Lisp-Stat裏可用的基本的數值與繪圖操做。

2.2.1 一維概要與圖形

    統計數據常常是由數字組組成的。做爲一個例子,下表展現了Minnea-St地區30年間3月份降水量數據記錄,單位爲英寸。

爲了在Lisp-Stat裏檢測這些數據,咱們能夠list函數將數據從新顯示爲數字列表,表達式以下:

> (list .77 1.74 .81 1.20 1.95 1.20 .47 1.43 3.37 2.20
      3.00 3.09 1.51 2.10 .52 1.62 1.31 .32 .59 .81
      2.81 1.87 1.18 1.35 4.75 2.48 .96 1.89 .90 2.05)

 

 

這些數字必須用空格分開,不是逗號。解釋器容許你跨多行表示數據,它不會求值直到你完成表達式的書寫並鍵入回車爲止。

"mean"函數可用來計算數字列表的平均值。咱們能夠用它和list函數組合在一塊兒,來求取降水量樣本數據的平均值:

> (mean (list .77 1.74 .81 1.20 1.95 1.20 .47 1.43 3.37 2.20
      3.00 3.09 1.51 2.10 .52 1.62 1.31 .32 .59 .81
      2.81 1.87 1.18 1.35 4.75 2.48 .96 1.89 .90 2.05))
1.675

降水量的中位數能夠這麼計算:

> (median (list .77 1.74 .81 1.20 1.95 1.20 .47 1.43 3.37 2.20
      3.00 3.09 1.51 2.10 .52 1.62 1.31 .32 .59 .81
      2.81 1.87 1.18 1.35 4.75 2.48 .96 1.89 .90 2.05))
1.47

    每次咱們想對數據樣本進行統計計算的時候,都不得不敲入30組數字,這固然是一件使人惱火的事情。爲了不作那要的工做,咱們能夠Lisp-Stat的一個特殊形式"def"來給數據列表一個變量名。

> (def precipitation
       (list .77 1.74 .81 1.20 1.95 1.20 .47 1.43 3.37 2.20
             3.00 3.09 1.51 2.10 .52 1.62 1.31 .32 .59 .81
             2.81 1.87 1.18 1.35 4.75 2.48 .96 1.89 .90 2.05))
PRECIPITATION

如今,符號precipitation有一個值與它綁定了,咱們的30個數字的列表能夠這樣得到:

> precipitation
(.77 1.74 .81 1.20 1.95 1.20 .47 1.43 3.37 ...)

如今,咱們再對這些數據記性各類數值描述性統計及容易得多了:

> (mean precipitation)
1.675

> (median precipitation)
1.47

> (standard-deviation precipitation)
1.00062

> (interquartile-range precipitation)
1.145

    函數histogram和boxplot能夠用來得到樣本數據的圖形化展現,見圖2.1和圖2.2。

> (histogram precipitation)
#<Object: 143cd68, prototype = HISTOGRAM-PROTO>

> (boxplot precipitation)
#<Object: 145ef04, prototype = SCATTERPLOT-PROTO>

圖2.1 數據柱狀圖

圖2.2 數據盒圖

這兩個表達式是會將繪圖窗口顯示在顯示器上。這兩個函數返回兩條相似這樣的信息:

#<Object: 143cd68, prototype = HISTOGRAM-PROTO>

這個結果以後會被用來識別包含圖形的窗口,但這會兒你能夠忽略它。

    Lisp-Stat也支持數字列表數據的對應元素相乘操做(點乘)。例如,咱們能夠在降水量數據的每一個元素上加1:

> (+ 1 precipitation)
(1.77 2.74 1.81 2.2 2.95 2.2 1.47 2.4299999999999997 4.37 3.2 4.0 4.09 2.51 3.1 1.52 2.62 2.31 1.32 1.5899999999999999 1.81 3.81 2.87 2.1799999999999997 2.35 5.75 3.48 1.96 2.8899999999999997 1.9 3.05)

或者計算它們的天然對數,

> (log precipitation)
(-0.2613647641344075 0.5538851132264376 -0.21072103131565253 0.1823215567939546 0.6678293725756554 0.1823215567939546 -0.7550225842780328 0.3576744442718159 1.2149127443642704 0.7884573603642703 1.0986122886681098 1.128171090909654 0.412109650826833 0.7419373447293773 -0.6539264674066639 0.4824261492442928 0.2700271372130602 -1.1394342831883648 -0.527632742082372 -0.21072103131565253 1.0331844833456545 0.6259384308664954 0.16551443847757333 0.30010459245033816 1.55814461804655 0.9082585601768908 -0.040821994520255166 0.636576829071551 -0.10536051565782628 0.7178397931503168)

或者平方根:

(sqrt precipitation)
(0.8774964387392122 1.3190905958272918 0.9 1.0954451150103321 1.396424004376894 1.0954451150103321 0.6855654600401044 1.1958260743101399 1.835755975068582 1.4832396974191326 1.7320508075688772 1.7578395831246945 1.2288205727444508 1.449137674618944 0.7211102550927979 1.2727922061357855 1.1445523142259597 0.565685424949238 0.7681145747868608 0.9 1.676305461424021 1.3674794331177345 1.0862780491200215 1.161895003862225 2.179449471770337 1.5748015748023623 0.9797958971132712 1.374772708486752 0.9486832980505138 1.4317821063276353)

    降水量數據的分佈有些右偏,平均值和中位數是分離的。你可能想要試着作一些簡單的變換來看看可否對稱化數據,好比說平方根或者對數。你能夠查看圖形和概要統計來肯定是否這些變換確實致使數據更加對稱化一些。例如,能夠用下邊的表達式計算數據樣本的平方根的均值:

> (mean (sqrt precipitation))
1.2401878297708169

也能夠得到樣本數據的平方根的柱狀圖:

> (histogram (sqrt precipitation))

    boxplot函數還能用來並列地產生兩個或者更多的數據樣本。只須要傳遞給boxplot函數列表型參數,這個列表型參數是由多個列表數據組成的,其中的每條列表數據都是要顯示在統計圖形串口上的數據樣本。讓咱們用這個函數測試一些數據樣本,這些數據來源於危地馬拉的一些社會經濟學團體對血液化學的研究。高收入的城市區域和低收入的鄉村區域樣本的血清總膽固醇數據,能夠用下邊的表達式輸入:

> (def urban (list 206 170 155 155 134 239 234 228 330 284
                   201 241 179 244 200 205 279 227 197 242))
URBAN

> (def rural (list 108 152 129 146 174 194 152 223 231 131
                   142 173 155 220 172 148 143 158 108 136))
RURAL

並列的盒圖能夠用如下表達式表示:

> (boxplot (list urban rural))

結果見圖2.3。

圖2.3 危地馬拉城市與鄉村總膽固醇水平並列盒圖

練習 2.2

先預測,後實踐,並解釋其間的不一樣。

a) (mean (list 1 2 3))

b) (+ (list 1 2 3) 4)

c) (* (list 1 2 3) (list 4 5 6))

d) (+ (list 1 2 3) (list 4 5))

練習 2.3

略。

練習 2.4

略。

2.2.2 二維繪圖

    不少單同樣本是隨着時間的推移收集的。上面用到的降水量數據就是一個例子。在一些狀況下假設觀測數據之間相互獨立式合理的,但在有的狀況下是不合理的。檢查數據以得到序列相關性和趨勢的一種方式,就是繪製觀測數據相對於時間,或者相對於數據自身順序的圖形。咱們可使用plot-points函數來產生一幅降水量數據相對於時間的散點圖。plot-points函數使用下邊的表達式形式調用:(plot-points <x-variable> <y-variable>)。其中,<y-variable>能夠是precipitation,就是咱們之前定義的變量;<x-variable>咱們可使用從1到30的天然數序列,咱們能夠本身手動輸入,但有一個更簡單的辦法,使用函數integer-sequence的簡寫形式iseq來生成兩個指定值之間的連續天然數列表。調用iseq函數的通常形式是這樣的: (iseq <start> <end>)。爲了生成咱們須要的序列,咱們使用:

 

> (iseq 1 30)

那麼爲了生成散點圖,咱們須要鍵入如下表達式:

 

> (plot-points (iseq 1 30) precipitation)

結果見圖2.4。對數據來講,本圖彷佛過多的形式,獨立性的假設對該數據來講多是合理的。

圖2.4 降水量水平相對於時間的散點圖

    有時用線將點鏈接起來,就很容易在圖形裏看到時間形式。使用plot-lines函數能夠構造一個鏈接點數據的鏈接線。plot-lines也能夠被用來構造函數的曲線圖。假設你想獲得自變量在-π到+π範圍內的sin(x)函數的圖形,其中常量π由變量pi預約義,你可使用這個表達式:(rseq <a> <b> <n>)來構建數字a與b之間的n等間隔實數。那麼爲了繪製間隔數爲50的等間隔點的sin(x)函數的圖像,能夠這樣:

 

> (plot-lines (rseq (- pi) pi 50) (sin (rseq (- pi) pi 50)))

你也能夠先定義一個變量x來簡化這個表達式,

 

> (def x (rseq (- pi) pi 50)

而後從新構造表達式:

 

> (plot-lines x (sin x))

最終圖形見圖2.5。

圖2.5 sin(x)圖形

    散點圖在檢測兩組數值觀測量之間的關係上,固然是很是有用的,這兩組觀測量是針對同一目標的。一項針對機動車尾氣排放過程的研究給出了46兩汽車的HC和CO排放數據,結果放在兩個變量hc和co裏:

 

> (def hc (list .50 .65 .46 .41 .41 .39 .44 .55 .72 .64 .83 .38
                .38 .50 .60 .73 .83 .57 .34 .41 .37 1.02 .87 1.10
                .65 .43 .48 .41 .51 .41 .47 .52 .56 .70 .51 .52
                .57 .51 .36 .48 .52 .61 .58 .46 .47 .55))

> (def co (list 5.01 14.67 8.60 4.42 4.95 7.24 7.51 12.30
                14.59 7.98 11.53 4.10 5.21 12.10 9.62 14.97
                15.13 5.04 3.95 3.38 4.12 23.53 19.00 22.92
                11.20 3.81 3.45 1.85 4.10 2.26 4.74 4.29
                5.36 14.83 5.69 6.35 6.02 5.79 2.03 4.62
                6.78 8.43 6.02 3.99 5.22 7.47))

而後,能夠用plot-points函數將它們一一對應地繪製到一幅圖裏:

> (plot-points hc co)

結果見圖2.6。

圖2.6 機動車尾氣排放的CO相對HC圖形

練習2.6 

略。

2.2.3 繪圖函數

    上一節中,正弦函數的製圖函數顯得很是笨重。做爲一個代替品,咱們可使用plot-function來繪製有一個參數的函數,該參數有一個指定的範圍。爲了產生圖2.5中的圖形,咱們可使用表達式:

> (plot-function (function sin) (- pi) pi)

    表達式(function sin)用來提取與符號sin相關聯的Lisp函數,僅僅使用sin還不夠,緣由是Lisp裏的符號多是由def設定的數值,同時也多是一個函數定義。起初這一特性看起來有些奇怪,可是它也有一個很大的優點——鍵入一個像這樣的合法表達式:

> (def list '(2 3 4))

而不會銷燬list函數。

    從一個符號裏提取一個函數定義就像引用表達式同樣常見。因此又有一個新的可用的縮寫符號,表達式#'sin和表達式(function sin)是等價的(注:這點與Common Lisp裏類似,有興趣的朋友能夠參考Common Lisp裏的function和apply的用法)。#'的讀音是sharp-quote。使用這個縮寫,產生正弦圖形的表達式能夠寫成這樣:

> (plot-function #'sin (- pi) pi)

2.3 進一步瞭解解釋器

    Lisp-Stat解釋器容許保存建立的變量,能夠部分或所有地記錄與解釋器的會話。另外,它還提供了得到最近幾條計算結果的機制,還有得到在線幫助的機制。這幾個特性可能根據所運行的系統的不一樣稍微有些差別,尤爲是在線幫助系統,例如,一些系統可能提供了單獨的幫助窗口。個人相關工做是基於XLIST-STAT這一實現版本的。

2.3.1 保存你的工做

    若是你想要與解釋器的一個會話,你可使用dribble函數。表達式

> (dribble "myfile")

將開始一個記錄。全部你輸入的表達式和解釋器返回的結果都會輸入到這個叫myfile的文件,表達式

> (dribble)

會中止錄製工做。表達式(dribble "myfile")一般會開始一個名稱爲myfile的新文件,若是你已經有了這個名字的文件,其內容將所有丟失。所以你不能用dribble函數對單個文件進行開啓和關閉。

    dribble僅會記錄文本,不會記錄圖形。Macintosh的XLISP-STAT裏你能使用標準Macintosh快捷命令COMMAND-SHIFT-3來保存當前屏幕的圖像。也能夠到Edit菜單裏選擇Copy命令,或者其快捷鍵COMMAND-C,圖形窗口是一個活動窗口,用來將內容保存到剪切板。在X11版的XLISP-STAT,繪圖菜單包含一個選項,該選項能夠將圖形的PostScript版本保存成文件。

    你在Lisp-Stat裏定義的變量只能存在於當前會話,若是你從Lisp-Stat裏退出,或者程序崩潰了,你的數據就會丟失。爲了保存你的數據可使用savevar函數。這個函數容許你將一個或多個變量保存到一個文件裏。新文件建立時,任何同名文件都會被銷燬。爲了將變量precipitation保存到一個名爲precipitation.lsp的文件,鍵入:

> (savevar 'precipitation "precipitation")

不須要加.lsp後綴名,savevar會幫你加上的。爲了保存兩個變量hc和co到examples.lsp文件,鍵入:

(savevar '(hc co) "samples")

我使用一個帶引號的列表'(hc co)並將該符號列表傳遞給savevar函數。還有一個長點的表達式能夠代替它:(list 'hc 'co)。

    文件precipitation.lsp和samples.lsp得到了表達式集合,當使用read函數讀取這兩個文件的時候,它們將從新建立precipitation、hc和co變量。爲了加載precipitation.lsp文件,可使用以下表達式:

> (load "precipitation")

2.3.2 命令歷史機制

    Common Lisp提供了一個簡單的命令歷史機制。符號 -, *, **, ***, +, ++和+++就是爲實現這一意圖的。解釋器將這些符號與其意義綁定以下:

  • -             當前輸入的表達式
  • +            讀取最近一個表達式
  • ++          +返回的表達式的前一個表達式
  • +++        ++返回的表達式的前一個表達式
  • *             最近的一個表達式的結果
  • **            *所屬表達式的前一個表達式的結果
  • ***           **所屬表達式的前一個表達式的結果

    *, **, ***這三個變量也許是最有用的。例如,若是你想求某個變量的對數值,但忘了給它定義變量名,你就可使用那幾個歷史變量來避免重複計算對數:

> (log precipitation)
(-0.2613647641344075 0.5538851132264376 -0.21072103131565253 0.1823215567939546 0.6678293725756554 0.1823215567939546 -0.7550225842780328 0.3576744442718159 1.2149127443642704 0.7884573603642703 1.0986122886681098 1.128171090909654 0.412109650826833 0.7419373447293773 -0.6539264674066639 0.4824261492442928 0.2700271372130602 -1.1394342831883648 -0.527632742082372 -0.21072103131565253 1.0331844833456545 0.6259384308664954 0.16551443847757333 0.30010459245033816 1.55814461804655 0.9082585601768908 -0.040821994520255166 0.636576829071551 -0.10536051565782628 0.7178397931503168)
> (def log-precip *)

如今,變量log-precip就包含降水量數值的對數了。

    系統利用了2.2.3節討論的一個事實,那就是符號能夠同時又函數定義和數值。符號*的函數定義是乘法函數,它的值是解釋器返回的上一個求值結果。

2.3.3 獲取幫助

Lisp-Stat提供了大量的不一樣的函數。咱們不可能精確地記起每一個函數的用法。一個交互式的幫助應用能夠完成一些任務:找出須要的正確的函數,找出如何更容易地使用這個函數。例如,下邊就是如何得到median函數幫助的例子:

> (help 'median)
loading in help file information - this will take a minute ...done
MEDIAN                                                          [function-doc]
Args: (x)
Returns the median of the elements of X.
NIL
>

median函數前面的引號是必要的。help函數自己也是函數,它的參數是median函數的符號形式。爲了確保help函數可以接受到符號,而不是符號的值,你須要引用(quote)這個符號(symbol)。

    若是你對函數的名字不肯定,你可能須要使用help*函數來獲取幫助。假設你想找出與正太分佈相關的函數,它們的名字裏多數都包含"norm",表達式(help* 'norm)將打印函數名裏包含"norm"的全部符號的幫助信息:

(help* 'norm)
------------------------------------------------------------------------------
BIVNORM-CDF                                                     [function-doc]
Args: (x y r)
Returns the value of the standard bivariate normal distribution function 
with correlation R at (X, Y). Vectorized.
------------------------------------------------------------------------------
Sorry, no help available on NORM
------------------------------------------------------------------------------
Sorry, no help available on NORMAL
------------------------------------------------------------------------------
NORMAL-CDF                                                      [function-doc]
Args: (x)
Returns the value of the standard normal distribution function at X.
Vectorized.
------------------------------------------------------------------------------
NORMAL-DENS                                                     [function-doc]
Args: (x)
Returns the density at X of the standard normal distribution. Vectorized.
------------------------------------------------------------------------------
NORMAL-QUANT                                                    [function-doc]
Args (p)
Returns the P-th quantile of the standard normal distribution. Vectorized.
------------------------------------------------------------------------------
NORMAL-RAND                                                     [function-doc]
Args: (n)
Returns a list of N standard normal random numbers. Vectorized.
------------------------------------------------------------------------------
Sorry, no help available on NORMALREG-MODEL
------------------------------------------------------------------------------
Sorry, no help available on NORMALREG-PROTO
------------------------------------------------------------------------------
NIL
>

符號norm, normal, normalreg-model和normalreg-proto沒有函數定義,所以沒有可用的幫助信息。函數文檔裏的術語vectorized意味着這個函數能夠被應用到形式爲列表的參數上;它的結果是這個函數做用到參數的每個元素後返回的結果列表。文檔中出現的一個相關的術語是vector reducing。若是一個函數能夠遞歸地操做它的參數,直到只剩一個單獨的數字爲止,那麼這個函數就是vector reducing(矢量減小)的。函數sum, prod, max和min都是矢量減小的。

    做爲使用help*的替代品,你可使用apropos函數來得到包含字符串"norm"做爲其名字一部分的那些符號的列表:

> (apropos 'norm)
NORMALREG-PROTO
NORMALREG-MODEL
NORM
BIVNORM-CDF
NORMAL-QUANT
NORMAL-RAND
NORMAL
NORMAL-CDF
NORMAL-DENS
>

而後使用help函數詢問每一個符號的更多信息。

    讓我來簡要地介紹一下help函數打印的信息所用到的表示方法,來描述一個函數指望參數有的功能。表示方法對應的Lisp函數定義的參數列表的規格將在4.2節描述。多數函數指望參數列表是固定的,想在幫助信息裏的這行同樣: Args: (x y z)。有的參數可能會接收一個或多個可選參數,這種函數的參數可能會這樣表示: Args: (x &optional y (z t)),這意味着x是必選參數,y和z是可選參數。若是這個函數命名是f,能夠這樣調用:(f x-val), (f x-val y-val),或者(f x-val y-val z-val)。列表(z t)意味着沒有提供參數z,那它的默認值就是t;y是沒有強制指定默認值的參數,所以它的默認值是nil。實參必須按照形參規定的順序和規則來提供,所以若是你要是想傳遞參數z,那就必須先對y賦值。

    另外一個形式的可選參數是關鍵字參數。例如,histogram函數的參數是這樣的:Args: (data &key (title "Histogram")),data參數是必選的,title就是可選的關鍵字參數,其默認值是"Histogram"。若是你用2.2.1節用到的降水量數據建立一張表,它的名字是"Precipitation"能夠用這個表達式:

> (histogram precipition :title "Prcipitation")

圖2.7 圖形title修改示意圖

所以,爲了給關鍵字參數一個值,你將須要顯式地給出這個關鍵字名,也就是一個由冒號緊跟着參數名的一個符號,而後是這個參數的值。若是一個函數帶多個關鍵字參數,在調用時這些關鍵字參數之間的順序能夠任意指定,前提是他們需呀在必選參數可可選函數以後。(注:關於這點建議參考《實用Common Lisp編程》第4.2~4.5節,認真實踐,方可領悟!)

    最後,一些函數可能帶不定數量的參數,可用下式表示:Args: (x &rest ars)。該參數列表表示參數x是必選參數,同時能夠帶0個到多個附加參數。(注:至關於C語言參數列表裏的 「...」)

    除了函數,help函數還能針對數據類型和一些變量給出幫助信息。例如:

> (help 'complex)
loading in help file information - this will take a minute ...done
COMPLEX                                                         [function-doc]
Args: (realpart &optional (imagpart 0))
Returns a complex number with the given real and imaginary parts.
COMPLEX                                                             [type-doc]
A complex number
NIL

這個例子展現了complex複數數據類型的功能和類型信息。下邊這個例子

> (help 'pi)
PI                                                              [variable-doc]
The floating-point number that is approximately equal to the ratio of the
circumference of a circle to its diameter.
NIL

展現了pi這個變量的文檔。

2.3.4 列舉變量和取消變量定義

    在工做一段時間以後,你可能想要查明使用def宏定義了哪些變量,variables函數將以列表的方式返回這些變量:

> (variables)
(CO HC PRECIPITATION RURAL URBAN)

    你可能偶爾對不在使用的變量進行釋放以得到一些空間。要達到這個目的你可使用undef函數:

> (undef 'co)
CO

> (variables)
(HC PreCIPITATION RURAL URBAN)

undef的參數是一個符號(symbol),本例總咱們必須引用它以免解釋器對其求值。undef也能接受符號列表當作參數,所以想要取消定義當前全部已經定義的變量,能夠這樣:

> (undef (variables))
(HC PRECIPITATION RURAL URBAN)
> (variables)
NIL

2.3.5 中斷計算

偶爾你可能在系統裏開始了這樣一種計算,它耗用了太長的時間或者產生了一個看起來不合理的結果。每一個系統都提供了一個機制來中斷這樣的計算過程,可是這些機制依據系統的不一樣而不一樣。在Macintosh操做系統下的XLISP-STAT裏,你能夠按住COMMON鍵和PERIOD鍵直到系統將解釋器交由你來輸入。在UNIX版本的XLISP-STAT裏,你能夠按下爲你的系統準備的中斷鍵組合來中斷這種計算,一般這個中斷鍵組合爲CONTROL+C。

2.4 一些數據處理函數

    目前爲止咱們都是將數據直接敲到解釋器裏的,數據集較小,因此這也不是什麼大問題,然而,有其它的替代方法會更有用。Lisp-Stat提供了一些函數,能夠生成系統的數據和隨機數據,能夠提取數據集的一部分,修改數據集,而且能夠從文件中讀取數據。

2.4.1 生成有系統的數據

    函數iseq用來生成連續的整數序列,rseq用來生成等間隔的實數序列,它們在2.2.2節提到了。函數iseq也能夠經過傳遞一個單獨參數來調用,這種狀況下參數應該是一個正整數n,結果是從0到n-1的整數列表。例如:

> (iseq 10)
(0 1 2 3 4 5 6 7 8 9)

    函數repeat在生成具備必定格式的序列時是頗有用的,調用repeat函數的通常格式是(repeat <vals> <pattern>),參數<vals>多是簡單的數據項或者數據項列表,參數<pattern>必須是一個單獨的正整數或者一組正整數列表。若是<pattern>是一個單獨的正整數,repeat就會將<vals>重複<pattern>次,例如:

> (repeat 2 3)
(2 2 2)

> (repeat (list 1 2 3) 2)
(1 2 3 1 2 3)

若是<pattern>是一個列表,那麼<vals>應該是一個與<pattern>相同長度的列表,本例中,<vals>和<pattern>兩個列表中相同位置的元素市一一對應的,<pattern>中每一個元素的值對應<vals>列表中相對元素重複的次數。例如:

> (repeat (list 1 2 3) (list 3 2 1))
(1 1 1 2 2 3)

    函數repeat在設計實驗中的編碼處理層面上是極其有用的。例如,在一個測試植物種植密度對馬鈴薯生長的影響的實驗裏,有3個品種的馬鈴薯,種植在4個不一樣種植密度的地方,每一種測試組合重複3次,結果在下表中:

咱們能夠將這些數據按行,輸入到產量變量裏,如:

> (def yield (list  9.2 12.4  5.0  8.9  9.2  6.0 16.3 15.2  9.4
                    12.4 14.5  8.6 12.7 14.0 12.3 18.2 18.0 16.9
                    12.9 16.4 12.1 14.6 16.0 14.7 20.8 20.6 18.7
                    10.9 14.3  9.2 12.6 13.0 13.0 18.3 16.0 13.0))

使用repeat,咱們能夠這樣生成表示不一樣品種和密度水平的變量:

> (def variety (repeat (repeat (list 1 2 3) (list 3 3 3)) 4))
VARIETY
> (def density (repeat (list 1 2 3 4) (list 9 9 9 9)))
DENSITY
>

練習2.7

略。

練習2.8

略。

2.4.2 生成隨機數據

    一些函數能夠生成僞隨機數字,例如,表達式(uniform-random 50)將生成一個的數字列表,列表的元素是50個在0和1之間的獨立隨機均勻分佈的數字。函數normal-rand和cauchy-rand也相似地生成標準正態的或標準柯西隨機變量。表達式:

> (gamma-rand 50 4)
(1.65463390351748 2.8773826027572174 3.324805876512575 3.7107737590852885 2.5322960063057325 2.3917079623070143 2.5738732854681405 5.114579884756656 4.503560922291987 3.8212691381896735 4.748609795710945 2.460930129538628 5.1583129719078125 2.8722892003673204 2.271632493827859 1.306689218830528 2.956587212984247 2.254496336799741 1.5432493948171886 1.0074517838362944 2.3221187903990685 2.9297357635208368 1.9111207190664348 4.086600231358658 8.498775494298609 3.309145002227253 8.43403656671126 2.9001355594385445 2.156028887326921 5.457111466807575 2.7180957670090233 6.727592025517272 3.2182935800579466 4.759355290676398 2.1463912302657 4.641326968027183 4.951716760175214 2.9685374305637593 5.934043534710339 6.847560103820065 5.612014958257557 4.4811213649680814 7.4465265497488975 2.6911820555097306 1.0087771420430949 1.4549088997831112 5.276729335543616 2.773807551655783 3.565519147704477 2.7611738530317305)

生成一個大小爲50的數據樣本,這些數據符合伽馬-分佈,單位尺度,指數爲4。

指數爲3.5~7.2的β分佈變量能夠這樣生成:

> (beta-rand 50 3.5 7.2)
(0.4565276604553964 0.2784641511528651 0.1992946182796908 0.15531331094658482 0.1524256312642991 0.26070619189857414 0.26797467262440017 0.22086809762954612 0.30818508575346343 0.31431320706862537 0.4767943879658169 0.3103165584027092 0.3099168733453603 0.48844354838756315 0.3850748530401529 0.5261394345014925 0.4729101627118159 0.34008228819566094 0.34728211784541685 0.11375668596438508 0.3177405260351405 0.15233790451255871 0.3247481923221466 0.16163693218582695 0.34148249992894264 0.26747115750175693 0.1917463536194362 0.3096184276126184 0.2815407126874605 0.27277307887112906 0.5364074613555034 0.5343432853990229 0.7113134763751001 0.2008014529296493 0.1610644249288089 0.14281693672275156 0.1994676889199915 0.3627351688059498 0.4983412412275205 0.4279902846854145 0.42865178072445864 0.32624949341033294 0.38525682966664915 0.3395031181658098 0.35242269440634416 0.3010577148339438 0.3497046925914022 0.12251980989006983 0.3635843569985974 0.12912742911247255)

大小爲50,自由度爲4的t-分佈隨機變量樣本:

> (t-rand 50 4)

自由度爲4的卡方-分佈隨機變量的生成表達式:

> (chisq-rand 50 4)

分子自由度爲3,分母自由度爲12的F-分佈的隨機變量的生成表達式:

> (f-rand 50 3 12)

    參數n=5,p=0.3的二項式分佈樣本生成表達式爲:

> (binomial-rand 50 5 .3)

大小爲50,指望值λ=4.3的泊松-分佈隨機數據樣本生成表達式:

> (poisson-rand 50 4.3)

    這些樣本生成函數的「大小」這一參數也能夠是一個整數列表,相應地,返回的結果是樣本列表的列表(即多條樣本列表)。例如:

> (normal-rand (list 10 10 10))
((0.24711492524397788 -0.08672761627290786 -0.3586542916304954 -1.2370666723700587 0.4973152707135397 1.2779437429528906 0.5433681505645361 -0.15811161781047764 -2.444951338731246 1.4391452075883056) 
(-0.024198159865474137 0.018494098689563702 -1.201497046490873 0.29215227642978814 -0.1504967544867967 0.8539957140604427 0.1866809359842353 -0.3235190277039044 1.2316501558225053 -0.3173822473522208) 
(2.6816999640212815 1.312819856883247 0.3381892163767421 0.20384488399908143 0.9360472535774436 -1.3160001465807794 0.033750173542764696 -0.4613961332004225 -0.9523674256520025 -2.2724093805908905))

該表達式產生了一個列表集,其中的每一個列表都是有十個正態分佈變量的列表。

    最後,你可使用sample函數從一個列表裏隨機裏選出一個數據樣本,表達式(sample (iseq 1 20) 5)的返回值就是,在1, 2, ..., 20這個列表裏,隨機地選擇5個元素,組成一個大小是5的新的沒有重複元素的樣本,其結果多是:

> (sample (iseq 1 20) 5)
(2 19 1 3 5)

sample函數能夠經過啓用第三個參數(一個可選參數)來指定採樣元素是否能夠重複。那麼下個表達式

> (sample (iseq 1 20) 5 t)
(11 11 5 4 16)

將容許sample函數在採樣的時候容許元素重複(注:原文使用replacement表示的,但我沒法翻譯成合適的中文,若是有更精確的翻譯,歡迎指正)。

第三個參數爲nil時,sample函數採樣時不容許元素重複,這與忽略這個參數的狀況是同樣的。

練習2.9

略。

2.4.3 創建子集和刪除

select函數容許你從一個列表裏選擇一個單獨的元素或者一組元素。例如,你這樣定義了一個x變量:

> (def x (list 3 7 5 9 12 3 14 2))

那麼

> (select x i)

將返回變量x的第i個元素。Lisp語言列表的下標是從0開始的,這點與C語言相同,可是和FORTRAN不一樣。所以它的下標應該是0, 1, 2, 3, 4, 5, 6, 7,因此有:

> (select x 0)
3
> (select x 2)
5

爲了獲取一組元素,咱們可使用下標數組來代替單個的下標,即:

> (select x (list 0 2))
(3 5)

    若是你想選擇除了下標爲2的元素以外的全部元素的話,使用如下表達式

> (remove 2 (iseq 8))

來產生下標序列,並把它做爲select函數的第二個參數:

> (select x (remove 2 (iseq 8)))
(3 7 9 12 3 14 2)

表達式(iseq 8)生成了從0到7的整數列表。remove函數會返回一個移除指定元素的列表,那麼

> (remove 3 (list 1 3 4 3 2 5))
(1 4 2 5)

剩下的元素按它們在原始列表裏的順序返回。

    另外一個返回除了2號下標元素的方法是使用對照函數"/="(意思是「不等於」)來告訴你哪些下標不等於2:

> (/= 2 (iseq 8))
(T T NIL T T T T T)

函數which可以返回,其函數的參數的不是nil的哪些下標的元素的列表,即

> (which (/= 2 (iseq 8)))
(0 1 3 4 5 6 7)

結果能夠傳遞給select函數:

> (select x (which (/= 2 (iseq 8))))
(3 7 9 12 3 14 2)

對於刪除一個元素來講,這個方法有一點太笨重了,然而它倒是一個通用的方法。例如以下表達式

> (select x (which (< 3 x)))
(5 7 9 12 14)

將返回全部比3大的元素。

    因爲咱們在上邊幾個例子裏用到的x變量比較短,經過數數咱們也能夠很容易地肯定他的長度。對於再長點兒的列表就不太容易了,咱們可使用length函數來實現。使用length函數,咱們也能實現"選擇除了下邊爲2的其它元素"這一功能:

> (select x (remove 2 (iseq (length x))))

練習2.10

略。

2.4.4 合併列表

    有時你可能想將幾個短的列表合併爲一個長的列表,能夠用append和combine函數來完成。例如,

> (append (list 1 2 3) (list 4) (list 5 6 7 8))
(1 2 3 4 5 6 7 8)

或者

> (combine (list 1 2 3) (list 4) (list 5 6 7 8))
(1 2 3 4 5 6 7 8)

這兩個函數的不一樣在於它們對列表的列表的處理方式上。append函數只合並外部列表,而combine函數會遞歸地調用到子列表,並返回一個沒有複合元素的列表(即不包含子列表):

> (append (list (list 1 2) 3) (list 4 5))
((1 2) 3 4 5)

> (combine (list (list 1 2) 3) (list 4 5))
(1 2 3 4 5)

combine函數容許他的參數是數字或者其它簡單數據項,好比:

> (combine 1 2 (list 3 4))
(1 2 3 4)

而append函數的參數必須是列表(list)形式,不然會報錯。

> (append 1 2 (list 3 4))
Error: bad argument type - 1
Happened in: #<Subr-APPEND: #14644cc>

2.4.5 修改數據

    setf能夠用來修改列表的單個元素或者列表的元素集合。再一次假設咱們設置了一個列表x:

> (def x (list 3 7 5 9 12 3 14 2))

咱們想把12改爲11,使用以下表達式(setf (select x 4) 11)來作這個替換:

> (setf (select x 4) 11)
11
> x
(3 7 5 9 11 3 14 2)

    setf的通常形式爲(setf <form> <value>),這裏<form>是你想改變的列表裏的那個元素或元素組,<value>是你想要改變成的那個值或多個值得列表。那麼表達式(setf (select x (list 0 2)) (list 15 16))將下標爲0和2的元素的值改爲15和16,即:

> (setf (select x (list 0 2)) (list 15 16))
(15 16)
> x
(15 7 16 9 11 3 14 2)

    這裏須要有一個小提示,Lisp符號只不過是不一樣項的標籤而已,當咱們使用def命令爲一個項分配名字時,咱們並無產生一個新項,所以,表達式(def x (list 1 2 3 4))求值以後,在進行(def y x),符號x和符號y是針對相同列表的兩個不一樣的名字。結果,咱們使用setf對x的元素進行操做時,咱們同時也改變了y的值,由於x和y在引用一個相同的列表。若是你想作一份x的拷貝,而且在x改變以前將之保存到y裏,那麼須要你強制執行,此時須要使用copy-list函數。表達式(def y (copy-list x))製做了x的一份拷貝,並將它設置給y變量。如今x與一引用了不一樣的項,改變x將不會影響到y的值。

    setf也能夠用來將一個值分配給一個符號,好比:

> (setf z 3)
3
> z
3

使用def而不使用setf來定義變量是有優點的,那就是def會記錄它定義過的變量的名稱,並容許使用variables函數找出可用的變量。setf還有一些其餘用法,那些咱們等到3.8節再講。

練習 2.11

略。

2.4.6 讀取數據文件

兩個可用的讀取數據的文件是read-data-columns和read-data-file。

第一個是read-data-columns,是被設計來讀取包含一些數據列的文件的。它帶兩個參數:表示文件名的字符串和一個整數,表示文件裏的列數。它返回一個列表集的列表,其中每一列都放到一個列表裏。數據裏的每項能夠由任意數目的空格分開,可是不能用逗號和其它標點符號。每項多是任何合法的Lisp表達式,特別地,他們能夠是數字、字符串或符號。全部的列必須是等長的。返回的結果列表裏的內容均可以使用select函數提取。

    假設咱們有一個叫作mydata的文件,包含兩列數字數據。咱們想要讀取這些數據並將數據列分別命名爲x和y,咱們能夠首先讀取整個文件並將結果保存到名爲mydata的變量裏:

> (def mydata (read-data-columns "mydata" 2))

如今,變量mydata的值是一個包含兩個數字列表的列表。爲了提取和命名這兩個數據列表,咱們能夠用以下表達式:

> (def x (select mydata 0))

> (def y (select mydata 1))

read-data-columns的第二個參數是可選的。若是沒提供,那麼數據的列數將由從第一個非空行開始的全部數目肯定。

    在有的系統上,文件名參數也多是可選的,若是忽略了該參數,將出現一個對話框以選擇要讀取的文件。

    若是數據文件不是按列形式組織的,你可使用read-data-file函數讀取文件的全部內容,並以單獨一個列表的形式返回。該列表可使用Lisp-Stat函數操做成爲合適的形式。文件條目是爲read-data-columns準備的,函數read-data-file帶單一參數,是一個要讀取的文件名的字符串,有的系統裏你能夠忽略這個參數,並使用對話框選擇要讀取的文件。

    這兩個函數可以知足大多數的需求,若是你必需要讀取非指定格式的文件,你可使用lisp的原始的文件處理函數,這些將在第4.5節描述。

2.5 動態繪圖

    在第2.2節裏,咱們使用直方圖和散點圖來測試單變量和兩個變量之間關係的分佈狀況。本節咱們將介紹一些動態繪圖方法和能夠探頭3個或多個變量之間關係的繪圖技術。前4個小節的繪圖技術和方法很簡單,僅須要一些鼠標操做和簡單的命令。第5小節將介紹繪圖對象的概念,還有從解釋器向一個對象發送消息的方法。這些想法是用來向散點圖裏加入直線和曲線。對象和消息的基本概念在2.6節迴歸模型一節裏還要用到。最後一小節將展現如何使用Lisp迭代器來構造動態模擬動畫,一副運動的圖形。

2.5.1 旋轉圖

    在探索多個變量之間的關係的時候,咱們很天然會想到構造一個數據的三維散點圖。固然,咱們只能在顯示屏上看到一個該圖形的二維投影,任何關於該數據的你能經過線性模型感知到的圖形深度信息都會丟失。一個試圖恢復圖形深度感知的方法就是繞着某一座標軸旋轉這些點。spin-plot函數容許你構造一個可以旋轉的三維圖形。

    舉個例子,讓咱們看一些收集好的數據,來檢測一下材料損耗量(當橡膠樣本與研磨材料一塊兒摩擦式發生的)與標本兩個屬性之間的關係,硬度和抗拉強度。 硬度用邵氏度來度量,抗拉強度用kg/cm²度量,研磨損失用g/馬力小時度量。數據能夠這樣輸入:

> (def hardness
       (list 45 55 61 66 71 71 81 86 53 60 64 68 79 81 56 68 75
             83 88 59 71 80 82 89 51 59 65 74 81 86))
HARDNESS
> (def tensile-strength
       (list 162 233 232 231 231 237 224 219 203 189 210 210 196
             180 200 173 188 161 119 161 151 165 151 128 161 146
             148 144 134 127))
TENSILE-STRENGTH
> (def abrasion-loss
       (list 372 206 175 154 136 112  55  45 221 166 164 113  82
              32 228 196 128  97  64 249 219 186 155 114 341 340
             283 267 215 148))
ABRASION-LOSS

表達式

> (spin-plot (list hardness tensile-strength abrasion-loss))
#<Object: 136b4f8, prototype = SPIN-PROTO>

將產生一幅圖形。

圖2.7 磨損數據的旋轉圖

spin-plot函數的參數是一個含有3個list或vector類型的list,表示x, y, z三個變量。初始條件下z軸垂直屏幕並指向屏幕外部。你能夠將鼠標放到Pitch、Roll和Yaw前的小方框裏,而後點擊鼠標來旋轉圖形。經過旋轉圖形,你能夠看到全部點幾乎落到同一個平面上。

圖2.8 磨損數據的第二個視圖

圖2.8展現的數據基本在一個平面上。所以,線性模型應該能提供合理的修正。事實上,在執行一段時間的旋轉操做以後,你會注意到數據點會有輕微的扭曲,暗示了存在一些二階項成分。不幸的是,在打印圖裏這部分信息很難傳達給用戶。

    spin-plot函數的一些選項能夠被指定爲關鍵字參數,給:title關鍵字傳遞一個字符串就能夠爲繪圖窗體制定一個可替換的標題,關鍵字:variable-labels能夠用來提供一個能夠替換的座標軸標籤,如使用一個字符串列表替換x、y和z。例如,對於磨損數據你可使用一下表達式來構建圖形來代替上邊咱們給出的表達式:

> (spin-plot (list hardness tensile-strength abrasion-loss)
             :title "Abrasion Loss Data"
             :variable-labels (list "Hardness"
                                    "Tensile Strength"
                                    "Abrasion Loss"))
#<Object: 142b29c, prototype = SPIN-PROTO>

函數spin-plot還能夠接收關鍵字參數:scale,若是:scale賦值爲t(t也是其默認值),那麼數據將三個變量的均值進行中心化,這三個變量都將被尺度化以適應圖形尺寸;若是:scale賦值爲nil,數據將被假定尺度化爲-1到1之間,圖形將繞原點旋轉。若是你想要以變量的均值來中心化你的圖形,而且用相同的量尺度化全部的觀測值的話,可使用該表達式:

> (spin-plot 
   (list (/ (- hardness (mean hardness)) 140)
         (/ (- tensile-strength (mean tensile-strength)) 140)
         (/ (- abrasion-loss (mean abrasion-loss)) 140))
   :title "Abrasion Loss Data"
   :variable-labels (list "Hardness"
                          "Tensile Strength"
                          "Abrasion Loss")
   :scale nil)
#<Object: 145f7c0, prototype = SPIN-PROTO>

    每個Lisp-Stat繪圖窗口都提供了一個菜單欄,用來與圖形進行通訊。使菜單 精確的可用的方式依賴於窗體管理系統。在Macintosh操做系統的XLISP-STAT裏,當圖形窗口是前置窗口時,圖形的菜單是放在菜單欄裏的;在XLISP-STAT的X11版本里,窗體上有一個按鈕來彈出菜單。可旋轉的圖形上的菜單容許你改變速度與角度,或者是否使用深度感知,是否顯示座標軸等。默認狀況下,圖形是使用深度感知的:離觀察者近的點總比離觀察者遠的點繪製的大些。

    大多數窗口系統都有提供一個修改鼠標點擊的慣例,目的是區別本次點擊是想要開始一次文本選擇作點擊,仍是想要擴展已存在的選擇而作的點擊。Macintosh操做系統的處理方式是點擊鼠標的時候按下shift鍵,來發出一個擴展信號。換句話說,shift鍵就是Macintosh的擴展修改器。在旋轉的圖形裏,不使用擴展修改器的話,在Pitch、Roll或者Yaw方框裏點擊鼠標將引發圖形的旋轉。在合適的擴展修改器上一直按着鼠標不放將引發圖形不停地旋轉甚至是釋放鼠標釋放以後,旋轉將在你再次點擊旋轉控制面板的時候纔會中止(注:經本人測試,該操做在Windows操做系統的XLISP-STAT一樣適用)。

練習2.12

略。

練習2.13

略。

2.5.2 散點圖矩陣

另外一個繪製變量集的方法是查看變量間全部可能散點圖組合矩陣。scatterplot-matrix函數將提供這樣一幅圖形。使用上節提到的磨損數據,

> (scatterplot-matrix
   (list hardness tensile-strength abrasion-loss)
   :variable-labels
   (list "Hardness" "Tensile Strength" "Abrasion Loss"))
#<Object: 142b78c, prototype = SCATMAT-PROTO>

表達式產生散點圖矩陣,如圖2.9。

圖2.9 磨損數據的散點圖矩陣

    "磨損"相對"抗拉強度"那幅圖形給了咱們關於這兩個變量的聯合變量的想法,可是此時"硬度"這個變量也隨着點的不一樣而不一樣。爲了得到這3個變量之間關係的理解,咱們最好將「硬度」這個變量固定在一個可變的水平,而後觀察當咱們改變硬度這個變量時,"磨損"和"抗拉強度"這兩個變量之間的關係是怎樣變化的。經過使用兩個高亮技術:selecting和brushing,你能夠在散點圖矩陣裏進行諸如此類的探索。

    選擇模式. 當鼠標時一個箭頭的時候,表示你的圖形處於選擇模式,這是默認的設置。在這個模式裏你能夠移動箭頭鼠標到一個點上而後點擊鼠標來選擇一個點。爲了選擇一組點,能夠在目標點上拖出一個選擇框。若是某些點不能包含到選擇框了,你能夠經過「拖選擴展」操做加入那些點,即按住shift鍵再拖選或點選(Macintosh和Windows操做系統一樣適用)。若是你不用擴展修改器(即shift鍵),任何已選的元素都會變爲非選擇狀態;當你適用擴展修改器時,應經選擇的點仍然是選擇狀態。

    刷模式.你能夠經過繪圖菜單的Mouse Mode項,在彈出的對話框裏選擇Brushing來進入刷模式。在這個模式裏,光標看起來像一把刷牆的刷子和一個附在它上邊的虛線框。當你在圖形裏移動刷子的時候,在刷子內部的點是高亮的。刷子外的點沒有高亮,除非是已經被選中的。爲了在刷模式下選擇點(讓它們永久高亮),你須要在移動鼠標時保持鼠標按下。

    在圖2.10的統計圖例,「硬度」變量的中間部分的點已經被用一把長且窄的刷子高亮了。經過使用繪圖菜單的Resize Brush命令,能夠改變刷子的大小。在磨損對抗拉強度那幅圖裏,你能夠看到高亮的點看起來是沿着一條曲線分佈的。若是你想擬合一個模型到這些數據上,本次觀測結果建議根據曲率來擬合模型。這與觀測這些數據旋轉獲得的結論是類似的。

    經過選擇而高亮數據在旋轉的圖形中是有用的。經過選定一個變量的一小片具備類似值得數據點,而後旋轉它,你能夠檢測給定的一個固定變量狀況下,其它的兩個變量的條件分佈的三維結構。

    散點圖矩陣在檢測一個定量變量和幾個分類變量之間的關係時是頗有用的。在第2.4.1節裏咱們看到了一些關於馬鈴薯種植的實驗數據,這些數據用到了3個馬鈴薯品種和4種種植密度。圖2.11展現了這些數據的散點圖。該圖使用一個長且窄的刷子來高亮第三個品種的數據點。若是品種和密度之間沒有相互做用,那麼當刷子從一個品種移動另外一個品種上的時候,高亮點的形狀應該大至是以平行的方式移動的。

與spin-plot函數類似,scatterplot-matrix函數也接受可選的關鍵字參數:title, :variable-labels和:scale。

練習2.14

略。

練習2.15

略。

2.5.3 與單個圖形交互

    旋轉的圖形和散點圖矩陣都是交互性的統計繪圖。簡單的散點圖也容許一些交互。若是你選擇菜單裏的Show Labels選項,挨着高亮點將出現一個標籤。你可使用選擇模式或刷模式來高亮一個點,默認的標籤多是"0", "1", ...這樣的形式,大多數繪圖函數都容許你使用:point-labels關鍵字參數來提供一個可替換的標籤列表。

    在查看大量數據集的時候有一個選項是頗有用的,那就是移除一部分點,你能夠這麼作,選擇要移除的點,而後選擇菜單裏的Remove Selection菜單項。使用菜單裏的Remove Selection能夠完成窗口尺寸的從新調整。

    當選中繪圖窗口裏的點集時,你可使用Selection Symbol菜單項來改變顯示點的符號。若是你的系統支持使用顏色,你可使用Selection Color菜單項來設置選中點的顏色。

    經過繪圖窗口菜單的Selection... 菜單項,你能夠保存一個變量裏當前選中點的指標。而後將彈出對話框讓你輸入一個名字。當沒有選擇點時,你可使用Selection...菜單項來指定隨後選擇的點的指標(注:該功能能夠將任意點集編成一組,並經過Selection...菜單項,使用名稱來選中、高亮。)。將彈出一個對話框讓你輸入一個表達式,用來肯定所選點的指標(注:即所選點的符號),這個表達式能夠是任意Lisp表達式,能夠是單獨指標或者指標列表。(注:我的認爲這裏的指標譯爲符號更容易理解,實際上就是起個名字)

    圖2.12鏈接瞭如下兩種數據,乙醇相對氧氣水平的散點圖和針對發酵數據的糖類型指標直方圖。來自於同一組的糖是高亮的。

2.5.4 鏈接圖

    當你在散點圖矩陣裏刷取或選擇的時候,你能夠看下這組散彈圖的其它圖形的相互做用。經過選擇你想要鏈接的那個圖形的菜單的Link View選項,你能夠構建本身的交互圖集。例如,2.5.2節練習2.14的數據,咱們能夠把乙醇和氧氣放到散點圖裏,把糖放到直方圖裏。若是咱們鏈接這兩個圖形,而後選擇直方圖裏兩個糖數據中的一個,散點圖裏對應的點就會高亮,就像圖2.12中顯示的那樣。

圖 2.12 乙醇相對氧氣水平的散點圖和針對發酵數據的糖類型指標直方圖

    若是你想用特定標籤來選擇數據點,你可使用name-list函數來生成帶標籤的窗體。這個窗體能夠被鏈接到任何圖形裏,選擇在命名列表裏的一個標籤,與之鏈接的那個圖形上的相對應的點就會高亮。你可使用帶一個數值參數的name-list函數,例如,(name-list 10)產生帶有標籤"0", ..., "9"的命名列表,你也能夠給出一個你本身的命名列表標籤。

練習 2.16

略。

2.5.5 修改散點圖

    在產生一個數據集的散點圖以後,你可能想要加一條線到圖形上,例如一條迴歸線。做爲例證,一項自行車道對車手的影響的研究,自行車手到馬路中間線的距離,自行車手與正經過的汽車的距離,記錄10個自行車手。數據輸入以下:

> (def travel-space
       '(12.8 12.9 12.9 13.6 14.5 14.6 15.1 17.5 19.5 20.8))
TRAVEL-SPACE
> (def separation
       '(5.5 6.2 6.3 7.0 7.8 8.3 7.1 10.0 10.8 11.0))
SEPARATION

伴隨隨機變量separation,能夠擬合這些數據的迴歸線有以下特徵:截距爲-2.18,斜率爲0.66。讓咱們看一下如何將這條線加到數據的散點圖的。

    咱們可使用以下表達式來產生這些數據點的散點圖:

> (plot-points travel-space separation)
#<Object: 1376d30, prototype = SCATTERPLOT-PROTO>

爲了能向圖形中加線,咱們必須可以在Lisp-Stat中引用到它。爲了完成這個任務,咱們能夠將plot-points函數返回的結果賦值給一個符號,即:

> (def myplot (plot-points travel-space separation))
MYPLOT

    plot-points函數返回的結果就是一個Lisp-Stat對象。爲了使用該對象,你須要向它發送一個消息,這要使用send函數來完成,其通常形式爲: (send <object> <message selector> <argument1> ... )

以下表達式將通知myplot添加一條a+bx的線,其中a=-2.18, b=0.66,。<message selector>是:abline這一關鍵字,數字-2.18和0.66是參數。消息由選擇器和參數組成。消息選擇器通常是一個Lisp關鍵字,也就是說它們是一個帶冒號的符號。結果見圖2.13。

圖2.13 帶擬合直線的自行車數據散點圖

散點圖對象還能理解一些其它消息,包括:help消息:

> (send myplot :help)
loading in help file information - this will take a minute ...done
SCATTERPLOT-PROTO
Scatterplot
Help is available on the following:

ABLINE ACTIVATE ADD-FUNCTION ADD-LINES ADD-METHOD ADD-MOUSE-MODE ADD-POINTS ADD-SLOT ADJUST-POINTS-IN-RECT ADJUST-SCREEN ADJUST-SCREEN-POINT ADJUST-TO-DATA ALL-POINTS-SHOWING-P ALL-POINTS-UNMASKED-P ALLOCATE ANY-POINTS-SELECTED-P APPLY-TRANSFORMATION BACK-COLOR BRUSH BUFFER-TO-SCREEN CANVAS-HEIGHT CANVAS-TO-REAL CANVAS-TO-SCALED CANVAS-WIDTH CENTER CHOOSE-MOUSE-MODE CLEAR CLEAR-LINES CLEAR-MASKS CLEAR-POINTS CLICK-RANGE CLIP-RECT CLOSE CONTENT-ORIGIN CONTENT-RECT CONTENT-VARIABLES COPY-TO-CLIP CURRENT-VARIABLES CURSOR CUT-TO-CLIP DELETE-DOCUMENTATION DELETE-METHOD DELETE-MOUSE-MODE DELETE-SLOT DISPOSE DO-CLICK DO-IDLE DO-KEY DO-MOTION DOC-TOPICS DOCUMENTATION DRAG-GREY-RECT DRAW-BITMAP DRAW-BRUSH DRAW-COLOR DRAW-DATA-LINES DRAW-DATA-POINTS DRAW-LINE DRAW-MODE DRAW-POINT DRAW-STRING DRAW-STRING-UP DRAW-SYMBOL DRAW-TEXT DRAW-TEXT-UP ERASE-ARC ERASE-BRUSH ERASE-OVAL ERASE-POLY ERASE-RECT ERASE-SELECTION ERASE-WINDOW FIND FIXED-ASPECT FOCUS-ON-SELECTION FRAME-ARC FRAME-LOCATION FRAME-OVAL FRAME-POLY FRAME-RECT FRAME-SIZE GET-METHOD H-SCROLL-INCS HAS-H-SCROLL HAS-METHOD HAS-SLOT HAS-V-SCROLL HELP HIDE-WINDOW IDLE-ON INTERNAL-DOC ISNEW LINE-TYPE LINE-WIDTH LINESTART-CANVAS-COORDINATE LINESTART-COLOR LINESTART-COORDINATE LINESTART-MASKED LINESTART-NEXT LINESTART-TRANSFORMED-COORDINATE LINESTART-TYPE LINESTART-WIDTH LINKED LOCATION MARGIN MASK-SELECTION MENU METHOD-SELECTORS MOUSE-MODE MOUSE-MODE-TITLE MOUSE-MODES MOVE-BRUSH NEW NUM-LINES NUM-POINTS NUM-VARIABLES OWN-METHODS OWN-SLOTS PAINT-ARC PAINT-OVAL PAINT-POLY PAINT-RECT PARENTS PASTE-FROM-CLIP PASTE-STREAM PASTE-STRING POINT-CANVAS-COORDINATE POINT-COLOR POINT-COORDINATE POINT-HILITED POINT-LABEL POINT-MASKED POINT-SELECTED POINT-SHOWING POINT-STATE POINT-SYMBOL POINT-TRANSFORMED-COORDINATE POINTS-HILITED POINTS-IN-RECT POINTS-SELECTED POINTS-SHOWING PRECEDENCE-LIST PRINT PROTO RANGE REAL-TO-CANVAS REDRAW REDRAW-CONTENT REMOVE REPARENT REPLACE-SYMBOL RESIZE RESIZE-BRUSH RETYPE REVERSE-COLORS ROTATE-2 SCALE SCALE-TO-RANGE SCALE-TYPE SCALED-RANGE SCALED-TO-CANVAS SCREEN-RANGE SCROLL SELECTION SELECTION-DIALOG SELECTION-STREAM SET-MODE-CURSOR SET-OPTIONS SET-SELECTION-COLOR SET-SELECTION-SYMBOL SHIFT SHOW SHOW-ALL-POINTS SHOW-WINDOW SHOWING-LABELS SIZE SLICE-VARIABLE SLOT-NAMES SLOT-VALUE START-BUFFERING TEXT-ASCENT TEXT-DESCENT TEXT-WIDTH TITLE TRANSFORMATION UNDO UNMASK-ALL-POINTS UNSELECT-ALL-POINTS UNSHOW-ALL-POINTS UPDATE USE-COLOR V-SCROLL-INCS VARIABLE-LABEL VIEW-RECT VISIBLE-RANGE WHILE-BUTTON-DOWN X-AXIS Y-AXIS 
NIL

對全部的散點圖來講,它們的主題列表都是類似的,可是對旋轉圖、散點圖矩陣和直方圖來講有一些不一樣。

    :clear消息,就像它的名字暗示的那樣,它會清除繪製的圖形並容許從頭開始建立一個新的圖形。另外兩個有用的消息是:add-points和:add-lines。爲了查明如何使用它們,咱們可使用:help消息並將:add-points和:add-lines做爲參數:

> (send myplot :help :add-points)
ADD-POINTS
Method args: (points &key point-labels (draw t))
or:          (x y  &key point-labels (draw t))
Adds points to plot. POINTS is a list of sequences, 
POINT-LABELS a list of strings. If DRAW is true the new points
are added to the screen. For a 2D plot POINTS can be replaced
by two sequences X and Y.
NIL
> (send myplot :help :add-lines)
ADD-LINES
Method args: (lines &key type (draw t))
or:          (x y  &key point-labels (draw t))
Adds lines to plot. LINES is a list of sequences, the 
coordinates of the line  starts. TYPE is normal or dashed. If
DRAW is true the new lines are added to the screen. For a 2D
plot LINES can be replaced by two sequences X and Y.
NIL

描述裏的術語sequence意思是一個list或者一個vector。

    圖2.13裏的統計圖展現了數據的曲率。在travel-sapce裏,separation一次項和二次項的迴歸計算將產生一些估計係數——常數估值爲-16.42,一次項係數估值爲2.433,二次項係數估值爲-0.05339。讓咱們使用:clear, :add-points和:add-lines消息改變myplot對象,用來顯示咱們擬合後的二次模型。首先咱們使用表達式定義x來表示從點12到點22之間的50等份的變量,定義y做爲這些值的擬合響應。

> (def x (rseq 12 22 50))
X
> (def y (+ -16.42 (* 2.433 x (* -0.05339 (* x x)))))
Y

圖2.14 帶擬合曲線的自行車數據散點圖

而後用下邊的表達式將圖形改爲圖2.14的樣子。

> (send myplot :clear)
NIL
> (send myplot :add-points travel-space separation)
NIL
> (send myplot :add-lines x y)
NIL

咱們已經使用plot-points來得到一個新的圖形,而後又用:add-lines修改了圖形,可是這裏咱們用到的方法容許咱們嘗試全部這3個消息。

2.5.6 動態模擬

做爲「經過修改現有統計圖形咱們能幹什麼」的有一個例證,讓咱們來構建一個動態模擬來檢測變量,這些變量來自於標準正態分佈樣本的直方圖,這個動態模擬就是一個小電影。爲了開始,使用以下表達式構建一個單變量的直方圖而後保存它的繪圖對象爲h。

> (def h (histogram (normal-rand 20)))
H

:clear消息對直方圖也是可用的,就像咱們在它的幫助信息裏看到的那樣:

> (send h :help :clear)
CLEAR
Message args: (&key (draw t))
Clears the plot data. If DRAW is nil the plot is redrawn; otherwise its
current screen image remains unchanged.

它帶一個可選的關鍵字函數。若是這個參數是nil,那麼統計圖不會重畫直到一些其它消息致使它重畫,這對動態模擬式很重要的。倒換和從新調整窗體大小,直到當你輸入命令到解釋器時仍能看到直方圖窗口爲止。(注:其實不用那麼麻煩,在菜單欄的windows->Tile直接將全部窗口平鋪就好了)而後輸入表達式:

> (dotimes (i 50)
           (send h :clear :draw nil)
           (send h :add-points (normal-rand 20)))

這應該能產生一個50個直方圖的序列,每個都基於一個大小爲20的正態分佈樣本。經過輸入關鍵字參數:draw並賦值nil給:clear消息,你就能保證每一個直方圖都在屏幕上待一下子直到下一個直方圖準備繪製的時候。再次運行這個例子,不帶這個參數,看看有什麼不一樣。

    編程實現動態模擬提供了另外一個查看一些變量之間關係的方法。做爲一個簡單的例子,假設咱們想檢測標量abrasion-loss和hardness之間的關係,這些變量來自第2.5.2節。咱們用下邊的表達式生成一個abrasion-loss變量的直方圖:

> (def h (histogram abrasion-loss))
H

對於動態模擬來講,:point-selected, :point-hilited和:point-showing消息是頗有用的。如下是直方圖裏對:point-selected消息的信息:

> (send h :help :point-selected)
POINT-SELECTED
Method args: (point &optional selected)
Sets or returns selection status (true or NIL) of POINT. Sends 
:ADJUST-SCREEN message if states are set. Vectorized.

如此你可使用這個消息檢測一個點當前是否是選中狀態,也能夠選中或者取消選中。再次倒換窗口到輸入命令時可以看到直方圖,而後鍵入表達式:

> (dolist (i (order hardness))
          (send h :point-selected i t)
          (send h :point-selected i nil))

表達式(order hardness)將產生一個hardness的有序值列表,從而返回結果的第一個元素是hardness的第一個元素,第二個元素是hardness的第二小的元素,以此類推。對相應數據的循環移動進行選擇與取消選擇操做。

    屏幕上的結果與從左到右刷hardness變量直方圖的結果是類似的,hardness直方圖是與abrasion-loss變量的直方圖相鏈接的。這種方法的缺陷就是與用鼠標實現相比,寫出表達式很難。換句話說,當你用鼠標刷圖形數據的時候,你趨向於將經歷集中在你刷的圖形上而不是另外一個鏈接圖上。當你運行一個動態模擬圖時,當模擬程序運行的時候你什麼也不用作,所以能夠將全部精力集中在結果上。

    一箇中間方案是可行的:你能夠設置一個帶滾動條的對話框,滾動條的值就是(order hardness)列表的全部值,當它滾動的時候就會選擇相應的點,2.7節將展現這種方案如何工做的例子。

    在不少Lisp系統裏,就像咱們如今描述的這個模擬器同樣,偶爾會停一下子,緣由就是須要垃圾回收機制來動態地回收已經分配的存儲單元。雖然說是這樣,在一些簡單的模擬器裏,好比這些迭代器對你來講,可能仍然太快了而不能捕捉任何處理的形式。爲了下降速度,你能夠用函數pause來加入一個短的延時,經過使用一個整數n,pause將致使可執行程序的掛起,大約要n/60秒鐘,例如,你可使用下式代替前邊的表達式:

> (dolist (i (order hardness))
          (send h :point-selected i t)
          (pause 10)
          (send h :point-selected i nil))

這將確保循環將以大概每秒6個的速度穿過數據點。

練習 2.17

略。

練習 2.18

略。

2.6 迴歸

    迴歸模型是使用Lisp-Stat的對象和消息傳遞機制實現的。這在2.5.5節已經介紹了。再繼續閱讀以前,你可能須要簡單複習一下那節內容。

    讓咱們來對2.5.5節的自行車數據作一個簡單的擬合。因變量是separation,自變量是tranvel-space。爲了造成迴歸模型,可使用regression-mode函數:

> (regression-model travel-space separation)

Least Squares Estimates:

Constant                  -2.18247      (1.05669)
Variable 0                0.660342      (6.747931E-2)

R Squared:                0.922901    
Sigma hat:                0.582108    
Number of cases:                10
Degrees of freedom:              8

#<Object: 141cbbc, prototype = REGRESSION-MODEL-PROTO>

    regression-model函數的基本語法是(regression-model <x> <y>)。對於一個簡單的迴歸模型,<x>能夠是一個列表的列表、矢量的列表或者一個矩陣。regression-model函數也包含一些可選參數,包括:intercept, :print和:weights,:intercept和:print的默認參數都是t,想要得到一個沒有截距的模型,使用如下表達式:

> (regression-model x y :intercept nil)

爲了造成一個帶權重的迴歸模型,使用如下表達式:

> (regression-model x y :weights w)

這裏的w是一個與y長度相等的權重列表或權重矢量,假設偏差方差與權重w成反比。

    regression-model函數打印擬合模型的簡單概述,而且返回一個模型對象。爲了可以進一步檢測模型,用如下表達式將返回的模型對象賦給一個變量:

> (def bikes (regression-model travel-space separation :print nil))
BIKES

值爲nil的關鍵字參數:print壓縮了咱們看到的模型的概述。爲了找出那些消息時可用的,咱們使用:help消息:

> (send bikes :help)
loading in help file information - this will take a minute ...done
REGRESSION-MODEL-PROTO
Normal Linear Regression Model
Help is available on the following:

ADD-METHOD ADD-SLOT BASIS CASE-LABELS COEF-ESTIMATES COEF-STANDARD-ERRORS COMPUTE COOKS-DISTANCES DELETE-DOCUMENTATION DELETE-METHOD DELETE-SLOT DF DISPLAY DOC-TOPICS DOCUMENTATION EXTERNALLY-STUDENTIZED-RESIDUALS FIT-VALUES GET-METHOD HAS-METHOD HAS-SLOT HELP INCLUDED INTERCEPT INTERNAL-DOC ISNEW LEVERAGES METHOD-SELECTORS NEW NUM-CASES NUM-COEFS NUM-INCLUDED OWN-METHODS OWN-SLOTS PARENTS PLOT-BAYES-RESIDUALS PLOT-RESIDUALS PRECEDENCE-LIST PREDICTOR-NAMES PRINT PROTO R-SQUARED RAW-RESIDUALS REPARENT RESIDUAL-SUM-OF-SQUARES RESIDUALS RESPONSE-NAME RETYPE SAVE SHOW SIGMA-HAT SLOT-NAMES SLOT-VALUE STUDENTIZED-RESIDUALS SUM-OF-SQUARES SWEEP-MATRIX TOTAL-SUM-OF-SQUARES WEIGHTS X X-MATRIX XTXINV Y 
NIL

這些消息裏不少都是自解釋的,不少已經被:display消息使用了,regression-model函數將:display消息發送給新的模型並打印其概述信息。做爲一個例子,讓咱們試一下:coef-estimates和:coef-standard-errors消息:

> (send bikes :coef-estimates)
(-2.182471511502903 0.6603418619651689)
> (send bikes :coef-standard-errors)
(1.0566881307586837 0.06747931359188968)
>

按理說,若是模型裏有截距的話,由這些消息返回的列表裏的實體是與截距對應的,當提供給regression-model函數的時候,該截距後邊跟着因變量。然而,若是在計算過程當中探測到了退化現象,在擬合過程當中一些變量就不會被使用,在打印的概述信息裏它們將被用別名的方式標記出來。已經使用的變量的標度將由:basis消息獲取,若是合適的話,由:coef-estimates返回的列表項將與截距對應,與截距一塊兒的還有各元素的係數。消息:x-matrix和:xtxinv是類似的。

    :plor-residuals消息產生一個殘留的統計圖,爲了查明殘留的統計圖針對什麼,讓咱們看一下它的幫助信息吧:

> (send bikes :help :plot-residuals)
loading in help file information - this will take a minute ...done
PLOT-RESIDUALS
Message args: (&optional x-values)
Opens a window with a plot of the residuals. If X-VALUES are not supplied 
the fitted values are used. The plot can be linked to other plots with the 
link-views function. Returns a plot object.
NIL

使用如下兩個表達式,咱們能夠構架圖2.15所示的兩個數據統計圖:

> (plot-points travel-space separation)
#<Object: 1350994, prototype = SCATTERPLOT-PROTO>
> (send bikes :plot-residuals travel-space)
#<Object: 1357680, prototype = SCATTERPLOT-PROTO>

圖2.15 針對自行車實例鏈接原始數據和殘存圖的示意圖

經過鏈接這些統計圖(注:使用菜單欄裏的plot->link view菜單選項),咱們可使用鼠標在兩個統計圖裏同步識別數據點。圖中一個標號爲6的點脫穎而出。

    這兩幅統計圖代表數據是存在曲率的,若是咱們忽略觀測點6,殘留圖的曲率特別明顯。爲了容許這個曲率的出現,咱們能夠試圖用二階項對travel-space變量來擬合出一個模型。

> (def bikes2
       (regression-model (list travel-space (^ travel-space 2))
                         separation))

Least Squares Estimates:

Constant                  -16.4192      (7.84827)
Variable 0                 2.43267      (0.971963)
Variable 1                -5.339121E-2  (2.922567E-2)

R Squared:                0.947792    
Sigma hat:                0.512086    
Number of cases:                10
Degrees of freedom:              7

BIKES2

我使用冪指數函數「^」來計算travel-space變量的平方,即咱們造成了一個多元迴歸模型,該回歸模型的第一個參數是一個自變量的列表。

    從這個點開始你能夠向多個方向繼續前進,若是你想計算觀測點的庫克距離(注:統計學專有名詞,詳見鏈接),你能夠向bikes2對象發送:cooks-distance消息:

> (send bikes2 :cooks-distances)
(0.16667299257295504 0.009185959754912747 0.030268009774083452 0.011098970664932503 0.009584417522457393 0.12066535826085792 0.5819289720128115 0.04601789946507566 0.006404473564781199 0.09400811436381019)

爲了檢測這個距離,你能夠首先計算內部學生化殘差,

> (def studres (/ (send bikes2 :residuals)
                  (* (send bikes2 :sigma-hat)
                     (sqrt (- 1 (send bikes2 :leverages))))))
STUDRES

而後再獲取距離的方式,

> (* (^ studres 2)
     (/ (send bikes2 :leverages)
        (- 1 (send bikes2 :leverages))
        3))
(0.16667299257295504 0.009185959754912747 0.030268009774083452 0.011098970664932505 0.009584417522457393 0.12066535826085789 0.5819289720128115 0.046017899465075666 0.006404473564781199 0.09400811436381019)

來直接計算它們。第7個數據元素——下標從0開始的觀測點6,很明顯地脫穎而出。

    檢測殘留圖中可能的異常值的另外一個方法是就是使用貝葉斯殘留圖,該方法是由Chaloner和Brant倡導的,可使用:plot-bayes-residuals消息得到該功能。如下表達式將產生圖2.16所示的統計圖。

> (send bikes2 :plot-bayes-residuals)
#<Object: 131b5a0, prototype = SCATTERPLOT-PROTO>

圖2.16 自行車數據的貝葉斯殘留圖

圖2.16中條形顯示的意思是實際實現偏差的後驗分佈的±2SD,它基於對迴歸係數的一個不當的均勻先驗分佈。(注:這句話可能不太準確,請參照原文),y軸是以σ爲單位的。從而本圖展現了一個蓋然率,數據點6偏離平均值爲3或更大標準方差的蓋然率大約是3%;偏離平均值至少2個單位標準協方差的蓋然率在50%上下。

練習 2.19

略。

練習 2.20

略。

2.7 定義函數和方法

    本節將對Lisp-Stat的一個簡潔的介紹。最基礎的編程操做是定義函數。與其緊密相關的是爲對象定義一個方法。

2.7.1 定義函數

    用來定義函數的特殊形式叫defun,最簡單的defun語法形式是 (defun <name> <parameters> <body>),這裏的<name>是一個符號,用來做爲函數名,<parameters>是符號的一個列表,命名爲參數,<body>表示一個或多個表達式,組成了函數體。例如,假設你想要定義一個函數從一個列表裏刪除一個元素,函數應該將列表和你想刪除的元素的下標做爲參數。函數體能夠基於上邊2.4.3節描述的兩個方法中的一種。這裏是其中的一個:

> (defun delete-case (x i)
    (select x (remove i (iseq (length x)))))
DELETE-CASE

defun的參數都沒有被引用(即沒有在符號前加單引號):defun是一個特殊形式,不須要對參數求值。

    你能夠寫一個向對象傳遞消息的函數。這裏有一個函數,它帶兩個迴歸模型變量參數,假設這些模型是被嵌套的,爲了比較這些模型來計算F-統計:

> (defun f-statistic (m1 m2)
	   "Args: (m1 m2)
           Computes the F statistic for testing model m1 within model m2."
	   (let ((ss1 (send m1 :sum-of-squares))
		 (df1 (send m1 :df))
		 (ss2 (send m2 :sum-of-squares))
		 (df2 (send m2 :df)))
	     (/ (/ (- ss1 ss2) (- df1 df2)) (/ ss2 df2))))
F-STATISTIC

這個函數定義使用Lisp的let構造來創建一些局部變量。變量ss1, df1, ss2, df2是爲模型參數m1和m2定義的,以後用來計算F統計。參數列表後的字符串叫文檔字符串。當defun表達式出現文檔字符串的時候,defun將它安裝在函數內部供help函數獲取。

2.7.2 函數做爲參數

    在第2.2.3節中,咱們使用plot-funcion來繪製[-pi, pi]範圍內的正弦函數。你可使用相同的方法來繪製你本身的函數。例如,假設你想繪製-2到3範圍內的f(x)=2x+x^2的圖形,你能夠定義一個函數f:

> (defun f (x) (+ (* 2 x) (^ x 2)))
F

而後使用plot-function:

> (plot-function #'f -2 3)

回憶一下,表達式#'f是(function f)的縮寫,用來得到與符號f關聯的那個Lisp函數。

    使用spin-function你能夠構建一個有兩個變量的函數的旋轉圖形,例如,定義f爲:

(defun f (x y) (+ (^ x 2) (^ y 2)))
F

那麼如下表達式將構建一個自變量x在[-1, 1]內,自變量y在[-1, 1]內的函數f(x, y)=x^2+y^2的函數圖像,該圖像用6X6網格表示,網格點的數量可使用:num-points關鍵字改變。結果如圖2.17所示。

    contour-function函數的必選參數數量與spin-function函數相同,將產生一個函數參數的輪廓線圖。除了:num-points關鍵字,你也能夠賦予contour-function函數:levels關鍵字,緊跟着一個輪廓等級的列表。

    本例中,不得不定義一個函數僅僅是做爲另外一個函數的參數,這種方法有些笨拙。可以寫一個產生咱們所需功能的表達式,而後將這個表達式做爲plot-function和spin-function的參數,這樣可能更方便一些。第3.6.1節將展現如何完成這個任務。

練習 2.21

2.73 繪圖動畫控制

做爲函數使用的另外一種展現,假設在Box-Cox能量變換中Box-Cox變換是統計建模中經常使用的一種數據變換,用於連續的響應變量不知足正態分佈的狀況),咱們想檢測改變指數對正態機率圖形形成的影響。

第一步,定義函數以計算冪轉換,並正規化轉換結果,使它落在0和1之間:

> (defun bc (x p)
    (let* ((bcx (if (< (abs p) .0001) (log x) (/ (^ x p) p)))
           (min (min bcx))
           (max (max bcx)))
      (/ (- bcx min) (- max min))))
BC

這個定義裏使用了let*形式來創建一個局部變量的綁定。let*形式與上邊用到的let形式類似,但在一下幾點是不一樣的,let*順序地定義變量,容許使用在let*中定義的其它變量來定義一個新的變量,相反地,let並行地建立其分配符。在這種狀況下,變量min和max根據bcx來定義。

    而後,咱們須要一個正數集合來作變換,讓咱們來使用2.2.1節裏介紹的降水量數據樣本,咱們須要對觀測值進行排序:

> (def x (sort-data precipitation))
X

指望的統一順序統計的常態分佈分位數由下式給出:

> (def nq (normal-quant (/ (iseq 1 30) 31)))
NQ

未經轉換的數據的機率統計圖由下式構建:

> (def myplot (plot-points nq (bc x 1)))
MYPLOT

由於這裏的冪是1,因此函數bc僅僅是對數據進行了尺寸調整。

    爲了改變myplot對象的轉換中冪的值,咱們能夠定義一個change-power函數:

> (defun change-power (p)
    (send myplot :clear :draw nil)
    (send myplot :add-points nq (bc x p)))
CHANGE-POWER

想這樣對一個新表達式求值:

> (change-power .5)

而後對新的冪值,實例中的平方根轉換,重畫統計圖,

    咱們可使用不一樣的冪值重複這個處理過程,目的是獲取轉換的影響的感受,不是這個方法是很笨拙的。一個可代替的方法就是設置一個滑動對話框來控制這個冪的值。這個滑動的對話框是一個無模式的對話框,包含一個滑塊和一個顯示的值域。當滑塊移動的時候,顯示值改變,同時動做觸發。這裏的動做由動做函數定義,這個函數使用當前滑塊的值做爲參數進行調用,每次這個值由用戶改變,咱們能夠用以下表達式構建針對咱們的問題的一個合適滑塊:

> (sequence-slider-dialog (rseq -1 2 31) :action #'change-power)

這個滑塊能夠經過如下值進行移動, -1.0, -0.8, ..., 1.8, 2.0,每次調用的時候都使用當前的移動滑塊獲得的冪值。圖形和滑塊見圖2.18所示。

圖2.18 一個滑塊控制冪轉換的統計圖

2.7.4 定義方法

    當消息發送給對象的時候,對象系統對象和消息選擇器來查找最合適的代碼段執行。從而,不一樣的對象對相同的消息做出不一樣的響應,線性迴歸模型和非線性迴歸模型都會對:coef-estimates消息進行響應,可是他們會執行不一樣的代碼來計算該響應。

    對象使用的用來響應消息的代碼叫作方法。對象被組織在層級關係裏,在那裏對象從它們祖先那裏繼承其它對象。若是一個對象沒有它本身的用來響應消息的方法,它就會使用從祖先那裏繼承的方法。send函數將提高該對象的祖先的列表的優先級,直到找到針對消息的方法。

    大多數咱們目前接觸到的對象都是直接繼承自原型對象。Scatterplots(散點圖)繼承自scatterplot-proto,histograms(直方圖)繼承自histogram-proto,迴歸模型對象繼承自regression-model-proto。原型對象與其它任何對象都是類似的,它們是定義必定類別的對象的典型版本,這類對象定義了默認的行爲。可是任何對象都有一個方法,在調試新方法的過程當中,最好將這個方法鏈接到獨立構建的對象上,達到代替原型的目的。

    defmeth宏用來向對象添加方法。方法定義的通常形式是:(def <object> <selector> <parameters> <body>),<object>是一個將擁有新方法的對象,參數<selector>是對方法的消息選擇關鍵字,參數<parameters>是方法的參數名列表,<body>是組成方法體的一個或多個表達式。當方法被使用時,這些表達式將按順序求值。

    做爲一個簡單的展現,在第2.5.6節咱們使用一個循環運行一個模擬器,在那裏咱們在一幅直方圖裏展現了大小爲20的樣本序列,該序列來源於一個正態分佈。原始的直方圖由下式構建:

> (def h (histogram (normal-rand 20)))
H

每次迭代的時候,咱們首先清除直方圖,而後增長一個新的樣本,經過爲新消息定義一個方法,咱們能夠簡化一下咱們的直方圖:

> (defmeth h :new-sample ()
    (send self :clear :draw nil)
    (send self :add-points (normal-rand 20)))
:NEW-SAMPLE

運行50次動畫的循環如今能夠寫成這樣:

> (dotimes (i 50) (send h :new-sample))

    在爲:new-sample消息的方法定義時須要一點解釋,self變量的實用。對象的方法常常須要可以參考獲取消息的對象,由於方法能夠被繼承,在寫方法的時候,你可能對接受的對象的標識不是很肯定。所以對象系統使用了一個慣例,當表達式在被求值的消息對應的那個方法的裏時,接收對象就是變量self的值。

    如今咱們有一個簡單的消息,用來向咱們的直方圖發送消息,告訴他顯示一條新數據樣本。擁有一個發送消息的簡單的辦法,而不須要向解釋器敲入命令,這樣好極了。直方圖菜單提供了一個方便的工具來達到該目的,表達式是:

> (let ((m (send h :menu))
        (i (send menu-item-proto :new "New Sample"
                 :action #'(lambda () (send h :new-sample)))))
    (send m :append-item i))

向直方圖的菜單末端加入一個新的菜單項,每次選中該菜單項的時候都將改變顯示的數據樣本。在閱讀過第6和第7章以後,你將可以理解這個表達式。

2.8 更多模型和技術

    Lisp-Stat爲非線性迴歸模型、最大似然估計和近似貝葉斯推理提供了一些工具。這三個工具集的共同特色是描述一個特定問題所須要的部分信息是函數提供的,例如,一個爲非線性迴歸設立的函數。本節假設你熟悉非線性迴歸基礎,最大似然估計和貝葉斯推理:

2.8.1 非線性迴歸

    Lisp-Stat容許非線性的正態迴歸模型。舉個例子,Bates和Watts描述了一個變量y和變量x之間關係的試驗,變量y表示酶促反應的速度,變量x表示底物濃度。試驗中的濃度和觀測速度由下式給定,試驗中使用嘌呤黴素對酶促反應進行處理:

> (def x1 (list .02 .02 .06 .06 .11 .11 .22 .22 .56 .56 1.1 1.1))
X1
> (def y1 (list 76 47 97 107 123 139 159 152 191 201 207 200))
Y1
>

對於速度對底物濃度的依賴,米氏函數一般能夠提供一個好的模型。

假定在一個給定的濃度水平上,米氏函數表明平均速度,其函數f1定義爲:

 

> (defun f1 (theta)
    (/ (* (select theta 0) x1) (+ (select theta 1) x1)))
F1

f1函數將對參數列表theta的平均反應時間的列表進行計算,這是在 數據點x1處作的計算。

    使用這些定義,咱們可使用nreg-model函數構建一個非線性迴歸模型。首先,咱們須要對兩個模型的參數進行初始化估計,爲米氏模型檢測表達式這一操做顯示,當x增加時,米氏函數逼近一條漸近線,θ0。第二個參數θ1能夠解釋成當函數達到漸近線一半時x的值。利用這些對參數的解釋,下式將構造一個統計圖,圖示見圖2.19,咱們能夠對θ0讀出200,對θ1讀出0.1這些合理的初始化估計。

圖2.19 酶促反應試驗的底物濃度對反應速度的統計圖

nreg-model函數須要的參數有平均值函數、響應矢量和參數的初始估計列表。nreg-model函數經過使用一個交互式的算法來計算精確的估計值,這個交互式的算法以初始猜想凱斯,打印結構的概要,並返回一個非線性迴歸模型對象:

 

> (def puromycin (nreg-model #'f1 y1 (list 200 .1)))
Residual sum of squares: 7964.19
Residual sum of squares: 1593.16
Residual sum of squares: 1201.03
Residual sum of squares: 1195.51
Residual sum of squares: 1195.45
Residual sum of squares: 1195.45
Residual sum of squares: 1195.45
Residual sum of squares: 1195.45

Least Squares Estimates:

Parameter 0                212.684      (6.94715)
Parameter 1                6.412127E-2  (8.280941E-3)

R Squared:                0.961261    
Sigma hat:                 10.9337    
Number of cases:                12
Degrees of freedom:             10

PUROMYCIN

nreg-model函數也帶一些關鍵字參數,包括用:epsilon來制定一個收斂準則,用:count-limit作爲迭代次數的極限值,這兩個值得默認值分別爲0.0001和20。使用值爲nil的:verbose關鍵字來抑制每次迭代過程的的平方和殘差信息,使用值爲nil的:print關鍵字將約束概述信息,擬合模型的算法是一個帶回溯簡單高斯-牛頓算法,導數是經過數值計算的。

    爲了使你明瞭如何進一步分析模型,你能夠向puromycin對象發送:help消息。結果與線性迴歸模型的幫助信息類似。緣由很簡單:非線性迴歸模型是用對象實現的,而非線性迴歸模型的原型nreg-model-proro是繼承自線性迴歸模型的原型的。內部數據、計算估計值的方法和計算擬合值得方法都作了修改以適應非線性模型,可是其它大多數方法都沒有改變。一旦模型完成擬合,預估參數值得平均函數的雅克比矩陣將被定義爲X矩陣來使用。在此定義的術語中,大多數針對線性迴歸的方法,像:coef-standard-errors和:leverages等,仍然是有意義的,至少是一階逼近的。

    除了能夠相應線性迴歸模型的消息以外,非線性迴歸模型還能相應如下消息:

  •     :COUNT-LIMIT
  •     :EPSILON
  •     :MEAN-FUNCTION
  •     :NEW-INITIAL-GUESS
  •     :PARAMETER-NAMES
  •     :THETA-HAT
  •     :VERBOSE

練習 2.23

略。

練習 2.24

略。

2.8.2 最大和最小似然估計

    Lisp-Stat有兩個函數能夠求得一些變量的最大似然率。第一個是newtonmax函數,它帶一個函數和一個列表或者矢量值來表示初始猜想值,該初始猜想針對最大值的位置,函數試圖用基於帶回溯的牛頓法來找到最大值。該算法基於Dennis和Schnabel描述的無約束的最小系統。

    舉個例子,Proschan描述了定時收集的一些數據,這些數據源自幾個飛機場空調單元拓機的時間間隔。其中一個機場的數據輸入以下:

> (def x (list 90 10 60 186 61 49 14 24 56 20 79 84 44 59 29 118
               25 156 310 76 26 44 23 62 130 208 70 101 208))
X

這些數據的一個簡單的模型做出以下假設:拓機時間間隔是獨立隨機變量,知足通用gamma分佈,gamma分佈的密度能夠寫成這樣:

這裏的μ是拓機時間的平均值,β是gamma的形狀參數,那麼該數據樣本的log似然函數應該這樣給出:

咱們能夠定義一個Lisp函數來計算這個log似然函數:

> (defun gllik (theta)
    (let* ((mu (select theta 0))
           (beta (select theta 1))
           (n (length x))
           (bym (* x (/ beta mu))))
      (+ (* n (- (log beta) (log mu) (log-gamma beta)))
         (sum (* (- beta 1) (log bym)))
         (sum (- bym)))))
GLLIK

該定義使用函數log-gamma對log(Γ(β))。

    封閉形式的最大似然機率估計最這個分佈的形狀參數是不可用的,可是咱們可使用newtonmax函數來用數值方法進行計算。咱們須要一個廚師猜想值用來做爲最大似然的起始值。爲了得到初始估計,咱們能夠計算x的指望和標準方差:

> (mean x)
83.51724137931036
> (standard-deviation x)
70.80588186304823

而後使用矩估計方法估計 ,像這樣計算:

> (^ (/ (mean x) (standard-deviation x)) 2)
1.3912770125276501

使用這些起始值,咱們能夠最大化log似然函數:

> (newtonmax #'gllik (list 83.5 1.4))
maximizing...
Iteration 0.
Criterion value = -155.603
Iteration 1.
Criterion value = -155.354
Iteration 2.
Criterion value = -155.347
Iteration 3.
Criterion value = -155.347
Reason for termination: gradient size is less than gradient tolerance.
(83.51726798264603 1.670986791540375)

一些狀態信息打印出來用於繼續優化。經過使用值爲nil的:verbose關鍵字,你能夠關閉這些信息。

    你可能想要檢查這個函數的確接近0的梯度。若是對這個梯度沒有一個封閉的表達式的話,你能夠用numgrad函數計算一個用數值上的近似值。對於咱們的例子:

> (numgrad #'gllik (list 83.5173 1.67099))
(-4.072480396396883E-7 -1.2576718377326971E-5)

梯度的元素確實至關接近0.

    你也能夠計算二階導數,或者使用numhess函數計算Hessian矩陣。最大似然估計的近似標準差能夠這樣給定,先對最大值處的Hessian矩陣取負,再求逆矩陣,再去該逆矩陣的對角線元素,最後對其求平方根:

> (sqrt
   (diagonal
    (inverse (- (numhess #'gllik (list 83.5173 1.67099))))))
(11.997551713830456 0.40264778727213846)

    使用newtonmax代替求取最大值,而後分別計算二階導數。經過使用值爲t的關鍵字參數:return-derivs,你可讓newtonmax返回一個元素位置列表,該列表的元素包括最大值、最有函數值、梯度和Hessian矩陣。

    牛頓法假設函數是兩次可微的,若是你的函數不平滑,或者若是你由於其它緣由對使用newtonmax存在麻煩。你可能要嘗試第二個最大化函數,nelmeadmax函數,該函數接收一個函數和一個初始化猜想值,而後試圖用Nelder-Mead單純形法找打最大值,該方法是由Press、Flannery、Teukolsky和Vetterling描述的。初始猜想值由一個單點組成,它以n個數字的列表或矢量的形式展示。初始猜想值也多是一個單形,即一個n維問題的n+1個點的列表。若是你制定了一個單點,那你也應該使用:size關鍵字參數去指定一個列表或者矢量,其長度爲n,且須要在初始單形的每個維度裏指定。這應該展現初始範圍的大小,在那個範圍裏算法將開始它的搜索直到找到最大值,咱們能夠在咱們都額gamma例子裏使用這個方法:

 

> (nelmeadmax #'gllik (list 83.5 1.4) :size (list 5 0.2))
Value = -155.52959743506827
Value = -155.52959743506827
Value = -155.52959743506827
Value = -155.5294988579447
Value = -155.46938588692785
Value = -155.46938588692785
Value = -155.42198557523471
Value = -155.42198557523471
Value = -155.40996834939776
Value = -155.3866822915446
Value = -155.38103392707978
Value = -155.36281930571423
Value = -155.36281930571423
Value = -155.34979141245745
Value = -155.34979141245745
Value = -155.34935123123307
Value = -155.3470397238531
Value = -155.3470397238531
Value = -155.3469673747863
Value = -155.34689776692113
Value = -155.3468231506837
Value = -155.3468231506837
Value = -155.34680431458563
Value = -155.34680431458563
Value = -155.3467951205185
Value = -155.34679314554154
Value = -155.3467903355422
Value = -155.3467903355422
Value = -155.346789572537
Value = -155.3467892236572
Value = -155.3467892236572
Value = -155.3467892236572
Value = -155.3467892236572
(83.5156874768436 1.6711172274500135)

這個方法與牛頓法相比有些長,可是它確實得到了相同的結果。

練習 2.25

略。

2.8.3 近似貝葉斯計算

    本小節描述的函數能夠用來計算後驗矩的一階和二階近似值,還有一階邊際後驗密度的鞍點狀近似。基於文獻[62,63,64],近似值假定後驗密度是光滑的,而且以單模式爲主。

    讓咱們從一個簡單的例子開始,一份關於白血病患者存活時間和病人體內白細胞數量的關係的研究數據集,該數據以周爲單位。該數據由兩類病人組成,AG+和AG-型。17個AG+型病人的數據能夠輸入以下:

> (def wbc-pos (list 2300 750 4300 2600 6000 40500 10000 17000
                     5400 7000 9400 32000 35000 100000 100000
                     52000 100000))
WBC-POS
> (def times-pos (list 65 156 100 134 16 108 121 4 39 143 56
                       26 22 1 1 5 65))
TIMES-POS

高白細胞數量表示處在疾病的比較嚴重的階段,也就是存活的機會更低。

    針對這些數據會用到一個模型:假設存活時間是按指數分佈的,它帶一個平均值參數,該平均值參數使用log-linear方法將"白細胞數量的對數"這一非線性指標進行"對數變換"以獲得線性模型,而後它纔可用於線性迴歸模型之中。爲了方便,我將對白細胞數量數據擴大10,000倍。也就是說,一個病人白細胞數量WBCi對應的平均存活時間爲

這裏的xi=log(WBCi/10000)。log似然函數給定爲:

這裏的yi表明存活時間。計算轉換後的WBC變量爲:

> (def transformed-wbc-pos (- (log wbc-pos) (log 10000)))
TRANSFORMED-WBC-POS

以後能夠用如下函數計算log似然率:

> (defun llik-pos (theta)
    (let* ((x transformed-wbc-pos)
           (y times-pos)
           (theta0 (select theta 0))
           (theta1 (select theta 1))
           (t1x (* theta1 x)))
      (- (sum t1x)
         (* (length x) (log theta0))
         (/ (sum (* y (exp t1x)))
            theta0))))
LLIK-POS

    我將使用一個模糊的、不肯定的先驗分佈來看待這個問題,該分佈在θi>0範圍內十個常量;從而,log後驗密度與上邊構造的log似然函數式等價的,最多多出一個常數。第一步就是用bayes-model函數構造一個貝葉斯模型對象。該函數能夠接受一個計算log後驗分佈密度的函數和這個後驗模型的初始猜想數據,該函數經過這個初始猜想數據的交互方法計算後驗模型,並打印後驗分佈信息的簡短概述,而後返回一個模型對象。咱們可使用函數llik-pos來計算log後驗密度,因此全部咱們能作的就是對後驗模型的初始值估計。由於咱們使用的模型假設是一個線性關係,是表示平均存活時間的對數與轉換後的WBC變量的線性關係。,對transformed-wbc-pos變量的存活時間的對數的一個線性迴歸模型應該提供一個合理的估計,該線性迴歸模型這樣給出:

> (regression-model transformed-wbc-pos (log times-pos))

Least Squares Estimates:

Constant                   3.58281      (0.328336)
Variable 0               -0.736302      (0.225778)

R Squared:                0.414869    
Sigma hat:                 1.32459    
Number of cases:                17
Degrees of freedom:             15

#<Object: 13daa98, prototype = REGRESSION-MODEL-PROTO>

因此說模型的合理初始猜想評估值是θ0評估值=exp(3.5),θ1的評估值=0.8,如今咱們能夠將這些評估值放到bayes-model函數裏:

> (def lk (bayes-model #'llik-pos (list (exp 3.5) 0.8)))
maximizing...
Iteration 0.
Criterion value = -96.3803
Iteration 1.
Criterion value = -87.4612
Iteration 2.
Criterion value = -85.2634
Iteration 3.
Criterion value = -84.8285
Iteration 4.
Criterion value = -84.7883
Iteration 5.
Criterion value = -84.7877
Iteration 6.
Criterion value = -84.7877
Reason for termination: gradient size is less than gradient tolerance.


First Order Approximations to Posterior Moments:

Parameter 0                60.6309      (15.0481)
Parameter 1               0.390576      (0.175448)

LK

經過使用值爲nil的:print關鍵字參數,來壓縮概述信息是可能的。可使用:verbose關鍵字壓縮迭代信息。

    bayes-model函數打印的概述信息給出了參數的後驗均值和標準方差的一階近似值。也就是說均值是由後驗模型的元素近似獲得的,標準方差是由模型的log後驗值得負Hessian矩陣的逆矩陣的對角線元素的平方根獲得的。這些近似值也能夠經過向模型發送:lstmoments消息來得到:

> (send lk :1stmoments)
((60.63089554959938 0.3905764020242036) (15.048125818750082 0.1754475290888174))

結果是兩個列表值的列表結構,它們是平均值列表和標準方差列表。

    一個相對較慢可是通常狀況下更精確的二階近似值也是可用的。它可使用:moments得到:

> (send lk :moments)
((69.89109467010141 0.4033868081199018) (18.577055399066786 0.1829812572509589))
> (send lk :moments (list 0 1))
((69.89109467010141 0.4033868081199018) (18.577055399066786 0.1829812572509589))

    θ0的一階近似值和二階近似值有一點不一樣,特別是,看起來平均值比模型的大些。這代表參數的後驗分佈是向右傾斜的,經過查看近似邊際後驗密度圖咱們能夠確認這一點。消息:marginl帶一個指標參數,還有一個進行密度求值的點數據序列,最後以列表形式返回它的值,該列表的元素就是使用的序列和在這些數據點處的近似密度值。爲了產生一個邊際密度圖,該列表可使用plot-lines函數給出:

> (plot-lines (send lk :margin1 0 (rseq 30 120 30)))
#<Object: 1401734, prototype = SCATTERPLOT-PROTO>

結果見圖2.20所示,確實看到圖形有些向右傾斜。

圖2.20 θ0的近似邊際後驗分佈

    除了檢測單獨參數,也能夠查看參數的光滑函數的後驗分佈。例如,你可能想問一個WBC=50000的病人存活一年以上的機率,那麼該方法包含了什麼信息,機率以下:

其中

x=log(5),能夠由如下函數計算:

> (defun lk-sprob (theta)
    (let* ((time 52.0)
           (x (log 5))
           (mu (* (select theta 0)
                  (exp (- (* (select theta 1) x))))))
      (exp (- (/ time mu)))))
LK-SPROB

該函數可使用:1stmoments, :moments和:margin1方法來近似計算函數的後驗矩和邊際後驗密度。對於這些矩的結果是又有不一樣而且帶有正向的傾斜:

> (send lk :1stmoments #'lk-sprob)
((0.20026961686731173) (0.10757618146651714))
> (send lk :moments #'lk-sprob)
((0.24368499587543843) (0.11571130867154422))

邊際密度能夠這樣表示:

> (plot-lines (send lk :margin1 #'lk-sprob (rseq .01 .8 30)))
#<Object: 13d338c, prototype = SCATTERPLOT-PROTO>

結果以下圖:

圖2.21 WBC=50000的病人一年存活機率的近似邊際後驗密度

基於本圖,數據顯示其存活率幾乎確定低於0.5,可是很難作出一個更精確的說法了。

    在本小節描述的函數式基於前一小節的優化的代碼。默認狀況下,導數是用數值方法計算的。若是你想自行計算導數,你可讓你的log後驗函數返回這樣一個列表,它包含函數值和梯度,或者包含函數值、梯度和Hessian矩陣。在較大的問題裏,在作二階近似的時候可能須要精確的導數值。

練習2.26 

略。

練習2.27

略。

相關文章
相關標籤/搜索