R magic 通過 GCV 或 UBRE 進行穩定的多重平滑參數估計

R語言 magic 位於 mgcv 包(package)。

說明

通過 GCV 或 UBRE 有效估計具有多重(二次)懲罰的廣義嶺回歸問題中的平滑參數的函數。該函數使用 multi-dimensions 中的牛頓法，並通過最速下降法來迭代調整每個懲罰的平滑參數(一個懲罰的平滑參數可能固定為 1)。

為了獲得最大數值穩定性，該方法基於正交分解方法，並嘗試使用截斷奇異值分解方法來優雅地處理數值秩不足。

用法

magic(y,X,sp,S,off,L=NULL,lsp0=NULL,rank=NULL,H=NULL,C=NULL,
      w=NULL,gamma=1,scale=1,gcv=TRUE,ridge.parameter=NULL,
      control=list(tol=1e-6,step.half=25,rank.tol=
      .Machine$double.eps^0.5),extra.rss=0,n.score=length(y),nthreads=1)

參數

`y`	是響應數據向量。
`X`	是模型矩陣(允許的列數多於行數)。
`sp`	是平滑參數的數組。向量 `L%*%log(sp) + lsp0` 包含實際乘以 `S` 中存儲的懲罰矩陣的平滑參數的對數(如果 `NULL` 則將 `L` 視為恒等式)。任何負的 `sp` 值都會自動初始化，否則它們將被視為提供起始值。如果 GCV/UBRE 分數的梯度在提供的值處太小，則提供的起始值將重置為默認起始值。
`S`	是懲罰矩陣的列表。 `S[[i]]` 是第 i 個懲罰矩陣，但請注意，它不是存儲為完整矩陣，而是存儲為包括懲罰矩陣的所有非零元素的最小方陣。 `S[[i]]` 的元素 1,1 占據第 i 個懲罰矩陣的元素 `off[i]` 、 `off[i]` 。每個`S[[i]]` 必須是半正定的。如果沒有要估計的平滑參數，則設置為`list()`。
`off`	是一個數組，指示參數向量中的第一個參數，該參數受到涉及 `S[[i]]` 的懲罰。
`L`	是一個矩陣，將 `log(sp)` 映射到對數平滑參數，該參數實際上乘以 `S` 的元素定義的懲罰。視為身份，如果 `NULL` 。請參閱上麵`sp` 下的內容。
`lsp0`	如果`L` 不是`NULL`，則這是從`log(sp)` 到實際對數平滑參數的線性變換中的常量向量。因此，平滑參數乘以 `S[[i]]` 的對數由 `L%*%log(sp) + lsp0` 給出。如果 `NULL` 則視為 0。
`rank`	是一個指定懲罰等級的數組。這對於形成罰矩陣的平方根很有用，但不是必需的。
`H`	是可選的偏移懲罰 - 即平滑參數固定為 1 的懲罰。這對於允許估計過程的正則化、固定平滑懲罰等很有用。
`C`	是指定擬合問題的任何線性等式約束的可選矩陣。如果 `\bf b` 是參數向量，則參數強製滿足 `{\bf Cb} = {\bf 0}` 。
`w`	回歸權重。如果這是一個矩陣，則將其視為 `y` (特別是 `{\bf V}_y^{-1} = {\bf w}^\prime{\bf w}` )協方差矩陣的逆矩陣的平方根。如果 `w` 是一個數組，則將其視為該矩陣的對角線，或者隻是 `y` 每個元素的權重。請參閱下麵的使用此示例的示例。
`gamma`	是 GCV 或 UBRE 分數中模型自由度的膨脹因子。
`scale`	是與 UBRE 一起使用的比例參數。
`gcv`	如果要使用GCV，則應設置為`TRUE`，對於UBRE，應設置為`FALSE`。
`ridge.parameter`	有時對擬合問題應用嶺罰分很有用，直接懲罰約束空間中的參數。將此參數設置為大於零的值將導致使用此類懲罰，其大小由參數值給出。
`control`	是具有以下元素的迭代控製常量列表： tol 用於判斷收斂性的容差。 step.half 如果嘗試步驟失敗，則該方法會嘗試將其減半，最多 `step.half` 次。 rank.tol 是一個常量，用於測試問題的數值排名缺陷。本質上任何小於 `rank_tol` 乘以問題的最大奇異值的奇異值都設置為零。
`extra.rss`	是在計算 GCV、UBRE 和尺度參數估計時添加到殘差平方和(平方範數)項的常數。與 `n.score` 結合使用，這對於處理非常大的數據集的某些方法很有用。
`n.score`	在 GCV/UBRE 分數計算中用作數據數量的數字：通常是實際數據數量，但有一些處理非常大的數據集的方法可以改變這一點。
`nthreads`	如果設置為 >1，`magic` 可以使用多個線程。

細節

該方法是一種計算高效的方法，將 GCV 或 UBRE(通常近似 AIC)應用於以下形式的廣義嶺回歸問題中的平滑參數選擇問題：

可能受到約束 {\bf Cb}={\bf 0} 。 {\bf X} 是設計矩陣，\bf b 參數向量，\bf y 數據向量，\bf W 權重矩陣， {\bf S}_i 正半正定係數矩陣，定義帶有相關平滑參數的第 i 個罰分 \theta_i 、\bf H 是正半定偏移懲罰矩陣，\bf C 是定義問題的任何線性等式約束的係數矩陣。 {\bf X} 不必具有完整的列等級。

選擇 \theta_i 是為了最小化 GCV 分數：

或 UBRE 分數：

其中 \gamma 是 gamma 自由度的膨脹因子(通常設置為 1)，\phi 是 scale 尺度參數。 \bf A 是擬合問題的帽子矩陣(影響矩陣)(即將數據映射到擬合值的矩陣)。分數對平滑參數的依賴性是通過 \bf A 實現的。

該方法通過牛頓或 \theta_i 日誌的最速下降更新進行操作。該方法的一個關鍵方麵是對分數的一階和二階導數進行穩定且經濟的計算。對數平滑參數。因為 GCV/UBRE 分數與其他分數持平。非常大或非常小 \theta_i ，獲得良好的起始參數非常重要，並且要小心不要進入平滑參數空間的平坦區域。因此，該算法會重新調整任何可能導致 log(\theta_i) 變化超過 5 的牛頓步。僅當 GCV/UBRE 的 Hessian 矩陣為正定時才使用牛頓步，否則使用最速下降。如果牛頓步必須收縮得太遠(表明牛頓基礎的二次模型很差)，則使用類似的最速下降。所有初始最速下降步驟都經過縮放，以便其最大分量為 1。然而，計算步驟時，如果成功，則永遠不會擴展(以避免目標的平坦部分)，但如果步驟不減少GCV/UBRE 評分，直到他們這樣做為止，否則該方向被視為失敗。 (給定平滑參數，很容易找到最佳的 \bf b 參數。)

該方法使用 C 進行編碼，並使用 LINPACK 和 LAPACK 例程執行矩陣分解。

值

該函數返回一個包含以下項目的列表：

`b`	給定估計的平滑參數的最佳擬合參數。
`scale`	估計的 (GCV) 或提供的 (UBRE) 尺度參數。
`score`	最小化的 GCV 或 UBRE 分數。
`sp`	估計的平滑參數的數組。
`sp.full`	實際上與 `S` 的元素相乘的平滑參數數組(如果 `L` 是 `NULL` 則與 `sp` 相同)。這是`exp(L%*%log(sp))`。
`rV`	參數協方差矩陣的因式分解形式。參數 `b` 的(貝葉斯)協方差矩陣由 `rV%%t(rV)scale` 給出。
`gcv.info`	是有關方法性能的信息列表，包含以下元素： full.rank 問題的表觀等級：參數數量減去等式約束的數量。秩問題的估計實際等級(在方法的最終迭代時)。 fully.converged 如果該方法通過滿足收斂標準而收斂，則為`TRUE`；如果該方法因未能沿搜索方向降低分數而被覆蓋，則為`FALSE`。 hess.pos.def 如果 UBRE 或 GCV 分數的 hessian 矩陣在收斂時是正定的，則為 `TRUE`。 iter 是牛頓/最速下降迭代次數。 score.calls 是必須評估 GCV/UBRE 分數的次數。 rms.grad 是 UBRE/GCV 分數相對於梯度的均方根。平滑參數。 R 來自加權模型矩陣的 QR 分解的因子 R。這是 un-pivoted ，因此列順序對應於 `X` 。所以它可能不是上三角形。

請注意，可以使用 magic.post.proc 獲得一些進一步有用的量。

例子

## Use `magic' for a standard additive model fit ... 
   library(mgcv)
   set.seed(1);n <- 200;sig <- 1
   dat <- gamSim(1,n=n,scale=sig)
   k <- 30
## set up additive model
   G <- gam(y~s(x0,k=k)+s(x1,k=k)+s(x2,k=k)+s(x3,k=k),fit=FALSE,data=dat)
## fit using magic (and gam default tolerance)
   mgfit <- magic(G$y,G$X,G$sp,G$S,G$off,rank=G$rank,
                  control=list(tol=1e-7,step.half=15))
## and fit using gam as consistency check
   b <- gam(G=G)
   mgfit$sp;b$sp  # compare smoothing parameter estimates
   edf <- magic.post.proc(G$X,mgfit,G$w)$edf # get e.d.f. per param
   range(edf-b$edf)  # compare

## p>n example... fit model to first 100 data only, so more
## params than data...

   mgfit <- magic(G$y[1:100],G$X[1:100,],G$sp,G$S,G$off,rank=G$rank)
   edf <- magic.post.proc(G$X[1:100,],mgfit,G$w[1:100])$edf

## constrain first two smooths to have identical smoothing parameters
   L <- diag(3);L <- rbind(L[1,],L)
   mgfit <- magic(G$y,G$X,rep(-1,3),G$S,G$off,L=L,rank=G$rank,C=G$C)

## Now a correlated data example ... 
    library(nlme)
## simulate truth
    set.seed(1);n<-400;sig<-2
    x <- 0:(n-1)/(n-1)
    f <- 0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10
## produce scaled covariance matrix for AR1 errors...
    V <- corMatrix(Initialize(corAR1(.6),data.frame(x=x)))
    Cv <- chol(V)  # t(Cv)%*%Cv=V
## Simulate AR1 errors ...
    e <- t(Cv)%*%rnorm(n,0,sig) # so cov(e) = V * sig^2
## Observe truth + AR1 errors
    y <- f + e 
## GAM ignoring correlation
    par(mfrow=c(1,2))
    b <- gam(y~s(x,k=20))
    plot(b);lines(x,f-mean(f),col=2);title("Ignoring correlation")
## Fit smooth, taking account of *known* correlation...
    w <- solve(t(Cv)) # V^{-1} = w'w
    ## Use `gam' to set up model for fitting...
    G <- gam(y~s(x,k=20),fit=FALSE)
    ## fit using magic, with weight *matrix*
    mgfit <- magic(G$y,G$X,G$sp,G$S,G$off,rank=G$rank,C=G$C,w=w)
## Modify previous gam object using new fit, for plotting...    
    mg.stuff <- magic.post.proc(G$X,mgfit,w)
    b$edf <- mg.stuff$edf;b$Vp <- mg.stuff$Vb
    b$coefficients <- mgfit$b 
    plot(b);lines(x,f-mean(f),col=2);title("Known correlation")

作者

Simon N. Wood simon.wood@r-project.org

參考

Wood, S.N. (2004) Stable and efficient multiple smoothing parameter estimation for generalized additive models. J. Amer. Statist. Ass. 99:673-686

https://www.maths.ed.ac.uk/~swood34/

也可以看看

magic.post.proc、gam

相關用法

注：本文由純淨天空篩選整理自R-devel大神的英文原創作品 Stable Multiple Smoothing Parameter Estimation by GCV or UBRE。非經特殊聲明，原始代碼版權歸原作者所有，本譯文未經允許或授權，請勿轉載或複製。