R magic 通过 GCV 或 UBRE 进行稳定的多重平滑参数估计

说明

通过 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)

参数

细节

该方法是一种计算高效的方法，将 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 例程执行矩阵分解。

值

该函数返回一个包含以下项目的列表：

请注意，可以使用 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