Contrasts

多項式對比變量的計算

  • July 8, 2013

請告訴我如何有效地將分類變量(因子)重新編碼為一組正交多項式對比變量。

對於許多類型的對比變量(例如偏差、簡單、Helmert 等),傳遞是:

  1. 組成與類型相對應的對比度係數矩陣。
  2. 逆或泛化逆它以獲得代碼矩陣。

例如:

Suppose there is 3-group factor and we want to recode it into a set of deviation  contrast variables.
The last group is treated as reference. Then the contrast coefficients matrix L is

        Group1 Group2 Group3
  var1   2/3   -1/3   -1/3
  var2  -1/3    2/3   -1/3

and ginv(L) is then the sought-for coding matrix

        var1 var2
 Group1   1    0
 Group2   0    1
 Group3  -1   -1

(We might also use inv(L) instead if we add a row for constant, equal to 1/3, at the head of L.)

是否有相同或相似的方法來獲得多項式對比變量?如果是的話,矩陣 C 會是什麼樣子以及如何組合它?如果沒有,那麼還有什麼方法可以有效地計算多項式對比變量(例如通過矩陣代數)。

作為我之前關於這個主題的帖子的延續,我想分享一些關於線性代數和相關 R 函數背後的函數的初步(儘管不完整)探索。這應該是一項正在進行的工作。

部分功能的不透明性與 Householder 的“緊湊”形式有關分解。Householder 分解背後的思想是在由單位向量確定的超平面上反映向量如下圖所示,但有目的地選擇這個平面,以便投影原始矩陣的每個列向量到標准單位向量。標準化 norm-2向量可用於計算不同的 Householder 變換.

在此處輸入圖像描述

得到的投影可以表示為

向量表示列向量之間的差異在矩陣中我們要分解的向量對應於通過子空間或“鏡像”確定的反射.

LAPACK 使用的方法通過將 Householder 反射器中的第一個條目轉換為的。而不是規範化向量到和, 只是第一個條目被轉換為; 然而,這些新的載體——稱它們為仍然可以用作方向向量。

該方法的美妙之處在於 在一個分解是上三角的,我們實際上可以利用中的元素在對角線下方用這些填充它們反射器。幸運的是,這些向量中的前導條目都相等,防止矩陣的“有爭議”對角線出現問題:知道它們都是它們不需要被包括在內,並且可以產生對角線的條目.

函數中的“compact QR”矩陣qr()$qr可以大致理解為矩陣和“修改”反射器的下三角“存儲”矩陣。

Householder 投影仍然具有形式,但我們不會合作(),而是使用向量,其中只有第一個條目被保證為, 和

.

有人會認為存儲這些就可以了對角線下方的反射器或不包括第一個條目,然後收工。然而,事情從未如此簡單。相反,存儲在對角線下方的qr()$qr是Householder 變換中的係數表示為 (1),因此,定義 作為:

, 反射器可以表示為. 這些“反射器”向量是存儲在下面的向量在所謂的“緊湊”。

現在我們距離向量,第一個條目不再是, 因此 will 的輸出qr()需要包含恢復它們的密鑰,因為我們堅持排除“反射器”向量的第一個條目以適合所有內容qr()$qr。所以我們看到輸出中的值?嗯,不,那是可以預測的。相反,在qr()$qraux(存儲此密鑰的位置)的輸出中,我們發現.

所以下面用紅色框起來,我們看到“反射器”(),不包括他們的第一個條目。

在此處輸入圖像描述

所有的代碼都在這裡,但由於這個答案是關於編碼和線性代數的交集,我將粘貼輸出以方便起見:


options(scipen=999)
set.seed(13)
(X = matrix(c(rnorm(16)), nrow=4, byrow=F))
          [,1]      [,2]       [,3]       [,4]
[1,]  0.5543269 1.1425261 -0.3653828 -1.3609845
[2,] -0.2802719 0.4155261  1.1051443 -1.8560272
[3,]  1.7751634 1.2295066 -1.0935940 -0.4398554
[4,]  0.1873201 0.2366797  0.4618709 -0.1939469

現在我寫的函數House()如下:

  House = function(A){
   Q = diag(nrow(A))
   reflectors = matrix(0,nrow=nrow(A),ncol=ncol(A))
   for(r in 1:(nrow(A) - 1)){ 
       # We will apply Householder to progressively the columns in A, decreasing 1 element at a time.
       x = A[r:nrow(A), r] 
       # We now get the vector v, starting with first entry = norm-2 of x[i] times 1
       # The sign is to avoid computational issues
       first = (sign(x[1]) * sqrt(sum(x^2))) +  x[1]
       # We get the rest of v, which is x unchanged, since e1 = [1, 0, 0, ..., 0]
       # We go the the last column / row, hence the if statement:
       v = if(length(x) > 1){c(first, x[2:length(x)])}else{v = c(first)}
       # Now we make the first entry unitary:
       w = v/first
       # Tau will be used in the Householder transform, so here it goes:
       t = as.numeric(t(w)%*%w) / 2
       # And the "reflectors" are stored as in the R qr()$qr function:
       reflectors[r: nrow(A), r] = w/t
       # The Householder tranformation is:
       I = diag(length(r:nrow(A)))
       H.transf = I - 1/t * (w %*% t(w))
       H_i  = diag(nrow(A))
       H_i[r:nrow(A),r:ncol(A)] = H.transf
       # And we apply the Householder reflection - we left multiply the entire A or Q
       A = H_i %*% A
       Q = H_i %*% Q
   }
   DECOMPOSITION = list("Q"= t(Q), "R"= round(A,7), 
           "compact Q as in qr()$qr"=  
           ((A*upper.tri(A,diag=T))+(reflectors*lower.tri(reflectors,diag=F))), 
           "reflectors" = reflectors,
           "rho"=c(apply(reflectors[,1:(ncol(reflectors)- 1)], 2, 
               function(x) sum(x^2) / 2), A[nrow(A),ncol(A)]))
   return(DECOMPOSITION)
}

讓我們將輸出與 R 內置函數進行比較。首先是自製功能:

(H = House(X))
$Q
           [,1]        [,2]       [,3]       [,4]
[1,] -0.29329367 -0.73996967  0.5382474  0.2769719
[2,]  0.14829152 -0.65124800 -0.5656093 -0.4837063
[3,] -0.93923665  0.13835611 -0.1947321 -0.2465187
[4,] -0.09911084 -0.09580458 -0.5936794  0.7928072

$R
         [,1]       [,2]       [,3]      [,4]
[1,] -1.890006 -1.4517318  1.2524151 0.5562856
[2,]  0.000000 -0.9686105 -0.6449056 2.1735456
[3,]  0.000000  0.0000000 -0.8829916 0.5180361
[4,]  0.000000  0.0000000  0.0000000 0.4754876

$`compact Q as in qr()$qr`
           [,1]        [,2]       [,3]      [,4]
[1,] -1.89000649 -1.45173183  1.2524151 0.5562856
[2,] -0.14829152 -0.96861050 -0.6449056 2.1735456
[3,]  0.93923665 -0.67574886 -0.8829916 0.5180361
[4,]  0.09911084  0.03909742  0.6235799 0.4754876

$reflectors
           [,1]        [,2]      [,3] [,4]
[1,]  1.29329367  0.00000000 0.0000000    0
[2,] -0.14829152  1.73609434 0.0000000    0
[3,]  0.93923665 -0.67574886 1.7817597    0
[4,]  0.09911084  0.03909742 0.6235799    0

$rho
[1] 1.2932937 1.7360943 1.7817597 0.4754876

到 R 函數:

qr.Q(qr(X))
           [,1]        [,2]       [,3]       [,4]
[1,] -0.29329367 -0.73996967  0.5382474  0.2769719
[2,]  0.14829152 -0.65124800 -0.5656093 -0.4837063
[3,] -0.93923665  0.13835611 -0.1947321 -0.2465187
[4,] -0.09911084 -0.09580458 -0.5936794  0.7928072

qr.R(qr(X))
         [,1]       [,2]       [,3]      [,4]
[1,] -1.890006 -1.4517318  1.2524151 0.5562856
[2,]  0.000000 -0.9686105 -0.6449056 2.1735456
[3,]  0.000000  0.0000000 -0.8829916 0.5180361
[4,]  0.000000  0.0000000  0.0000000 0.4754876

$qr
           [,1]        [,2]       [,3]      [,4]
[1,] -1.89000649 -1.45173183  1.2524151 0.5562856
[2,] -0.14829152 -0.96861050 -0.6449056 2.1735456
[3,]  0.93923665 -0.67574886 -0.8829916 0.5180361
[4,]  0.09911084  0.03909742  0.6235799 0.4754876

$qraux
[1] 1.2932937 1.7360943 1.7817597 0.4754876

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

comments powered by Disqus