R

什麼是 R 中計算的多元正交多項式?

  • December 24, 2016

單變量點集中的正交多項式是以點積和成對相關性為零的方式在該點上產生值的多項式。R 可以產生具有函數poly的正交多項式。

相同的函數有一個變體多項式,它在多元點集上產生正交多項式。無論如何,在具有成對零相關性的意義上,得到的多項式不是正交的。實際上,由於一階多項式應該只是原始變量,因此除非原始變量不相關,否則一階多項式不會是正交的。

然後,我的問題是:

  • 由 R 中的 polym 計算的多元正交多項式是什麼?它們只是單變量正交多項式的乘積嗎?它們是用來做什麼的?
  • 可以存在真正的多元正交多項式嗎?有沒有一種簡單的方法來生產它們?在 R? 它們實際上是否用於回歸?

更新

為了回應 Superpronker 的評論,我舉了一個例子來說明我對不相關多項式的意思:

> x<-rnorm(10000)
> cor(cbind(poly(x,degree=3)))
             1             2             3
1  1.000000e+00 -6.809725e-17  2.253577e-18
2 -6.809725e-17  1.000000e+00 -2.765115e-17
3  2.253577e-18 -2.765115e-17  1.000000e+00

Poly 函數返回以點 x 計算的正交多項式(這裡每個多項式有 10,000 個點)。不同多項式上的值之間的相關性為零(有一些數值誤差)。

當使用多元多項式時,相關性不為零:

> x<-rnorm(1000)
> y<-rnorm(1000)
> cor(cbind(polym(x,y,degree=2)))
             1.0           2.0           0.1         1.1           0.2
1.0  1.000000e+00  2.351107e-17  2.803716e-02 -0.02838553  3.802363e-02
2.0  2.351107e-17  1.000000e+00 -1.899282e-02  0.10336693 -8.205039e-04
0.1  2.803716e-02 -1.899282e-02  1.000000e+00  0.05426440  5.974827e-17
1.1 -2.838553e-02  1.033669e-01  5.426440e-02  1.00000000  8.415630e-02
0.2  3.802363e-02 -8.205039e-04  5.974827e-17  0.08415630  1.000000e+00

因此,我不明白這些雙變量多項式在什麼意義上是正交的。

更新 2

我想澄清回歸中使用的“正交多項式”的含義,因為在連接區間中應用正交多項式的想法時,這種上下文可能會產生某種誤導 - 正如最後一個 Superpronker 的評論。

我使用 R第 101 和 102 頁引用 Julian J. Faraway 的Practical Regression and Anova :

正交多項式通過定義 z1=a1+b1x

z2=a2+b2x+c2x2
z3=a3+b3x+c3x2+d3x3
等等,其中選擇係數 a、b、c… 使得 zTi·zj=0 什麼時候 ij . z 稱為正交多項式。

通過輕微的語言濫用,作者在這裡使用 zi 對於多項式(作為函數)和多項式在集合點中的值的向量 x . 或者它甚至根本不是語言濫用,因為從本書的開頭 x 一直是預測器(例如,預測器採用的一組值)。

正交多項式的這種含義實際上與區間上的正交多項式沒有區別。我們可以以通常的方式(使用積分)在任何具有任何測量函數的可測量集合上定義正交多項式。這裡我們有一個有限集( x ) 並且我們使用的是點積而不是積分,但是如果我們將測量函數作為有限集點中的狄拉克增量,那仍然是正交多項式。

並且關於相關性:正交向量的點積 Rn (作為有限集上的正交向量的圖像)。如果兩個向量的點積為零,則協方差為零,如果協方差為零,則相關性為零。在線性模型的上下文中,關聯“正交”和“不相關”非常有用,如“正交實驗設計”。

讓我們探索發生了什麼。我相信您已經知道以下大部分材料,但是為了建立符號和定義並使想法清晰,我將在回答問題之前介紹多項式回歸的基礎知識。如果您願意,請跳到R這篇文章大約三分之二的標題“做什麼”,然後跳回您可能需要的任何定義。

那個設定

我們正在考慮一個 n×k 模型矩陣 X 某種回歸中的潛在解釋變量。這意味著我們正在考慮的列 X 作為存在 n -向量 X1,X2,,Xk 我們將形成它們的線性組合, β1X1+β2X2++βkXk, 預測或估計反應。

有時可以通過引入通過將不同列相乘而創建的附加列來改進回歸 X 彼此,一個係數一個係數。這樣的產品稱為“單項式”,可以寫成

Xd11Xd22Xdkk

每一個“權力” di 為零或更大,表示每個次數 X1 出現在產品中。請注意 X0 是一個 n - 常數係數向量 ( 1 ) 和 X1=X 本身。因此,單項式(作為向量)生成一個向量空間,其中包括 X. 它可能是一個更大的向量空間的可能性為這個過程提供了更大的空間來模擬線性組合的響應。

我們打算替換原來的模型矩陣 X 通過單項式的集合線性組合。當這些單項式中至少一個的次數超過 1, 這稱為多項式回歸。

多項式的分級

單項式的次數是其冪的總和, d1+d2++dk. 單項式的線性組合(“多項式”)的次數是具有非零係數的單項項中最大的次數。 **度數具有內在意義,**因為當你改變原始向量空間的基礎時,每個向量 Xi 由所有向量的線性組合新表示;單項式 Xd11Xd22Xdkk 從而成為相同次數的多項式;因此任何多項式的次數都是不變的。

度數為這個多項式代數提供了一個自然的“分級”:由單項式的所有線性組合生成的向量空間 X 程度不超過並包括 d+1, 稱為“[或最高]度的多項式 d+1X, " 將多項式的向量空間擴展到度 dX.

多項式回歸的使用

通常,多項式回歸是探索性的,因為我們一開始並不知道要包括哪些單項式。從單項式創建新模型矩陣並重新擬合回歸的過程可能需要重複多次,在某些機器學習設置中可能需要重複很多次。

這種方法的主要問題是

  1. 單項式經常在新模型矩陣中引入有問題的“多重共線性”,主要是因為單個變量的冪往往是高度共線性的。(兩個不同變量的冪之間的共線性是不可預測的,因為它取決於這些變量如何相關,因此難以預測。)
  2. 僅更改模型矩陣的一列,或引入新的一列,或刪除一列,可能需要回歸過程的“冷啟動”,這可能需要很長時間的計算。

多項式代數的分級提供了克服這兩個問題的方法。

一個變量中的正交多項式

給定一個單列向量 X, 一組“正交多項式” X 是一系列列向量 p0(X),p1(X),p2(X), 形成為單項式的線性組合 X 單獨——也就是說,作為 X –具有以下屬性:

  1. 對於每個學位 d=0,1,2,, 向量 p0(X),p1(X),,pd(X) 生成相同的向量空間 X0,X1,,Xd. (請注意 X0 是個 n -向量和 X1 只是 X 本身。)
  2. pi(X) 是相互正交的,因為對於 ij, pi(X)pj(X)=0.

通常,替換模型矩陣P=(p0(X)p1(X)pd(X))

通過將其列標準化為單位長度,選擇由這些單項式形成的正交:$$ \mathbb{P}^\prime \mathbb{P} = \mathbb{I}{d+1}. $$ 因為倒數 PP 出現在大多數回歸方程和單位矩陣的逆 $ \mathbb{I}{d+1} $ 是本身,這代表了巨大的計算增益。

正交性幾乎決定了 pi(X). 您可以通過構造看到這一點:

  • 第一個多項式, p0(X), 必須是的倍數 n -向量 1=(1,1,,1) 單位長度。只有兩種選擇, ±1/n1. 通常選擇正平方根。
  • 第二個多項式, p1(X), 必須正交於 1. 可以通過回歸得到 X 反對 1, 其解是平均值的向量 ˆX=ˉX1. 如果殘差 ϵ=XˆX 不完全為零,它們只給出兩種可能的解決方案 p1(X)=±(1/||ϵ||),ϵ.

  • 一般來說, pd+1(X) 通過回歸得到 Xd+1 反對 p0(X),p1(X),,pd(X) 並將殘差重新縮放為單位長度的向量。當殘差不全為零時,有兩種符號選擇。否則,該過程結束:查看任何更高的權力都是徒勞的 X. (這是一個很好的定理,但它的證明不需要在這里分散我們的注意力。)

這是應用於向量內在序列的Gram-Schmidt 過程 X0,X1,,Xd,. 通常它是使用QR 分解來計算的,這幾乎是相同的,但以數值穩定的方式計算。

這種構造產生了一系列附加列,以考慮包含在模型矩陣中。因此,一個變量的多項式回歸通常通過依次添加該序列的元素來進行,直到不再獲得回歸的改進。因為每個新列都與之前的列正交,包括它不會改變任何之前的係數估計。這使得程序高效且易於解釋。

多變量多項式

探索性回歸(以及模型擬合)通常首先考慮模型中包含哪些(原始)變量;然後評估這些變量是否可以通過包括它們的各種變換來增加,例如單項式;然後引入由這些變量及其重新表達的產物形成的“相互作用”。

那麼,執行這樣的程序將首先在 的列中形成**單變量正交多項式 X 分別地。 在為每列選擇合適的度數後,您將引入交互作用。

在這一點上,部分單變量程序崩潰了。在確定合適的模型之前,您將應用什麼樣的交互順序?此外,既然我們已經真正進入了多變量分析的領域,可用選項的數量及其日益增加的複雜性表明,在構建一系列多變量正交多項式時收益可能會遞減。 但是,如果您有這樣的序列,您可以使用 QR 分解來計算它。


R做什麼

因此,多項式回歸軟件傾向於專注於計算單變量正交多項式序列。R將這種支持盡可能自動地擴展到單變量多項式組的特點是。這是做什麼poly的。(它的伴侶polym本質上是相同的代碼,花里胡哨的東西更少;這兩個函數做同樣的事情。)

具體來說,poly將在給定單個向量時計算一系列單變量正交多項式 X, 停在指定的程度 d. (如果 d 太大 - 並且很難預測太大 - 不幸的是會引發錯誤。)當給定一向量時 X1,,Xk 以矩陣的形式 X, 它會回來

  1. 正交多項式的序列 p1(Xj),p2(Xj),,pd(Xj) 對於每個 j 達到要求的最大程度 d. (由於常數向量 p0(Xi) 對所有變量都是通用的,而且非常簡單——它通常由回歸中的截距來適應——R不需要包含它。)
  2. 那些正交多項式之間的所有相互作用,直到並包括度數 d.

步驟(2)涉及幾個微妙之處。 通常,變量之間的“相互作用”是指“所有可能的產品”,但其中一些可能的產品的度數大於 d. 例如,與 2 變量和 d=2, R計算

p1(X1),p2(X1),p1(X2),p1(X1)p1(X2),p2(X2).

R不包括更高程度的交互 p2(X1)p1(X2), p1(X1)p2(X2) (3次多項式)或 p1(X2)p2(X2) (4 次多項式)。(這不是一個嚴重的限制,因為您可以輕鬆地自己計算這些產品或在回歸對formula像中指定它們。)

另一個微妙之處是,沒有對任何多變量產品應用任何標準化。 在示例中,唯一的此類產品是 p1(X1)p1(X2). 但是,即使它的平均值為零也不能保證,而且幾乎肯定不會有單位範數。從這個意義上說,這是一種真正的“互動” p1(X1)p1(X2) 因此可以解釋為交互通常在回歸模型中。

一個例子

讓我們看一個例子。 我隨機生成了一個矩陣X=(13 56 24).

為了使計算更容易理解,所有內容都四捨五入為兩位有效數字以進行顯示。

第一列的正交多項式序列 X1=(1,5,2) 從規範化開始 X01=(1,1,1) 到單位長度,給出 p0(X1)=(1,1,1)/3(0.58,0.58,0.58). 下一步包括 X11=X1 本身。使其正交於 p0(X1), 回歸 X1 反對 p0(X1) 並設置 p1(X1) 等於該回歸的殘差,重新調整為具有單位長度。結果是通常的標準化 X1 通過將其重新置中並除以其標準偏差獲得, p1(X1)=(0.57,0.79,0.23). 最後, X21=(1,25,4) 回歸到 p0(X1)p1(X1) 並且這些殘差被重新調整為單位長度。我們不能走得更遠,因為 X1 不能生成大於 n=3 方面。(我們之所以能走到這一步,是因為 X1,(t1)(t5)(t4), 有學位 3, 證明所有度的單項式 3 或更大是較低功率的線性組合,而這些較低功率是線性獨立的。)

結果矩陣表示一個正交多項式序列 X1

P1=(0.580.570.59 0.580.790.20 0.580.230.78)

(至兩位有效數字)。

以同樣的方式,一個正交多項式矩陣 X2

P2=(0.580.620.53 0.580.770.27 0.580.150.80).

The interaction term is the product of the middle columns of these matrices, equal to (0.35,0.61,0.035). The full matrix created by poly or polym, then, is

P=(0.570.590.620.350.53 0.790.200.770.610.27 0.230.780.150.0350.80).

Notice the sequence in which the columns are laid out: the non-constant orthonormal polynomials for X1 are in columns 1 and 2 while those for X2 are in columns 3 and 5. Thus, the only orthogonality that is guaranteed in this output is between these two pairs of columns. This is reflected in the calculation of PP, which will have zeros in positions (1,2),(2,1),(3,5), and (5,3) (shown in red below), *but may be nonzero anywhere else, and will have ones in positions (1,1),(2,2),(3,3), and (5,5) (shown in blue below), but is likely not to have a one in the other diagonal positions ( (4,4) in this example). Indeed,

P,P=(1010.280.091 010.0910.31 10.09110.250 0.280.30.250.50.32 0.091100.321).

When you inspect the P matrix shown in the question, and recognize that multiples of 1017 are really zeros, you will observe that this pattern of zeros in the red positions holds. This is the sense in which those bivariate polynomials are “orthogonal."

引用自:https://stats.stackexchange.com/questions/253123