廣義極小殘差法演算法

1. 機器學習中的損失函數

機器學習中的損失函數
損失函數（loss function）是用來估量你模型的預測值f(x)與真實值Y的不一致程度，它是一個非負實值函數,通常使用L(Y, f(x))來表示，損失函數越小，模型的魯棒性就越好。損失函數是經驗風險函數的核心部分，也是結構風險函數重要組成部分。模型的結構風險函數包括了經驗風險項和正則項，通常可以表示成如下式子：

其中，前面的均值函數表示的是經驗風險函數，L代表的是損失函數，後面的是正則化項（regularizer）或者叫懲罰項（penalty term），它可以是L1，也可以是L2，或者其他的正則函數。整個式子表示的意思是找到使目標函數最小時的值。下面主要列出幾種常見的損失函數。
一、log對數損失函數（邏輯回歸）
有些人可能覺得邏輯回歸的損失函數就是平方損失，其實並不是。平方損失函數可以通過線性回歸在假設樣本是高斯分布的條件下推導得到，而邏輯回歸得到的並不是平方損失。在邏輯回歸的推導中，它假設樣本服從伯努利分布（0-1分布），然後求得滿足該分布的似然函數，接著取對數求極值等等。而邏輯回歸並沒有求似然函數的極值，而是把極大化當做是一種思想，進而推導出它的經驗風險函數為：最小化負的似然函數（即max F(y, f(x)) —-> min -F(y, f(x)))。從損失函數的視角來看，它就成了log損失函數了。
log損失函數的標准形式：
L(Y,P(Y|X))=?logP(Y|X)L(Y,P(Y|X))=?log?P(Y|X)剛剛說到，取對數是為了方便計算極大似然估計，因為在MLE中，直接求導比較困難，所以通常都是先取對數再求導找極值點。損失函數L(Y, P(Y|X))表達的是樣本X在分類Y的情況下，使概率P(Y|X)達到最大值（換言之，就是利用已知的樣本分布，找到最有可能（即最大概率）導致這種分布的參數值；或者說什麼樣的參數才能使我們觀測到目前這組數據的概率最大）。因為log函數是單調遞增的，所以logP(Y|X)也會達到最大值，因此在前面加上負號之後，最大化P(Y|X)就等價於最小化虛鏈L了。
邏輯回歸的P(Y=y|x)表達式如下：P(Y=y|x)=11+exp(?yf(x))P(Y=y|x)=11+exp(?yf(x))

將它帶入到上式，通過推導可以得到logistic的損失函數表達式，如下：

L(y,P(Y=y|x))=log(1+exp(?yf(x)))L(y,P(Y=y|x))=log?(1+exp(?yf(x)))
邏輯回歸最後得到的目標式子如下：

如果是二分類的話模譽仿，則m值等於2，如果是多分類，m就是相應的類別總個數。這里需要解釋一下：之所以有人認為邏輯回歸是平方損失，是因為在使用梯度旦纖下降來求最優解的時候，它的迭代式子與平方損失求導後的式子非常相似，從而給人一種直觀上的錯覺。
這里有個pdf可以參考一下：Lecture 6: logistic regression.pdf.
二、平方損失函數（最小二乘法, Ordinary Least Squares ）
最小二乘法是線性回歸的一種，OLS將問題轉化成了一個凸優化問題。在線性回歸中，它假設樣本和雜訊都服從高斯分布（為什麼假設成高斯分布呢？其實這里隱藏了一個小知識點，就是中心極限定理，可以參考【central limit theorem】），最後通過極大似然估計（MLE）可以推導出最小二乘式子。最小二乘的基本原則是：最優擬合直線應該是使各點到回歸直線的距離和最小的直線，即平方和最小。換言之，OLS是基於距離的，而這個距離就是我們用的最多的歐幾里得距離。為什麼它會選擇使用歐式距離作為誤差度量呢（即Mean squared error， MSE），主要有以下幾個原因：
簡單，計算方便；
歐氏距離是一種很好的相似性度量標准；
在不同的表示域變換後特徵性質不變。
平方損失（Square loss）的標准形式如下：
(Y,f(X))=(Y?f(X))2L(Y,f(X))=(Y?f(X))2
當樣本個數為n時，此時的損失函數變為：

Y-f(X)表示的是殘差，整個式子表示的是殘差的平方和，而我們的目的就是最小化這個目標函數值（註：該式子未加入正則項），也就是最小化殘差的平方和（resial sum of squares，RSS）。
而在實際應用中，通常會使用均方差（MSE）作為一項衡量指標，公式如下：
MSE=1n∑i=1n(Yi~?Yi)2MSE=1n∑i=1n(Yi~?Yi)2
上面提到了線性回歸，這里額外補充一句，我們通常說的線性有兩種情況，一種是因變數y是自變數x的線性函數，一種是因變數y是參數的線性函數。在機器學習中，通常指的都是後一種情況。
三、指數損失函數（Adaboost）
學過Adaboost演算法的人都知道，它是前向分步加法演算法的特例，是一個加和模型，損失函數就是指數函數。在Adaboost中，經過m此迭代之後，可以得到:

Adaboost每次迭代時的目的是為了找到最小化下列式子時的參數 和G：

而指數損失函數(exp-loss）的標准形式如下

可以看出，Adaboost的目標式子就是指數損失，在給定n個樣本的情況下，Adaboost的損失函數為：

關於Adaboost的推導，可以參考Wikipedia：AdaBoost或者《統計學習方法》P145.
四、Hinge損失函數（SVM）
在機器學習演算法中，hinge損失函數和SVM是息息相關的。在線性支持向量機中，最優化問題可以等價於下列式子：

下面來對式子做個變形，令：

於是，原式就變成了：

如若取，式子就可以表示成：

可以看出，該式子與下式非常相似：

前半部分中的就是hinge損失函數，而後面相當於L2正則項。
Hinge 損失函數的標准形式
可以看出，當|y|>=1時，L(y)=0。
更多內容，參考Hinge-loss。
補充一下：在libsvm中一共有4中核函數可以選擇，對應的是-t參數分別是：
0-線性核；
1-多項式核；
2-RBF核；
3-sigmoid核。
五、其它損失函數
除了以上這幾種損失函數，常用的還有：
0-1損失函數

絕對值損失函數

下面來看看幾種損失函數的可視化圖像，對著圖看看橫坐標，看看縱坐標，再看看每條線都表示什麼損失函數，多看幾次好好消化消化。

OK，暫時先寫到這里，休息下。最後，需要記住的是：參數越多，模型越復雜，而越復雜的模型越容易過擬合。過擬合就是說模型在訓練數據上的效果遠遠好於在測試集上的性能。此時可以考慮正則化，通過設置正則項前面的hyper parameter，來權衡損失函數和正則項，減小參數規模，達到模型簡化的目的，從而使模型具有更好的泛化能力。

2. 二維網格數值插值技術

大部分油氣藏的數據是散亂分布的，因此稱為散亂數據。散亂數據指的是在二維平面上或三維空間中，無規則的、隨機分布的數據。利用散亂數據建模就要對散亂數據進行插值或擬合。

設在二維平面上有n個點 （x_i，y_i） （i ＝1，…，n），並有Z_i ＝f（x_i，y_i）。插值問題就是要構造一個連續的函數F （x，y），使其在 （x_k，y_k） （後＝1，…，n） 點的函數值為Z_k，即Z_k ＝f（x_k，y_k） （k＝1 ，…，n）。

早在20世紀60年代，散亂數據的插值問題就已引起人們的注意。近50年來，已經有多種演算法被提了出來。但是，由於應用問題千差萬別，數據量大小不同，對連續性的要求也不同等等，沒有一種演算法適用於所有的場合。而且大多數演算法只能適用於具有中、小規模數據量的散亂點插值問題。大規模散亂數據 （例如，10000個點以上） 的插值問題還正在研究之中。

據散亂數據的復雜程度，其可分為單自變數、雙自變數及多自變數3種類型。下面將主要討論雙自變數散亂數據的插值問題。

（一） 插值的一般概念

插值的概念最早可追溯到 「控制論之父」 諾伯特·維納的不朽著作 《平穩時間序列的外推、插值和光滑及其工程應用》 （Wiener，1949）。隨著計算機技術的發展，插值的概念已廣泛地應用於數據處理、自動控制、數值分析、地球物理及數學地質等領域。

1. 插值方法

計算機插值方法一般可分為兩大類：擬合函數方法和加權平均方法。它們的原理都是來自手工方法。Crain （1970） 把這兩種方法得到的曲面分別稱為數學曲面和數值曲面。Alfeld ＆ Barnhill （1984） 分別稱它們為分片方法和點方法，而Cuyt （1987） 則把這兩種方法分別稱為系數問題和數值問題。

利用擬合函數方法進行插值，就是利用二元多項式來表示一個插值曲面，插值問題化為確定這個二元多項式的系數的問題。一般來說，往往可通過求解一個線性代數方程組來獲取這些系數。這個方程組的系數可由觀測數據來確定，它們代表了觀測數據的影響，而其方程的次數則表示了控制多項式擬合程度。方程次數越高，擬合的程度就越高。當這個多項式的系數確定以後，將空間某一點處的坐標代入該多項式，即可求得該點處的值。

這個方法的特點是可以制服畸變的原始數據或帶有雜訊的原始數據。所以，用一個函數進行擬合是一個光滑的過程，一些局部的細節可能消失。所得的插值曲面的復雜程度取決於多項式的次數，即所求解的線性方程的數目。

加權平均插值方法把求插值的問題化為求取觀測數據的加權平均。每個觀測數據點對應的加權系數恰恰反映了該數據點對插值點的影響大小。為了求取一個插值，必須要計算出一組加權系數。

加權平均方法的一個主要優點是，可以獲取在觀測數據點附近變數的小尺度趨勢。而利用一個適當次數的多項式是無法獲取曲面的這種局部細節的。

從原則上講，這兩種方法的差別就在於：加權平均方法強調了曲面的局部細節，而擬合函數方法則概括了曲面的整體性質。從計算時間上看，前者花費的時間比後者要多得多。

2. 插值效果評判

從理論上講，一個插值方法的效果如何應通過插值結果和客觀存在的原始曲面的比較，按以下3條標准來進行判斷：

（1） 原始曲面和插值結果之間差異的最大值為最小。

（2） 原始曲面和插值結果之差的平方和為最小。

（3） 在每個觀測數據點處，插值結果本身的數值及其1階到k階導數和原始曲面的相等。

由於原始曲面本身是未知的，所以在以上3個標准中，第一個和第二個標準是無法檢驗的，僅有第三個標準是在一定的模型假設之下可以進行檢驗。在一些實際應用中，當原始曲面可用解析函數來表達時，僅利用第三個標准來檢驗插值的效果也是可行的。然而，在地質建模中，觀測數據點往往不夠多，且還有一定的觀測誤差，不可能斷定原始曲面是否可用解析函數表達。這時，插值技術的合理性必須從直觀的幾何和人們的經驗等方面進行評價。

利用計算機進行插值所遇到的困難，主要來自觀測數據點數目不足和觀測誤差。如果觀測數據充分多且精確，那麼幾乎所有的插值方法都會給出良好的效果。另一方面，對於圓形或狹長的隆起，凹陷和鞍點等變數的空間變化幾何特徵，在數據點分布較稀的情況下，用任何插值方法都是難以推斷出它們的存在的。所以，插值方法需要考慮曲面的局部斜率的影響。

下面主要介紹加權平均方法和擬合函數方法中最常用的插值演算法。這些演算法能解決大部分油氣藏建模問題。

（二） 與距離成反比的加權法

距離成反比加權插值方法是基於如下的模型：每個數據點都有局部影響，這個影響隨著數據點和插值點距離的增加而減弱，且在一定的范圍以外，可以忽略不計；這個影響是以該數據點為中心；而在任一點處的插值恰是各數據點影響之和。

1. 與距離成反比加權插值公式

這一方法首先是由氣象學及地質學工作者提出來的，後來由於D. Shepard的工作被稱為Shepard方法。其基本思想是將插值函數F （x，y）定義為各數據點函數值f i的加權平均，即：

油氣田開發地質學

在 （x_k，y_k） 點處函數值可寫成：

油氣田開發地質學

式中： 表示由第i個（x_i，y_i）點到插值點（x_k，y_k）的距離；W_i（x_k，y_k）——權函數；μ——功率因素，通過改變值來調整權函數與距離的關系，與距離成反比加權和與距離成平方反比加權分別是μ＝1和μ＝2時的特殊情形。

與距離成反比加權插值方法是最早使用的計算機插值方法，至今仍被廣泛地應用著。在大多數商業性的等值線圖繪制軟體包中被用來形成網格化數據。這種方法較為直觀：一個數據點對於插值點的影響模型化為與這兩點之間的距離成反比。

2. 與距離成反比加權插值改進

距離成反比加權演算法中的功率因素μ應該取為μ≥0，否則表明距離越遠的點作用越突出，這違背了普通常識。考慮以下最極端的情況是：

油氣田開發地質學

假設有n個數據點，一個插值點 （x_k，y_k） 位於第i個數據點附近，相應的權函數可寫成：

油氣田開發地質學

式中：d_j （x_k，y_k），d_i （x_k，y_k）——插值點 （x_k，y_k）到各數據點的距離。

當插值點 （x_k，y_k） 和第i個數據點很靠近時，可以認為其他數據點對W_D的影響是一個常數，即C為常數。當該插值點和第i個數據點的距離趨於零時，對於不同的μ，W_D會有不同的性質。首先看W_D對D的導數，有：

油氣田開發地質學

然後再有：

油氣田開發地質學

可見，當μ＝1時，隨著D趨於0，W′_D近似為不隨D變化的常數，這可用圖6-8中左端的圖形來表示。

當μ>1時，隨著D趨於零，W′_D趨於零。這說明在該數據點附近，加權系數的變化為零，即可用圖6-8中間的圖形來表示。

當μ<1時，隨著D趨於零，W′_D趨於無窮大。這說明加權系數在該數據點附近還有一個尖點，如圖6-8右端的圖形所示。

圖6-8 與距離反比加權的權數隨參數μ的變化

（1）功率因素μ越小時，近距離點和遠距離點的作用越接近，生成的平面網格數據越平滑。隨著μ增大，平面網格數據的光滑性越差。同時，μ直接影響網格數據的極值和均值，μ越小網格數據的均值越接近原始數據的均值，但極值相差越大。因此，當要求插值結果盡可能接近原始數據的均值時，功率因素不能選擇過大。例如對於開發早期的油氣藏，因為僅有少數探井控制，此時網格化得到的各類物性參數應該在總體上符合井點的統計結果，均值是比極值更有價值的參數，因而通常將功率因素取為1。相反，μ取值越大，網格數據越能恢復原始數據的極值，但也容易使均值誤差增大。原因是當μ取較大的值時，近距離點的作用越突出，原始數據點分布的不均勻性使得部分點在網格節點上發揮了更大的作用，而另外一些點的作用則受到屏蔽。因此，當要求突出數據的局部特徵，體現儲層的非均質性特徵時，功率因素應選擇得大些，一般取為2。

（2） 利用與距離成反比加權法進行插值時，當增加、刪除或改變一個點時，權函數W_i （x_k，y_k） 均需重新計算，因而該方法是一個全局插值演算法。

為了克服Shepard方法的上述缺陷，Franke及Nielson提出了MQS （Modified QuadraticShepard） 方法，它仍然是一個與距離成反比的加權方法。對它的改進如下：

插值點 （x_k，y_k） 到已知數據點的距離d_i作適當修改，使其只能在局部范圍內起作用，以改變Shepard方法的全局插值性質。這時重新定義距離函數：

油氣田開發地質學

式中：r_w為一個常數。而

油氣田開發地質學

因此，當 （x_k，y_k） 點與某一點的距離大於r_w時，權值就為零。

（3） 與距離成反比加權法的權函數W_i（x_k，y_k）始終滿足W_i（x_k，y_k）≤1/n，因此插值結果不會大於或小於原始數據的最大值、最小值。當已知數據點過少時 （這種情況在早期地質研究中是最為常見），使用具有外推能力的曲面樣條或趨勢面分析，得到的結果往往背離實際。其原因是這兩種方法在遠點不具有控制能力，它將沿趨勢無限發展下去。特別是在僅有少數井資料可用的情況下，與距離成反比加權法應是優先選擇的方法。

但是，隱含在原始數據中的尖峰會被淹沒而無法顯示出來。因為距離成反比加權進行插值時，每個數據點所發揮的作用基本上是中心對稱的，因此對山脊和山谷等非各向同性的幾何形狀的顯示不利。為了克服這種狀況，需要考慮變數的局部變化趨勢，為此需要對梯度進行估計。

當已知數據點用與距離成反比加權方法形成數據插值曲面時，該曲面被稱數據曲面。與距離成反比加權方法也可用於各個數據點處的切平面，把任一點處的插值值取成各切平面在該點處取值的一個加權平均，所形成的曲面稱為與距離成反比加權梯度插值曲面，簡稱梯度曲面。如果在數據點以外的一個點是變數的局部高點，那麼梯度曲面在該點的值容易大於變數的真實值，即呈現 「過估計」 的狀態。如果數據曲面在該點的值小於變數的真實值，則呈現 「欠估計」 的狀態。因此，數據曲面可以通過和梯度曲面的相互結合來克服本身的缺陷。

可以用數據曲面和梯度曲面之差乘以一個系數作為一個修正量，對數據曲面進行修正，這樣，可以用下述的曲面來代替單純的數據曲面：

油氣田開發地質學

式中：L （x，y） ——距離反比加權插值；W_i——和第i個數據點的距離反比加權系數；τ_i——曲面在該點處的粗糙度指數；S_i（x，y）——第i個數據點 （x_i，y_i） 處的切平面在點（x，y） 處的值；H（W_i，τ_i）——混合函數。

如此得到的曲面稱為混合曲面。它通過所有的數據點，具有連續的坡度，其變化在空間的分布更均勻。數據曲面和梯度曲面是混合曲面的兩種極端情況。由於數據曲面和梯度曲面之差在各數據點處為零，還因為混合函數的變化范圍為0～1，且當混合函數等於0或1時，其一階導數為零，故混合曲面和梯度曲面相切於各數據點處，且其高階導數在各數據點處亦為零。

（4） 與距離成反比加權法僅考慮了插值點與數據點之間距離的影響，沒有考慮到各數據點之間的關系，物性參數分布的趨勢性沒有得到充分的體現。為此，Franke及Nielson進行了改進。

用節點函數Q_i （x，y） 代替f_i，Q_i （x，y） 是一個插值於 （x_i，y_i） 點的二次多項式，即有Q_i （x_i，y_i） ＝f，i＝1，…，n。Q_i可由下式表示：

Q_i（x，y）＝f_i＋a₁（x-x_i）＋a₂（y-y_i）＋a₃（x-x_i）2＋a₄（x-x_i）（y-y_i）＋a₅（y-y_i）²

式中：a₁，a₂，…，a₅是按下式最小二乘法得出的優化解：

油氣田開發地質學

式中：f_i，f_j分別為 （x_i，y_i）和 （x_j，y_j）點的函數值，而ρ_j可按下式選取：

油氣田開發地質學

其中r_q為一常數，而

油氣田開發地質學

求出Q_i （x，y） 後，插值函數可表示為：

油氣田開發地質學

上述方法消除了Shepard方法中的一些缺陷，因而在散亂點插值中得到廣泛的應用。但是，為了求得Q_i （x，y） （i＝1，…，n），需要多次求解線性方程組，計算量大，因此，一般只用於中、小規模散亂點的插值運算。

（三） 多項式趨勢面法

由計算機產生的曲面一般不會總是和原始的觀測數據一致。如果兩者的差別在給定的尺度之下不是很明顯，那麼產生的曲面可被認為是插值曲面，否則就被認為是近似曲面。如果觀測數據含有明顯的觀測誤差，近似曲面就顯得更合理。這時，和插值曲面相比，近似曲面由於數據的各種誤差所產生的擾動不太容易看得清，但是近似曲面空間變化的一些主要性質還是能清晰地被體現出來的。

近似曲面和每個數據點之間的差稱為殘差，可視為每一個數據點上的一種誤差表示。然而，計算出來的這種殘差是意味著對未知的觀測誤差的一種度量，還是意味著一種允許的插值誤差，或者意味兩者都是，這要依變數的空間性質和觀測數據的獲取方法而定。

確定近似曲面的方法可分為3種。第一種方法是以殘余的平方和最小為條件，確定多項式的系數，以獲取曲面。第二種方法是利用觀測數據誤差的附加信息，並滿足最小曲率的原則以確定曲面。最後一種方法是利用觀測誤差和插值誤差的附加信息，以滿足最小平方差或最小曲率為條件確定曲面。以下主要討論多項式構造趨勢面法。

多項式構造趨勢面是目前最常用的方法，一次多項式表示的趨勢面是空間的一個平面，二次趨勢面是拋物面，橢球面或雙曲面，三次及三次以上的趨勢面是形態復雜的空間曲面，隨著趨勢面的次數增高，曲面的形態就越復雜。

如果有一組總共n個觀測數據，其觀測點的平面坐標為 （x_i，y_i），地質變數的觀測值為f_i（x_i，y_i）。對於這組觀測數據的多項式趨勢面方程表示成如下形式：

油氣田開發地質學

式中： ——第i個觀測點的趨勢值；a₁，a₂，…，a₅——待定系數，它們的個數m與所選用趨勢面方程的多項式次數n存在下列關系式：

m＝［n（n＋3）＋2］/2

為使趨勢面最大限度地逼近原始觀測數據，可採用最小二乘法使每個觀測點的觀測值與趨勢值之差 （殘差） 的平方和最小，即：

油氣田開發地質學

得到需求解的m階正規方程組。當系數矩陣滿秩時，趨勢面方程也就被唯一確定了。

由於趨勢面分析不具有過點性，使得局部井點上誤差可能很大。同時參數場在三維空間中的分布過於復雜，無論從理論上還是實驗中都無法確證某類參數場能較好地符合某確定次數的曲面，多項式次數過高，會導致趨勢面發生頻繁振動，多項式次數過低，得到的趨勢面又過於光滑，喪失許多細節。再者，趨勢面方程在外推過程中容易使參數場發生畸變，產生無意義的結果。因而現代地質建模研究中已經很少將趨勢面分析單獨作為插值方法使用。但對於下列兩種情況趨勢面分析仍然能達到較為理想的效果，一是對小范圍內具有顯著趨勢性分布的數據點；二是數據點分布過於密集，而且可能存在若干異常數據點時，趨勢面分析會自動削弱異常數據點的影響。為提高趨勢面分析的精度，可採用殘差來校正趨勢面分析的結果。具體過程如下：

（1） 由趨勢面分析得到任一網格節點 （x_k，y_k）處的趨勢值

油氣田開發地質學

（2）計算n個已知點 （x_i，y_i）處的殘差△f_i（x_i，y_i）＝f_i（x_i，y_i）

油氣田開發地質學

（3） 調用某種插值方法將殘差分配到每個網格節點上，對網格節點 （x_k，y_k） 有△f_k（x_k，y_k）；

（4） 網格節點 （x_k，y_k）經校正後的最終結果為

油氣田開發地質學

特別地，如果殘差分配的插值演算法也是趨勢面分析，就形成所謂的多級趨勢面分析。此時，網格節點上的值是多次趨勢面分析的結果，多級趨勢面分析在不斷減少觀測點插值誤差的同時，整張參數場曲面仍然保持其連續性和光滑性，原因是該演算法同樣符合線性迭加原理。

（四） 徑向基函數插值法

徑向基函數的名字來源於這樣一種情況，即基函數是由單個變數的函數構成的。一個點 （x，y） 的這種基函數的形式往往是h_k（x，y）＝h（d_k），這里的d_k表示由點 （x，y） 至第k個數據點的距離。一般說來，這種方法不具有多項式精度，但只要稍加改進，即可獲得具有多項式精度的插值公式：

油氣田開發地質學

式中：q_k（x，y）是一個多項式基，其階次小於m。

上式中的系數a_k和b_k應滿足下面的聯立方程組：

油氣田開發地質學

油氣田開發地質學

第一式中的n個方程式滿足了插值要求，而第二式中的m個方程式則保證了多項式精度。兩式中共有m＋n個未知數，同時存在m＋n個方程式，聯立求解，即可得出待定系數。

下面，介紹兩種主要的徑向基函數插值法。

1. Multiquadric方法

Multiquadric方法是由R. L. Hardy在1971年提出來的。它是最早被提出並且應用得最為成功的一種徑向基函數插值法。它採用的插值函數，即 （x，y）處的值F （x，y）：

油氣田開發地質學

式中： 為基函數；e_i——非負常數；a_i——加權系數，滿足如下方程組：

MV_a＝V_z

式中：V_a=(a₁,a₂,…,a_n)^T,V_z=[f(x₁,y₁),f(x₂,y₂),…,f(x_n,y_n)]^T,f(x_i,y_i)是(x_i,y_i)處的數據點的值。

油氣田開發地質學

由於M和V_z都不依賴於插值點的坐標 （x，y），所以a_i也不依賴於 （x，y）。然而，基函數C（X-X_i）則是以Xi為參數的 （x，y）的函數。所以說，插值曲面F（x，y）是n個基函數C（X-X_i）所構成的n個空間曲面配置而成的。

Arther提出了如下形式的基函數：C（d）＝1-d²/e²，其中d是數據點到插值點之間距離，而e則是一個常數。顯然，基函數C（d）是d的一個衰減函數。當d＝0時，C（d）取得最大值1，而當d≤e時，有0≤C（d）≤1。這時，基函數呈現為橢圓拋物面，而加權系數a_i（i＝1，2，…，n）所滿足的線性代數方程組的矩陣M應作相應的改動，其對角線元素應改成1。Hardy（1971）引入如下的基函數：C（d）＝（d²＋e²）^1/2，該基函數呈現為橢圓雙曲面。Hardy還建議將e²取成0.815乘以數據點間距離的平均值。

2. 薄板樣條法

樣條 （Spline）本來是繪圖員用來繪制光滑曲線的工具，是一種用木材或金屬等彈性材料做成的細條。在繪圖時，沿著通過圖紙各已知點的樣條，便可繪出一條光滑曲線。數學上所說的樣條 （多項式樣條） 實質上是分段多項式曲線的光滑連接。當函數為分段的m次多項式，在分段點上有直至m-1階連續導數，那麼該函數則稱為m次樣條函數，簡稱為樣條。一般來說，研究和應用得比較多的是三次樣條。零次和一次樣條函數分別是台階狀函數和折線狀函數。以上所述的是關於一維樣條函數，對於二維樣條函數也可作為類似的考慮。

樣條函數的主要功能是進行插值，其主要優點在於，能在插值多項式的次數盡可能低的條件下，使插值曲線或插值曲面取得較高的光滑度，且只需要利用函數本身的值，而不需要提供函數的各階導數的值。

三次樣條曲麵包含有三種不同的類型：雙三次樣條、偽三次樣條及薄片樣條。這些樣條曲面以m和s為其兩個參數，使得希氏空間Hs中元素的m階導數的范數所構成的一個泛函達到最小。此外，這個泛函具有旋轉不變性。

對於薄片樣條曲面，m＝2，s＝0。這一方法是由R.L. Harder及R. N. Desmarais在1972年提出來的，後來由J. Duchon及J. Meinguet等人予以發展。薄板樣條法得名於如下事實，即用此方法求出的散亂點的插值函數使下面這一泛函表達式具有最小值：

油氣田開發地質學

在這里，I（F）表示受限於插值點的無限彈性薄板的彎曲能量。因此，這一方法的實質從力學觀點看是使插值函數所代表的彈性薄板受限於插值點，並且具有最小的彎曲能量。這是一個泛函求極值的問題。這一變分問題的解即為我們所需要的插值函數，具有徑向基函數插值法的一般形式。

R.L. Harder及R. N. Desmarais提出解析形式如下：

油氣田開發地質學

且有 其中t（x，y）和t_i（x_i，y_i）是二維空間中的點，而f_i是t_i處的觀測值。還有，K（t_i，t） 這里，r_i代表點t和t_i之間的距離，

由解析表達式及其約束條件，可給出用以確定系數的線性代數方程組：

油氣田開發地質學

其中，K=[k_ij]_n×n,k_ij=K(t_i,t_j),k_ii=0,F^T[f₁,f₂,…,f_n],α^T=[b,a₁,a₂],A^T［λ₁，λ₂，…，λ_n］，且有：

油氣田開發地質學

求解上述n＋3階方程組則得到待定系數b，a₁，a₂，λ_i，然後可插值出平面任一位置的函數值F（x，y）。

上述方程與Enriguez等人給出的方程相類似，其微小的差異就是基函數中的自然對數（In） 變成了常用對數 （lg）。

在構造薄片樣條曲面的過程中，Franke （1982）提出了以r²lgr作為基函數，Sandwell（1987）則提出了以雙調和格林函數r²（lgr-1）作為基函數。另外，Ayeni （1979）也對不同的基函數進行了討論。

使用平面上n個已知點進行曲面樣條插值時，實質上是求解一個n＋3階線性方程組以確定n＋3個系數。為了保證解的存在性和唯一性，系數矩陣應該是滿秩的。對下列3種情況必須避免：

（1） 在給定的n個已知點中存在著距離過近的兩點 （x_i，y_i） 與 （x_j，y_j） 極端的情形是同一個數據點的重復輸入，此時系數矩陣的第i行與第j行對應各元素非常接近，導致線性方程組的系數矩陣是奇異的。因此在數據預處理過程中必須消除沉余數據。

（2） 給定的已知數據點過少，此時系數矩陣的後3行線性相關，矩陣是不滿秩的。

（3） n個已知點數據分布在一條直線上，顯然由這樣的n個點不能唯一決定一張曲面，系數矩陣表現為不滿秩。

針對情況 （1），通常在做曲面樣條插值時，首先對數據進行預處理，通過給定一個適當的距離下限r_min來濾掉那些相距過近的點，研究發現r_min＝（△x＋△y）/8是一個較合理選擇 （△x和△y分別表示x和y方向的步長）。情況 （2） 和 （3） 實際上意味著不能進行曲面樣條插值，除非通過數據均整來改變數據分布狀態。

曲面樣條插值方法是一種嚴格的過點插值法，即由生成的樣條曲面必定通過給定的n個已知數據點，這樣井點數據的控製作用自然得到體現。同時，曲面樣條方法充分考慮了數據間的相對位置，其插值精度很高，在外推過程中，總是沿數據點的分布趨勢外推，因此曲面樣條法是具有一定外推能力的插值方法。同時曲面樣條方程得到的是一張連續光滑的曲面。

曲面樣條插值方法特別適合於地層層面的生成和地層厚度的插值。如果從曲面樣條法嚴格的過點性、良好的光滑性及外推性看，適用於那些光滑、趨勢性明顯、變化連續的儲層物性參數諸如油氣飽和度的插值。

曲面樣條法插值的精度很大程度上取決於數據點分布的均勻程度，稀疏區域主要由鄰近區域的外推得到，其插值結果可能偏差較大。同時，曲面樣條法的嚴格過點性使得它不能分別對待不同精度點的數據，因此它不具備數據的校正能力。

3. 回歸直線方程的計算方法

要確定回歸直線方程①，只要確定a與回歸系數b。回歸直線的求法通常是最小二乘法：離差作為表示xi對應的回歸直線縱坐標y與觀察值yi的差，其幾何意義可用點與其在回歸直線豎直方向上的投影間的距離來描述。數學表達:Yi-y^=Yi-a-bXi.總離差不能用n個離差之和來表示，通常是用離差的平方和即(Yi-a-bXi)^2計算。即作為總離差，並使之達到最小，這樣回歸直線就是所有直線中除去最小值的那一條。這種使「離差平方和最小」的方法，叫做最小二乘法。用最小二乘法求回歸直線方程中的a,b有圖一和圖二所示的公式進行參考。其中，

(3)廣義極小殘差法演算法擴展閱讀

回歸直線方程指在一組具有相關關系的變數的數據（x與Y）間，一條最好地反映x與y之間的關系直線。

離差作為表示Xi對應的回歸直線縱坐標y與觀察值Yi的差，其幾何意義可用點與其在回歸直線豎直方向上的投影間的距離來描述。數學表達:Yi-y^=Yi-a-bXi.

總離差不能用n個離差之和來表示，通常是用離差的平方和，即(Yi-a-bXi)^2計算。

4. 《R語言實戰》自學筆記71-主成分和因子分析

主成分分析
主成分分析（(Principal Component Analysis，PCA）是一種數據降維技巧，它能將大量相關變數轉化為一組很少的不相關變數，這些無關變數稱為主成分（原來變數的線性組合）。整體思想就是化繁為簡，抓住問題關鍵，也就是降維思想。
主成分分析法是通過恰當的數學變換，使新變數——主成分成為原變數的線性組合，並選取少數幾個在變差總信息量中比例較大的主成分來分析事物的一種方法。主成分在變差信息量中的比例越大，它在綜合評價中的作用就越大。

因子分析
探索性因子分析法（Exploratory Factor Analysis，EFA）是一系列用來發現一組變數的潛在結構的方法。它通過尋找一組更小的、潛在的或隱藏的結構來解釋已觀測到的、顯式的變數間的關系。

PCA與EFA模型間的區別
參見圖14-1。主成分（PC1和PC2）是觀測變數（X1到X5）的線性組合。形成線性組合的權重都是通過最大化各主成分所解釋的方差來獲得，同時還要保證個主成分間不相關。相反，因子（F1和F2）被當做是觀測變數的結構基礎或「原因」，而不是它們的線性組合。

R的基礎安裝包提供了PCA和EFA的函數，分別為princomp()和factanal()。
最常見的分析步驟
(1)數據預處理。PCA和EFA都根據觀測變數間的相關性來推導結果。用戶可以輸入原始數據矩陣或者相關系數矩陣到principal()和fa()函數中。若輸入初始數據，相關系數矩陣將會被自動計算，在計算前請確保數據中沒有缺失伍世拿值。
(2)選擇因子模型。判斷是PCA（數據降維）還是EFA（發現潛在結構）更符合你的研究目標。如果選擇EFA方法，你還需要選擇一種估計因子模型的方法（如最大似然估計）。
(3)判斷要選擇的主成分/因子數目。
(4)選擇主成分/因子。
(5)旋轉主成分/因子。
(6)解釋結果。
(7)計算主成分或因子得分。

PCA的目標是用一組較少的不相關變數代替大量相關變數，同時盡可能保留初始變數的信息，這些推導所得的變數稱為主成分，它們是觀測變數的線性組合。如第一主成分為：

它是k個觀測變數的加權組合，對初始變數集的方差解釋性最大。第二主成分也是初始變數的線性組合，對方差的解釋性排第二，同時與第一主成分正交（不相關）。後面每一個主成分都最大化它對方差的解釋程度，同時與之前所有的主成分都正交。理論上來說，你可以選取與變數數相同的主成分，但從實用的角度來看，我們都希望能用較少的腔搭主成分來近似全變數集。

主成分與原始變數之間的關系
（1）主成分保留了原始變數絕大多數信息。
（2）主成分的個數大大少於原始變數的數目。
（3）各個主成分之間互不相關。
（4）每個主成分都是原始變數的線性組合。

數據集USJudgeRatings包含了律師對美國高等法院法官的評分。數據框包含43個返州觀測，12個變數。

用來判斷PCA中需要多少個主成分的准則：
根據先驗經驗和理論知識判斷主成分數；
根據要解釋變數方差的積累值的閾值來判斷需要的主成分數；
通過檢查變數間k × k的相關系數矩陣來判斷保留的主成分數。
最常見的是基於特徵值的方法。每個主成分都與相關系數矩陣的特徵值相關聯，第一主成分與最大的特徵值相關聯，第二主成分與第二大的特徵值相關聯，依此類推。
Kaiser-Harris准則建議保留特徵值大於1的主成分，特徵值小於1的成分所解釋的方差比包含在單個變數中的方差更少。Cattell碎石檢驗則繪制了特徵值與主成分數的圖形。這類圖形可以清晰地展示圖形彎曲狀況，在圖形變化最大處之上的主成分都可保留。最後，你還可以進行模擬，依據與初始矩陣相同大小的隨機數據矩陣來判斷要提取的特徵值。若基於真實數據的某個特徵值大於一組隨機數據矩陣相應的平均特徵值，那麼該主成分可以保留。該方法稱作平行分析。

圖形解讀：線段和x符號組成的圖（藍色線）：特徵值曲線；
紅色虛線：根據100個隨機數據矩陣推導出來的平均特徵值曲線；
綠色實線：特徵值准則線（即：y=1的水平線）
判別標准：特徵值大於平均特徵值，且大於y=1的特徵值准則線，被認為是可保留的主成分。根據判別標准，保留1個主成分即可。

fa.parallel函數學習
fa.parallel(data,n.obs=,fa=」pc」/」both」,n.iter=100,show.legend=T/F)
data：原始數據數據框；
n.obs：當data是相關系數矩陣時，給出原始數據（非原始變數）個數，data是原始數據矩陣時忽略此參數；
fa：「pc」為僅計算主成分，「fa」為因子分析，「both」為計算主成分及因子；
n.iter：模擬平行分析次數；
show.legend：顯示圖例。

principal(r, nfactors = , rotate = , scores = )

r：相關系數矩陣或原始數據矩陣；
nfactors：設定主成分數（默認為1）；
rotate：指定旋轉的方法，默認最大方差旋轉（varimax）。
scores：設定是否需要計算主成分得分（默認不需要）。

PC1欄包含了成分載荷，指觀測變數與主成分的相關系數。如果提取不止一個主成分，那麼還將會有PC2、PC3等欄。成分載荷（component loadings）可用來解釋主成分的含義，解釋主成分與各變數的相關程度。
h2欄為成分公因子方差，即主成分對每個變數的方差解釋度。
u2欄為成分唯一性，即方差無法被主成分解釋的部分（1-h2）。
SS loadings包含了與主成分相關聯的特徵值，其含義是與特定主成分相關聯的標准化後的方差值，即可以通過它來看90%的方差可以被多少個成分解釋，從而選出主成分（即可使用nfactors=原始變數個數來把所有特徵值查出，當然也可以直接通過eigen函數對它的相關矩陣進行查特徵值）。
Proportion Var表示每個主成分對整個數據集的解釋程度。
Cumulative Var表示各主成分解釋程度之和。
Proportion Explained及Cumulative Proportion分別為按現有總解釋方差百分比劃分主成分及其累積百分比。

結果解讀：第一主成分（PC1）與每個變數都高度相關，也就是說，它是一個可用來進行一般性評價的維度。ORAL變數99.1%的方差都可以被PC1來解釋，僅僅有0.91%的方差不能被PC1解釋。第一主成分解釋了11個變數92%的方差。

結果解讀：通過碎石圖可以判定選擇的主成分個數為2個。

結果解讀：從結果Proportion Var： 0.58和0.22可以判定，第一主成分解釋了身體測量指標58%的方差，而第二主成分解釋了22%，兩者總共解釋了81%的方差。對於高度變數，兩者則共解釋了其88%的方差。

旋轉是一系列將成分載荷陣變得更容易解釋的數學方法，它們盡可能地對成分去噪。旋轉方法有兩種：使選擇的成分保持不相關（正交旋轉），和讓它們變得相關（斜交旋轉）。旋轉方法也會依據去噪定義的不同而不同。最流行的正交旋轉是方差極大旋轉，它試圖對載荷陣的列進行去噪，使得每個成分只是由一組有限的變數來解釋（即載荷陣每列只有少數幾個很大的載荷，其他都是很小的載荷）。 結果列表中列的名字都從PC變成了RC，以表示成分被旋轉。

當scores = TRUE時，主成分得分存儲在principal()函數返回對象的scores元素中。

如果你的目標是尋求可解釋觀測變數的潛在隱含變數，可使用因子分析。
EFA的目標是通過發掘隱藏在數據下的一組較少的、更為基本的無法觀測的變數，來解釋一
組可觀測變數的相關性。這些虛擬的、無法觀測的變數稱作因子。（每個因子被認為可解釋多個
觀測變數間共有的方差，因此准確來說，它們應該稱作公共因子。）

其中 是第i個可觀測變數（i = 1…k）， 是公共因子（j = 1…p），並且p<k。 是 變數獨有的部分（無法被公共因子解釋）。 可認為是每個因子對復合而成的可觀測變數的貢獻值。

碎石檢驗的前兩個特徵值（三角形）都在拐角處之上，並且大於基於100次模擬數據矩陣的特徵值均值。對於EFA，Kaiser-Harris准則的特徵值數大於0，而不是1。
結果解讀：PCA結果建議提取一個或者兩個成分，EFA建議提取兩個因子。

fa(r, nfactors=, n.obs=, rotate=, scores=, fm=)
 r是相關系數矩陣或者原始數據矩陣；
 nfactors設定提取的因子數（默認為1）；
 n.obs是觀測數（輸入相關系數矩陣時需要填寫）；
 rotate設定旋轉的方法（默認互變異數最小法）；
 scores設定是否計算因子得分（默認不計算）；
 fm設定因子化方法（默認極小殘差法）。
與PCA不同，提取公共因子的方法很多，包括最大似然法（ml）、主軸迭代法（pa）、加權最小二乘法（wls）、廣義加權最小二乘法（gls）和最小殘差法（minres）。統計學家青睞使用最大似然法，因為它有良好的統計性質。

結果解讀：兩個因子的Proportion Var分別為0.46和0.14，兩個因子解釋了六個心理學測試60%的方差。

結果解讀：閱讀和詞彙在第一因子上載荷較大，畫圖、積木圖案和迷宮在第二因子上載荷較大，非語言的普通智力測量在兩個因子上載荷較為平均，這表明存在一個語言智力因子和一個非語言智力因子。

正交旋轉和斜交旋轉的不同之處。
對於正交旋轉，因子分析的重點在於因子結構矩陣（變數與因子的相關系數），而對於斜交旋轉，因子分析會考慮三個矩陣：因子結構矩陣、因子模式矩陣和因子關聯矩陣。
因子模式矩陣即標准化的回歸系數矩陣。它列出了因子預測變數的權重。因子關聯矩陣即因子相關系數矩陣。

圖形解讀：詞彙和閱讀在第一個因子（PA1）上載荷較大，而積木圖案、畫圖和迷宮在第二個因子（PA2）上載荷較大。普通智力測驗在兩個因子上較為平均。

與可精確計算的主成分得分不同，因子得分只是估計得到的。它的估計方法有多種，fa()函數使用的是回歸方法。

R包含了其他許多對因子分析非常有用的軟體包。FactoMineR包不僅提供了PCA和EFA方法，還包含潛變數模型。它有許多此處我們並沒考慮的參數選項，比如數值型變數和類別型變數的使用方法。FAiR包使用遺傳演算法來估計因子分析模型，它增強了模型參數估計能力，能夠處理不等式的約束條件，GPArotation包則提供了許多因子旋轉方法。最後，還有nFactors包，它提供了用來判斷因子數目的許多復雜方法。

主成分分析

1.數據導入
數據結構：對10株玉米進行了生物學性狀考察，考察指標有株高，穗位，莖粗，穗長，禿頂，穗粗，穗行數，行粒數。

結果解讀：選擇2個主成分即可保留樣本大量信息。

3.提取主成分

結果解讀：主成分1可解釋44%的方差，主成分2解釋了26%的方差，合計解釋了70%的方差。

4.獲取主成分得分

5.主成分方程

PC1 = 0.27 株高 - 0.04 穗位 + 0.29 莖粗 - 0.01 穗長 - 0.21 禿頂 - 0.13 穗粗 + 0.16 穗行數 + 0.24 行粒數

PC2 = -0.01 株高 + 0.36 穗位 - 0.10 莖粗 + 0.41 穗長 - 0.08 禿頂 + 0.43 穗粗 - 0.15 穗行數 + 0.01 行粒數

圖形解讀：此圖反映了變數與主成分的關系，三個藍點對應的RC2值較高，點上的標號2，4，6對應變數名穗位，穗長，穗粗，說明第2主成分主要解釋了這些變數，與這些變數相關性強；黑點分別對應株高，莖粗，穗行數，行粒數，說明第一主成分與這些變數相關性強，第一主成分主要解釋的也是這些變數，而5號點禿頂對於兩個主成分均沒有顯示好的相關性。

因子分析

圖解：可以看到需要提取4個因子。

2.提取因子

結果解讀：因子1到4解釋了80%的方差。

3.獲取因子得分

圖解：可以看出，因子1和因子2的相關系數為0.4，行粒數，株高，莖粗，禿頂在因子1的載荷較大，穗長，穗位在因子2上的載荷較大；因子3隻有穗行數相關，因子4隻有穗粗相關。

參考資料：