Lisp-Stat翻譯 —— 第五章 統計函數

第五章 統計函數

Lisp-Stat裏的一些統計函數已經在前幾章裏介紹了。這些函數的多數都是處理數值數據集的，表示爲Lisp-Stat組合數據項。本章的第一節將介紹一些額外的函數用來檢測複合數據和矢量運算系統。接下來的一節描述了用來處理數據、相關標準機率分佈的計算和數值線性代數計算的一些函數。這些函數中的不少模仿了S系統裏的類似的函數。本章以兩個列子做爲結束，這兩個例子說明了一些介紹過的工具的用法。

5.1 組合數據

5.1.1 組合數據屬性

組合數據包括列表、矢量、數組合組合數據對象(這個下一章介紹)。非組合數據叫作簡單數據。謂詞compound-data-p能夠用來確認組合數據項。

    組合數據項的主要特徵是它們包含數據項的序列。這個序列可使用函數compound-data-seq來提取。對於列表和矢量該函數僅返回它的參數。對於數組它返回一個元素爲行主序的替換數組。複合數據項的元素自己也能夠是組合的。函數element-seq遞歸地遍歷組合數據項而且返回它的簡單元素序列：

> (compound-data-seq '(a (b c)))
(A (B C))
> (element-seq '(a (b c)))
(A B C)

若是一個組合數據項的全部元素都是簡單數據，那麼element-seq函數與compound-data-seq函數的效果是等價的。

    函數count-elements返回簡單項的數量，即element-seq函數的結果的長度：

> (count-elements '(a (b c)))
3

5.1.2 矢量化計算

在前幾章中，咱們已經簡略點討論了矢量化計算的幾個要點。本小節討論一些與這個主題相關的問題和慣例。

    大多數統計操做都涉及到對數值集合的操做。爲了簡化函數編寫來表現這些操做，Lisp-Sstat重定義了基本Lisp計算函數+、-、*、/、log等等，目的是將組合數據做爲參數，而且返回組合數據項，這些數據項是經過將這些函數做用到組合數據的每個元素上獲得的。

    Lisp提供了一些映射函數來將函數逐元素地做用到列表項或者矢量項上。由於兩個緣由，這些函數不能直接用來實現矢量運算。首先，在表達式(+ '(1 2 3） 4)裏，它的意圖是將4加到列表(1 2 3)的每個元素上。由於參數4不是一個列表，一些像mapcar同樣的函數就不能直接做用在它們上邊。咱們想要一個常量被當作就像是一個常量列表或者合適長度或維度的數組，而不須要必定要強制構建它。

    對於第二個緣由，考慮到如下表達式(+ '(1 2 3) '#(4 5 6))，它的意圖就是逐元素地對兩個序列相加，而後返回結果。可是，一個序列是列表，而另外一個是矢量。那本例中的結果會是什麼呢？map函數須要結果序列的類型必須強制指定。由於這個可能不太方便，在參數類型不一樣時，咱們須要使用一個約定來肯定結果的類型。Lisp-Stat裏使用的約定就是從參數列表的第一個組合數據項裏提取任何須要的結構信息，包括類型、長度、維度等等。其它組合數據參數必須匹配第一個組合數據參數的形式。多個序列，若是它們之間長度相同，不管類型是否相同，都被視爲是有相同形式的。多維數組之間，若是它們的維度是相同的，那它們就是有相同形式的。

    用來處理常量和肯定結果結構的這些規則在函數map-elements的定義中用到過，那個函數是在3.7節介紹的。

    在實現矢量化計算函數的過程裏，還須要作一個決定：像(+ '(1 (2 3)) 4)這樣的表達式如何被求值？有兩個選擇，若是咱們僅容許一個層級的矢量化計算的話，那麼這個表達式將致使一個錯誤。換句話說，咱們能夠遞歸地矢量化加法函數，爲表達式產生(5 (6 7))這樣的結果。哪個選擇是正確的不太清晰，我是決定使用遞歸矢量化的。

5.2 數據處理函數

5.2.1 基礎操做

用來產生成系統的數據的函數，如repeat、iseq和rseq，已經在第二章裏介紹過了。和它們相關的還有函數difference，它做用於數字序列，返回相鄰兩個元素之間的差分：

> (difference '(1 3 6 10))
(2 3 4)
> (difference #(1 3 6 10))
#(2 3 4)

    函數split-list帶一個列表和一個整數長度n，而後將列表分裂成列表的列表，每一個子列表長度爲n：

> (split-list '(1 2 3 4 5 6) 3)
((1 2 3) (4 5 6))

若是列表的長度不是n的倍數，將發送一個錯誤信號。

    函數cumsum帶有一個序列做爲參數，而後返回原序列的累計和序列：

> (cumsum '(1 2 3 4))
(1 3 6 10)
> (cumsum #(1 2 3 4))
#(1 3 6 10)

累積和也可使用accumulate函數來計算，這個函數帶有一個二進制函數和一個序列爲參數，而後返回一個新的結果序列，這個新序列是將二進制函數累積地做用到原序列的每個元素後獲得的結果。那麼累加能夠這樣計算：

> (accumulate #'+ '(1 2 3 4))
(1 3 6 10)

累積則產生：

> (accumulate #'* '(1 2 3 4))
(1 2 6 24)

    一些函數對組合數據項計算簡單總結。除了在第二章介紹的函數mean和standard-deviation，這樣的函數還包括count-elements、sum、prod、max和min。這些函數叫作矢量減小函數，由於它們減小組合數據項成爲單個數字。

    函數min和max能夠帶多個參數，可是它們一般返回全部參數的最大值和最小值。函數pmin和pmax是矢量化的函數，能夠用來並行地計算多個參數的最大值與最小值。例如：

> (pmin '(1 2 3 4 5) '(5 4 3 2 1))
(1 2 3 2 1)

    最後，函數covariance-matrix帶序列或者矩陣爲參數，並將它們造成一個列的列表形式，而後對這個數據列返回樣本的協方差矩陣。全部的數據列都應該是等長度的。

5.2.2 排序函數

函數sort-data帶一個組合數據項參數，返回一個排好序的元素的序列。元素應該都是實數或者字符串。字符串排序是大小寫敏感的。

    order函數將element-seq函數做用到參數上，返回一個最小的和第二小的下標的列表。元素序列的元素。rank函數帶一個組合數據項做爲參數，返回一個相同大小的組合數據項，該數據的條目被它的階數代替：

> (rank '(14 10 12 11))
(3 0 2 1)
> (rank '#2a((14 10) (12 11)))
#2A((3 0) (2 1))

最低的階是0，聯繫被任意打斷。（注：此處不甚瞭解，可能翻譯不精確！！！）

    函數quantile帶兩個參數，分別是一個組合數據項x和一個數字或數字序列p，0≤p≤1，而後返回x的第p級分位數：

> (def x (uniform-rand 20))
X
> (quantile x .5)
0.6490398497045595
> (quantile x '(.2 .8))
(0.44269173085226765 0.8681219493754327)

函數median、interquartile-range和fivnum是使用術語quantile定義的。fivum返回參數的最小值、第一分位數、第二分位數、第三分位數和最大值，即五數歸納法。

5.2.3 插值和平滑

函數spline帶連個等長度的序列x和y做爲參數，而後針對這兩個序列裏指定的點返回一個天然三次樣條插值的x與y值的列表。參數x裏的值必須嚴格遞增。默認狀況下，插值函數將在x值的範圍內計算30個等距離點。這個能夠經過使用關鍵字:xvals來改變。若是使用了這個關鍵字，其值爲整數n的話，那麼結果會使用n個等距離點來代替30等距離點。所以表達式

(let ((x (rseq 0 pi 10)))
    (spline x (sin x) :xvals 50))

將使用50個點返回正弦函數的插值。若是使用了:xvals關鍵字，其值爲一個序列的話，那麼這個序列將被用做計算樣條時的那個x值。因此咱們的50點插值也能夠這麼構建：

(let ((x (rseq 0 pi 10)))
    (spline x (sin x) :xvals (rseq 0 pi 50)))

    lowess函數使用LOWESS平滑算法做用到兩個序列值x和y上。它返回x序列的列表和y序列的平滑值。因此表達式

(setf p (plot-points x y))
(setf p :add-lines (lowess x y))

將使用變量x和y的數據構建一個散點圖，而後使用LOWESS算法對圖形進行平滑處理。

    LOWESS算法的一些參數能夠經過向lowess函數船體關鍵字參數來控制。用來估計每一個點的擬合程度的數據的分數經過關鍵字:f來指定；默認值爲0.25。穩定迭代的數量可使用關鍵字:steps來指定，默認值爲2。最後，在給定點處的擬合使用線性插值肯定的。這個臨界值可使用關鍵字:delta來改變，默認值爲x值的範圍的2%。

    函數kernel-smooth提供了一個備用的平滑函數。它以x和y序列爲參數，返回x值的列表和平滑後的y值的列表的集合。默認地，結果的x值是在輸入的x值的範圍上的30份的等距離點。針對spline函數，關鍵字:xvals能夠用來指定點數量的備選值，或者是針對結果的x值的序列的備選值。核心的4個類型都是被支持的，指定符號G表明高斯，指定符號T表明三角形，指定符號U表明正態，指定符號B表明二項分佈。默認核是二項分佈，可是你可使用關鍵字:type來指定備用的核。核的寬度可使用關鍵字:keyword來指定。所以表達式

(send p :add-lines (kernel-smooth x y :type 'g))

將向上邊的繪圖結構里加入一個高斯核平滑。

    核密度估計量也是可用的。kernel-dens須要一個單獨的x值的序列做爲參數，返回x值的列表和密度估計值。表達式

> (plot-lines (kernel-dens (normal-rand 50)))
#<Object: 1431e04, prototype = SCATTERPLOT-PROTO>

針對50標準正態隨機變量組成的樣本，產生一個核密度預測的圖形。kernel-dens函數接受與kernel-smooth函數相同的關鍵字參數。

5.3 機率分佈

Lisp-Stat含有一些函數，用來評估一些連續單變量分佈的密度、累積分佈函數和分位數，包括正態分佈、柯西分佈、β-分佈、γ-分佈、x²-分佈、t-分佈和F-分佈。

    用來計算二項式分佈和泊松分佈的機率密度函數、累積密度函數和分位數函數的相似函數是可用的。針對這些離散分佈的quantile函數被定義成以下形式：累積分佈函數的左連續逆。機率密度函數僅被定義爲整型參數。最後，還有一個用來估計二元正態分佈的累積分佈函數的估值的函數，這些可用的函數見下表，全部的函數都是矢量化的：

針對正態分佈和柯西分佈的函數假設是標準分佈，所以只須要一個參數，這裏的值將來計算密度函數或者累積分佈函數，這裏的0和1之間的比例係數，用來估計分位數函數。針對γ-分佈、x²-分佈、t-分佈和泊松分佈的函數還須要一個額外的參數：分別是γ指數、x²和t分佈的自由度，還有泊松均值。假定γ-分佈有一個尺度爲1的參數。這有幾個例子：

> (chisq-cdf 7.5 8)
0.5162326184463124
> (t-quant .975 3)
3.1824463052837086
> (poisson-pmf 3 2.5)
0.2137630172497364

    針對β-分佈、F-分佈和二項式分佈的函數須要兩個額外的參數：分別是β-分佈的兩個指數、F-分佈的分子和分母的自由度，二項式分佈的樣本大小與成功的可能性。針對這個也有幾個例子：

> (beta-cdf .3 5.2 6.3)
0.14839195210538714
> (f-quant .9 3 7)
3.0740719939090018
> (binomial-pmf 3 5 .5)
0.3124999999999998

    函數bivnorm-cdf假設兩個邊際都是標準正態的，它帶三個參數，累積分佈函數須要的x和y參數，還有相關性係數：

> (bivnorm-cdf .3 .2 .6)
0.4551159427475644

    在第2.4.1節裏描述的這些分佈的隨機數字生成器是使用Common Lisp隨機數生成器實現的。Common Lisp隨機數接口是打算要跨不一樣系統的。咱們使用的特定的生成器可能隨着咱們使用系統的不一樣產生不一樣的數值。生成器的種子一般使用系統時鐘生成，能夠保存和恢復，也能夠生成新的種子。可是卻不能爲生成器指定一個任意種子，由於一個適當的選擇是依據生成器的實現、機器的字長等等因素的。

    狀態的當前值保存在全局變量*random-state*裏。函數make-random-state可用來設置和保存狀態。該函數帶有一個可選參數，吐過這個參數是nil或者被忽略的話，那麼make-random-state函數將返回變量*random-state*的當前值的一份拷貝。若是參數是一個隨機狀態對象，將返回它的拷貝。最後，若是參數是符號t的話，將產生一個新的隨機化的初始化狀態對象而後返回它。一般使用系統時鐘產生新的種子，可是精確的機制是隨之Lisp系統的不一樣而不一樣的。

    一個帶打印表達式的隨機狀態對象可使用read函數讀回。所以你能夠將當前狀態保存到一個文件裏，而後使用它重複一個帶有僞隨機數的相同流的模擬過程。例如，使用表達式將當前狀態寫入到文件以後：

(with-open-file (s "currentstate" :direction :output)
                  (print *random-state* s))

如下表達式

(with-open-file (s "currentstate")
                  (setf *random-state* (read s))
                  (run-simulation))

每求值一次，都會使用相同的隨機數字集合來重複模擬器。

5.4 數組和線性代數函數

5.4.1 基本矩陣和數組函數

矩陣是二維數組，在統計計算中扮演一個重要的角色，因此Lisp-Stat提供了一些附加函數來幫助管理矩陣。謂詞matrixp遇到一個矩陣時返回t，遇到其它類型數據將返回nil。

    有時候，可以將列表或者矢量集合拼接到成矩陣的行或列是頗有用的。可使用函數bind-rows和bind-columns：

> (bind-rows '#(1 2 3) '(4 5 6))
#2A((1 2 3) (4 5 6))
> (bind-columns '#(1 2 3) '(4 5 6))
#2A((1 4) (2 5) (3 6))

    傳遞給bind-rows和bind-columns函數的參數多是矩陣：

> (bind-rows '#2a((1 2 3) (4 5 6)) '(7 8 9))
#2A((1 2 3) (4 5 6) (7 8 9))

這個能夠被利用來編寫一個函數，該函數使用兩個矢量和一個常量造成矩陣的新邊界：

> (defun border-matrix (a b c d)
    (bind-rows (bind-columns a b)
               (concatenate 'list c (list d))))
BORDER-MATRIX

而後，舉個例子：

> (border-matrix '#2a((1 2) (3 4)) '(5 6) '(7 8) 9)
#2A((1 2 5) (3 4 6) (7 8 9))

    例如，爲了可以計算列的平均值，你須要將矩陣拆散成列。column-list函數就是執行這個慚怍的：

> (column-list '#2a((1 2 3) (4 5 6)))
(#(1 4) #(2 5) #(3 6))

如今，可使用mapcar函數定義column-mean函數了：

> (defun column-means (x)
    (mapcar #'mean (column-list x)))
COLUMN-MEANS

函數row-list也以一樣的方式運。

    還有四個函數也值得一提。transpose函數帶一個矩陣做爲參數，返回參數的轉置。transpose函數也能夠用來表示具備相同長度的列表集的列表的轉置，表示矩陣的行。而後它返回列表集的列表，表明轉置的行。函數identity-matrix帶一個整型參數k，返回一個k乘k的特徵矩陣。函數diagonal有兩種用法。傳遞給它一個矩陣做爲參數，它返回矩陣對角線項的列表(這個矩陣必須是方陣)：

> (diagonal '#2a((1 2) (3 4)))
(1 4)

若是傳遞給它一個序列參數，diagonal返回一個方陣，該方陣的對角線上的元素就是序列的元素，其它元素爲0。

> (diagonal '(a b c))
#2A((A 0 0) (0 B 0) (0 0 C))

    最後，print-matrix函數能夠用來打印一個比傳統打印機方法更好的矩陣表達式：

> (print-matrix #2a((1 2 3) (4 5 6) (7 8 9)))
#2a(
    (1     2     3    )
    (4     5     6    )
    (7     8     9    )
   )
NIL

與print函數相似，函數print-matrix也接受可選參數。

    Lisp-Stat還提供了一個函數permute-array，它帶一個數組和一個列表做爲參數，而後根據列表指定的轉置對數組進行轉置：

> (permute-array '#2a((1 2 3) (4 5 6)) '(1 0))
#2A((1 4) (2 5) (3 6))
> (permute-array '#3a(((1 2 3) (4 5 6)) ((7 8 9) (10 11 12)))
                 '(2 0 1))
#3A(((1 4) (7 10)) ((2 5) (8 11)) ((3 6) (9 12)))

5.4.2 矩陣乘法

函數帶兩個矩陣爲參數，並對它們運行乘法操做。矩陣的元素必須是數字。

> (matmult '#2a((1 2 3) (4 5 6)) '#2a((1 2) (3 4) (5 6)))
#2A((22.0 28.0) (49.0 64.0))

參數中的一個也能夠是一個序列——一個列表或者一個矢量。若果是這種狀況，而且它是第一個參數的話，這個序列被當成一個行矢量；若是它是第二個參數的話，它將被當成列矢量。結果是與序列參數是相同類型的一個序列：

> (matmult '#(1 2) '#2a((1 2 3) (4 5 6)))
#(9.0 12.0 15.0)
> (matmult '#2a((1 2 3) (4 5 6)) '(1 2 3))
(14.0 32.0)

    若是傳遞給函數matmult的參數都是具備相同長度的序列，將返回一個數值，也就是兩個參數的內積。內積也可使用inner-product函數來計算求得：

> (matmult '(1 2 3) '#(4 5 6))
32.0
> (inner-product '(1 2 3) '#(4 5 6))
32.0

    matmult函數能夠經過傳遞多於兩個參數來調用。在這種狀況下，前兩個參數相乘，其結果再與右邊的第三個參數相乘，等等。

    函數outer-product函數帶兩個序列做爲參數，返回他們的外積：

> (outer-product '#(1 2) '(3 4 5))
#2A((3 4 5) (6 8 10))

outer-product提供了一個通常意義上的外積。它有一個可選參數，該可選參數是一個帶兩個參數的函數。若是outer-product函數使用了這個可選參數，那麼它的結果的(i, j)位置的元素是這樣的兩個數的乘積組成的：第一個序列的第i個元素和第二個序列的第j個元素：

> (outer-product '#(1 2) '(3 4 5) #'+)
#2A((4 5 6) (5 6 7))

    函數cross-product帶一個參數x，返回(matmult (transpose x) x)的值。

5.4.3 線性代數函數

Lisp-Stat提供了一些進行線性代數運算的函數。其中一個常常會被其它函數調用的函數就是lu-decomp函數。該函數帶一個數值方陣，而後運行帶局部定支點的LU分解，最後返回4個元素的列表。第一個元素是一個包含LU分解的方陣，該方陣定支點過程肯定的那個矩陣的以行爲主的排列的LU分解方陣。除對角線外，下三角矩陣是L部分；L的對角線元素全爲1。結果矩陣的上三角部分，包括對角線，是U部分。結果列表的第二個元素是一個整數的矢量，用來描述定支點過程裏用到的排列。第三個元素是一個數字，它等於±1，若是用來進行航互換的數目爲偶數，它就是正1,；是技術的話，該參數就取負1。結果矩陣的最後一個元素是一個標誌，若是輸入矩陣是數值上的奇異矩陣，該值爲t；若是是非奇異矩陣，該值爲nil。

> (lu-decomp '#2a((1 2) (3 4)))
(#2A((3.0 4.0) (-0.3333333333333333 0.6666666666666667)) #(1 1) -1.0 NIL)

    lu-solve函數帶兩個參數。第一個參數是由函數lu-decomp使用矩陣A返回的結果。第二個參數是一個列表或者矢量，它是表達式Ax = b的右值b。lu-solve函數使用帶LU分解功能的back-solving解法解該方程，而後返回結果。

> (lu-solve (lu-decomp '#2a((1 2) (3 4))) '(5 6))
(-3.9999999999999987 4.499999999999999)
> (matmult '#2a((1 2) (3 4)) '(-4 4.5))
(5.0 6.0)

    在back-solving算法處理期間，若是矩陣被發現是數值上奇異的，將發生一個錯誤信號。

    determinant、inverse和solve這幾個函數將使用LU分解來處理它們的矩陣參數。determinant函數帶有一個矩陣做爲參數，而後返回它的行列式。inverse返回參數矩陣的逆矩陣，若是這個矩陣是數值奇異的則發送一個錯誤信號。

> (inverse '#2a((1 2) (3 4)))
#2A((-1.9999999999999996 0.9999999999999998) (1.4999999999999998 -0.4999999999999999))
> (matmult '#2a((1 2) (3 4)) '#2a((-2 1) (1.5 -0.5)))
#2A((1.0 0.0) (0.0 1.0))

solve函數帶有一個矩陣A和一個參數b做爲參數，而後返回Ax = b的解。b能夠是一個序列或者一個矩陣。

> (solve '#2a((1 2) (3 4)) '(5 6))
(-3.9999999999999987 4.499999999999999)
> (solve '#2a((1 2) (3 4)) (identity-matrix 2))
#2A((-1.9999999999999998 0.9999999999999999) (1.4999999999999998 -0.4999999999999999))

lu-decomp、lu-solve、determinant、inverse和solve等函數的參數能夠包含實數或者複數。

    chol-decomp函數帶有一個對稱方陣A做爲參數，而後計算矩陣A+D的cholesky分解，這裏的D是一個非負的對角矩陣，它被加到A矩陣上，若是能確保A+D是數值上嚴格正值定義。由函數chol-decomp返回的結果是下三角矩陣L的一個列表，矩陣L知足LL^T= A+D，該返回值是D的元素裏的最大值。若是矩陣A是嚴格的正值定義，結果的第二個元素爲0：

> (setf x (chol-decomp '#2a((2 1) (1 2))))
(#2A((1.4142135623730951 0.0) (0.7071067811865475 1.224744871391589)) 0.0)
> (matmult (first x) (transpose (first x)))
#2A((2.0000000000000004 1.0) (1.0 1.9999999999999996))

    qr-decomp帶一個m行n列的矩陣A做爲參數，這裏m≥n，而後計算其QR分解，返回的結果是一個列表，該列表是由兩部分組成的：一個帶正交列的m行n列矩陣Q，一個上三角矩陣R，好比說QR = A：

> (setf x (qr-decomp '#2a((1 2) (3 4) (5 6))))
(#2A((-0.16903085094570325 0.8970852271450608) (-0.50709255283711 0.2760262237369412) (-0.8451542547285165 -0.34503277967117696)) #2A((-5.916079783099616 -7.4373574416109465) (0.0 0.8280786712108248)))
> (matmult (first x) (second x))
#2A((0.9999999999999997 1.9999999999999996) (3.0 4.0) (4.999999999999999 6.0))

咱們可使用print-matrix函數來使結果更有可讀性：

> (print-matrix (first x))
#2a(
    (-0.169031      0.897085    )
    (-0.507093      0.276026    )
    (-0.845154     -0.345033    )
   )
NIL
> (print-matrix (second x))
#2a(
    ( -5.91608      -7.43736    )
    ( 0.000000      0.828079    )
   )
NIL

    若是函數qr-decomp使用了值爲t的:pivot關鍵字，它將從新排列矩陣的列以確保矩陣R的對角線元素的絕對值是沒有增加的。在這種狀況下，返回的結果裏將包含第三個元素，也就是按它們被使用的順序的那些列的下標號。

    sv-decomp函數帶有一個m行n列的矩陣A做爲參數，這裏的m≥n，而後計算它的奇異值分解。該函數返回(U W V FLAG)形式的列表，這裏的U和V都是矩陣，它們的列是左奇異和右奇異矢量，W是遞減序列上的矩陣A的奇異值矢量。FLAG在算法收斂的時候，值爲t，不然值爲nil。若是FLAG是t，而且矩陣D是一個W爲對角線元素，其它位置爲0的一個對角線矩陣的話，那麼UDV^T= A：

> (setf x (sv-decomp '#2a((1 2) (3 4) (5 6))))
(#2A((-0.22984769640007144 0.8834610176985256) (-0.5247448187602937 0.2407824921325462) (-0.8196419411205158 -0.4018960334334318)) #(9.525518091565107 0.5143005806586441) #2A((-0.6196294838293404 -0.7848944532670524) (-0.7848944532670524 0.6196294838293404)) T)
> (matmult (first x)
           (diagonal (second x))
           (transpose (third x)))
#2A((0.9999999999999994 1.9999999999999993) (3.0 3.9999999999999996) (4.999999999999999 6.0))

    函數eigenvalues 以遞減的順序 返回一個特徵值的向量；函數eigenvectors返回一個對應於一個對稱方陣的特徵向量的列表。eigen函數返回一個列表，它的兩個元素是：一個特徵值的向量，一個對稱方陣的特徵向量的列表。函數backsolve帶有一個上三角矩陣A和一個序列b爲參數，返回Ax = b的解。函數rcondest返回一個上三角矩陣的條件數的倒數的估計值。（注：矩陣的條件數就是 一個矩陣的最大奇異值與其最小奇異值之比。）

    最後一個函數在構造某些動態圖形時是很是有用的，它就是make-rotation函數。這個函數帶有兩個等長度的序列和一個可選的實數做爲參數，該實數表示一個角的弧度數。它返回旋轉矩陣，該旋轉矩陣在由兩個序列定義的平面裏旋轉而成，並保持從第一個序列到第二個序列的正交補不變。例如：

> (make-rotation '(1 0 0 0) '(1 1 1 1) (/ pi 10))
#2A((0.9510565162951535 -0.17841104488654497 -0.17841104488654497 -0.17841104488654497) (0.17841104488654497 0.9836855054317178 -0.01631449456828216 -0.01631449456828216) (0.17841104488654497 -0.01631449456828216 0.9836855054317178 -0.01631449456828216) (0.17841104488654497 -0.01631449456828216 -0.01631449456828216 0.9836855054317178))

以上表達式返回一個矩陣，該矩陣在第一個座標軸和一個四維空間的對角線的平面裏以π/10的角度旋轉。若是make-rotation函數沒有給定角度參數，它將使用兩個序列中較小的角度。

    Lisp-Stat內部的線性代數函數使用double-float或者long-float類型的浮點數，使用哪一個依賴Lisp系統和硬件。全局變量machine-epsilon合適類型的最小的浮點數x，所以，表達式(not (= 1 (+ 1 x)))將返回false，即：

> (setf x machine-epsilon)
1.0842021724855044E-19
> (not (= 1 (+ 1 x)))
NIL

5.5 迴歸計算函數

迴歸計算可使用文獻65裏的問題2.7提到的消除操做來執行。消除是按行執行的。也就是說，在中心點k上清除一個矩陣A就是使用給定的元素將它轉換爲矩陣B：

清除操做的定義是累積的、自反相的：兩個中心點被清除的地方，它們的位置是可有可無的，在同一點上進行兩次清除操做與在該中心點沒有進行清除操做是等價的。（注：該段理解不清，翻譯可能不許確，後期研究清除再補上！）

    函數mak-sweep-matrix帶一個矩陣x、一個序列y和一個可選的權重的序列w，返回校訂的叉積矩陣，該叉積矩陣已經在常量行上被消除掉了。也就是說，若是A是由一個只有1的列，一個x的多列和列y組成的。矩陣W是一個對角線矩陣，該矩陣列對角線上的元素就是w。那麼，make-sweep-matrix函數計算叉積矩陣A^TWA，而後在常數行上消除它。

    函數sweep-operator帶一個矩陣，一個要消除的行座標的序列和一個可選的容插值序列做爲參數，它返回一個列表，該列表由在指定行消除的矩陣的一份複製和事實上要消除的行的下標的列表組成的。僅當中心點的絕對值比指定的偏差容差值大的時候才消除，默認容差值是0.000001。

    若是一個叉積矩陣已經在一箇中心點集合上被消除，那麼被消除的矩陣的最下邊一行的對應元素將是由消除的中心點索引的變量組成的那個模型的最小二乘係數。那麼一個帶x矩陣但排除其常數項，一個y向量和一個權值矩陣序列爲參數，而且用來計算最小二乘係數的函數能夠這樣定義：

> (defun reg-coefs (x y weights)
    (flet ((as-list (x) (coerce (compound-data-seq x) 'list)))
      (let* ((m (make-sweep-matrix x y weights))
             (p (array-dimension x 1)))
        (as-list (select (first (sweep-operator m (iseq 1 p)))
                         (+ p 1)
                         (iseq 0 p))))))
REG-COEFS

select表達式返回的結果是一個1行p+1列的矩陣，首先使用函數compound-data-seq將該矩陣轉換爲一個矢量，而後使用函數coerce轉化爲一個列表。該定義假設模型包含一個截距。對於沒有截距的模型，咱們須要首先消除下標0來移除常數項，僅僅選擇下標爲1, 2, 3, ...,p的項做爲結果。

5.6 一些例子

本節將給出兩個額外的例子，它們使用是使用本章裏介紹的一些工具完成的。第一個例子是咱們早前使用的make-projection函數的稍微顯得更加詳盡的版本；第二個例子是一個函數，該函數的用途是針對大量的權值函數計算其穩定的迴歸係數。

5.6.1 構造一個投影操做符

在第三章裏咱們看到了一些關於函數的變量，它們使用一個函數閉包創建了一個投影操做符的模型。那時咱們僅能較容易地處理一維空間上的投影，可是如今咱們擁有一些能夠處理高維投影的工具，一個方法使用列表的列表或者向量的方式生成一個空間，使用它們做爲矩陣X的列，而後造成X的特徵值分解矩陣UDV^T。矩陣U的列對應正特徵值，而後造成一個關於X的列空間的正交基。所以若是U₁ 是這些列的矩陣，那麼向量y在X列空間上的投影P_y能夠這樣表示：P_{y =}  U₁U₁^Ty 。

    在make-projection函數的新版本的定義裏，若是須要，咱們經過使用表達式(if (matrixp cols) cols (apply #'bind-columns cols))，容許指定的列以向量的列表或者一個矩陣的形式轉換成一個矩陣。

    在使用u計算矩陣U，而且使用s-vals來計算特徵值以後，咱們能夠這樣計算非零特徵值的列集列表：(select (column-list u) (which (> s-vals 0.0)))。這些列能夠經過使用bind-columns函數造成矩陣U₁。

    在最終的投影計算裏，咱們能夠避免在U₁的左側乘一個列表或者向量y來對U₁進行轉置，咱們可使用結果，該結果也多是一個列表或者向量，將它放到U₁的右側，就像下邊的表達式：(matmult u1 (matmult y u1))。

    合併這些步驟，咱們獲得如下定義：

> (defun make-projection (cols)
    (let* ((x (if (matrixp cols)
                  cols
                  (apply #'bind-columns cols)))
           (svd (sv-decomp x))
           (u (first svd))
           (s-vals (second svd))
           (basis (select (column-list u)
                          (which (> s-vals 0.0))))
           (u1 (apply #'bind-columns basis)))
      #'(lambda (y) (matmult u1 (matmult y u1)))))
MAKE-PROJECTION

練習 5.1

略。

5.6.2 穩健迴歸

對一個形如y = Xβ的線性模型，計算其穩健迴歸的方法是使用一個能夠交互地重定義權值的最小二乘算法，該算法裏的權值是使用一個文件的權值函數w計算的，該權值矩陣對有較大標準化殘差的觀測值賦予較小的權重。特別地，給定一個權重的集合，咱們能夠：

• 計算帶權值的最小二乘估計量β帽(上邊帶一個尖帽，表示β的估計量)
• 計算殘差r
• 計算一個範圍的估計量s帽(上邊帶一個尖帽，表示s的估計量)，做爲被0.6745劃分的中位數絕對殘差
• 計算權值的一個新的集合，w(r/(s帽))

而後這些步驟將重複運行，直到知足收斂準則。範圍估計量包括除以0.6745，以確保它是對正態分佈數據的作出的全局標準差估計。

    在文獻中，有不少不一樣的權值函數被提出，包括biwieight函數：

柯西權值函數：

還有Huber權值函數：

這三個參數都依賴一個調整參數K。對於biweight函數，K的建議值是4.685；對於柯西權值函數，K的建議值是2.385；對於Huber權值函數，K的建議值是1.345。

    在實現一個穩健迴歸函數過程當中，假設模型裏包含一個截距，還有一個在指定的X矩陣裏的列集，咱們從這裏開始。對於給定的權值集合，帶權值的最小二乘估計能夠像使用函數reg-coefs函數同樣計算，該函數是在5.5節裏定義的。若是咱們得到了參數集合，記作bate列表，咱們就能夠這樣計算它的擬合值：
(+ (first beta) (matmult x (rest beta)))。

    可使用一個do循環來控制迭代，就像這樣：

(do ((lase nil beta)
       (count 0 (+ count 1))
       (beta (reg-coefs weights) (improve-guess beta)))
      ((good-enough-p last beta count) beta))

迭代變量last用來處理係數集的前一個值，目的是與當前的猜想值做比較；nil的初始值用來講明沒有可用的前一個值了。迭代的總次數也被監測。

    若是變量wf引用了權值函數，那麼improve-guess函數的函數體能夠這樣編寫：

(let* ((resids (- y (fitvals beta)))
         (scale (/ (median (abs resids)) 0.6745))
         (wts (funcall wf (/ resids scale))))
    (reg-coefs wts))

最後，若是迭代次數超出給定範圍，good-enough-p函數應該打印一條消息，而且檢查參數估計量的改變是否足夠的小。能夠這樣編寫測試表達式：

(flet ((rel-err (x y) (mean (/ (abs (- x y)) (+ 1 (abs x))))))
    (or (> count count-limit)
        (and last (< (rel-err beta last) tol))))

這裏的and表達式首先檢查last變量是否爲nil，而後用參數列表裏的平均相對變化與容差值tol的比較。在改變計算時加1是要保證對大的係數以相對偏差量度，對小的係數以絕對偏差量度。其它的選擇是明顯可能的。

    如今咱們能夠將這些元素組合到一個單一函數裏，使用一個labels表達式來定義一些函數用以實現每一步。使用關鍵字參數指定初始權重、收斂法則的容差值和迭代次數的邊界值：

> (defun robust-coefs (xy wf &key
                          (weights (repeat 1 (length y)))
                          (tol 0.0001)
                          (count-limit 20))
    (let ((x (if (matrixp x) x (apply #'bind-columns x))))
      (labels ((as-list (x) (coerce (compound-data-seq x) 'list))
               (rel-err (x y)
                        (mean (/ (abs (- x y)) (+ 1 (abs x)))))
               (reg-coefs (weights)
                          (let* ((m (make-sweep-matrix x y weights))
                                 (p (array-dimension x 1)))
                            (as-list
                             (select (first (sweep-operator m (iseq 1 p)))
                                     (1+ p)
                                     (iseq 0 p)))))
               (firstvals (beta)
                          (+ (first beta) (matmult x (rest beta))))
               (improve-guess (beeta)
                              (let* ((resids (- y (fitvals beta)))
                                     (scale (/ (median (abs resids)) 0.6745))
                                     (wts (funcall wf (/ resids scale))))
                                (reg-coefs wts)))
               (good-enough-p (last beta count)
                              (if (> count count-limit)
                                  (format t "Iteration limit exceeded~%"))
                              (or (> count count-limit)
                                  (and last (< (rel-err beta last) tol)))))
        (do ((last nil beta)
             (count 0 (+ count 1))
             (beta (reg-coefs weights) (improve-guess beta)))
            ((good-enough-p last beta count) beta)))))
ROBUST-COEFS

    爲了使用該函數，你須要指定一個權值函數。爲了容易地使用上邊列出的標準權值函數中的一個，咱們可使用以下定義的函數：

> (defun make-wf (name &optional
                       (k (case name
                            (biweight 4.685)
                            (cauchy 2.385)
                            (huber 1.345))))
    #'(lambda (x)
        (leet ((u (abs (/ r k))))
              (case name
                (biweight (^ (- 1 (^ (pmin u 1) 2)) 2))
                (cauchy (/ 1 (+ 1 (^ u 2))))
                (huber (/ 1 (pmax u 1)))))))
MAKE-WF

    做爲說明，咱們能夠用一個例子測試一下這個函數，來自於布朗利的stack loss數據(注： stack loss指從煙道里排除的浪費的熱量、不可再生能源和過量的氣體，它們來自住宅區通風系統、工廠或者政府建築物設施)。該數據集由21組觀測值組成，它們是關於工廠進行氨氧化成硝酸的操做過程。三個獨立變量時AIR，即經過工廠的氣流，表明該過程的執行速度；TEMP，即經過吸取塔內部盤管裏的循環的製冷水的溫度；CONC，一個線性酸循環的濃度的編碼。獨立變量LOSS，即10次飄到到廠房的氨氣的百分比，也就是從吸取塔裏逸出的未吸取的部分。

    最小二乘擬合能夠這樣計算：

表5.1 煙道損失數據

> (def air '(80 80 75 62 62 62 62 62 58 58 58 58 58 58 50 50 50 50 50 56 70))
AIR
> (def temp '(27 27 25 24 22 23 24 24 23 18 18 17 18 19 18 18 19 19 29 29 29))
TEMP
> (def conc '(89 88 90 87 87 87 93 93 87 80 89 88 82 93 89 86 72 79 80 82 91))
CONC
> (def loss '(42 37 37 28 18 18 19 20 15 14 14 13 11 12 8 7 8 8 9 15  15))
LOSS
> (regression-model (list air temp conc) loss)

Least Squares Estimates:

Constant                  -34.8780      (15.8781)
Variable 0                 1.04201      (0.139894)
Variable 1                 8.517019E-2  (0.270213)
Variable 2               -0.144541      (0.206340)

R Squared:                0.851471    
Sigma hat:                 4.25193    
Number of cases:                21
Degrees of freedom:             17

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

使用Huber權值函數的穩健迴歸模型的係數是：

> (robust-coefs (list air temp conc) loss (make-wf 'huber))
(-41.0267 0.829344 0.926232 -0.127856)

    已經獲取這些係數了，咱們能夠返回去計算在最後的擬合裏要用到的權值，而且檢查是否存在某些點收到的權重特別的小。在這種狀況下，收到的權重小於0.5的點的下標分別是0、二、3和20，這些點一般被當成該數據集裏可能的異常值。

練習 5.2

略。

練習 5.3

略。