R XWXd 用于离散模型矩阵处理的内部函数

说明

使用离散模型矩阵进行计算的例程，如 Wood 等人所述。 (2017)以及李和伍德(2019)。

用法

XWXd(X,w,k,ks,ts,dt,v,qc,nthreads=1,drop=NULL,ar.stop=-1,ar.row=-1,ar.w=-1,
     lt=NULL,rt=NULL)
XWyd(X,w,y,k,ks,ts,dt,v,qc,drop=NULL,ar.stop=-1,ar.row=-1,ar.w=-1,lt=NULL)
Xbd(X,beta,k,ks,ts,dt,v,qc,drop=NULL,lt=NULL)
diagXVXd(X,V,k,ks,ts,dt,v,qc,drop=NULL,nthreads=1,lt=NULL,rt=NULL)

参数

X	包含完整模型矩阵项的模型矩阵的唯一行或项边距的模型矩阵的矩阵列表。如果术语子集参数 lt 和 rt 为非 NULL，则需要 "lpip" 属性：请参阅详细信息。 X 的元素可能是类 "dgCMatrix" 的稀疏矩阵，在这种情况下，列表需要定义反向索引的属性 "r" 和 "off"(参见详细信息)。
w	n-vector 权重
y	n-vector 数据。
beta	系数向量。
k	一个矩阵，其列的索引为n-vectors，每个列选择创建完整矩阵所需的 X[[i]] 的行。
ks	第 i 项具有索引向量 ks[i,1]:(ks[i,2]-1) 。相应的完整模型矩阵被求和。
ts	X 中每个模型项开始的元素。
dt	X 的每个项有多少个元素。
v	v[[i]] 是第 i 个学期的 Householder 向量，如果 qc[i]>0 。
qc	如果 qc[i]>0 则术语有约束。
nthreads	使用的线程数
drop	要删除的模型矩阵/参数的列列表
ar.stop	消极的忽略。否则，对 ar.row 选择的行中的 (ar.stop[i-1]+1):ar.stop[i] 行进行求和，并按 ar.w 进行加权，以获得要使用的模型矩阵的第 i 行。
ar.row	提取这些行...
ar.w	对这些权重进行加权，并根据 ar.stop 求和。用于实现AR模型。
lt	仅使用与这些模型矩阵项相对应的 X 列(对于 XWXd 中的左手 X )。如果 NULL 设置为 rt 。
rt	作为右手 lt X 。如果 NULL 设置为 lt 。如果lt 和rt 是NULL，则使用所有列。
V	系数协方差矩阵。

细节

这些函数实际上是内部函数，但被导出，以便它们可以毫无问题地在系列的初始化代码中使用。它们主要由 bam 用于实现参考文献中给出的方法。 XWXd 生成 X^TWX ， XWy 生成 X^TWy ， Xbd 生成 X\beta ，diagXVXd 生成 XVX^T 的对角线。

X 的"lpip" 属性是每一项的系数索引的列表。如果通过 lt 和 rt 进行子集化，则为必需。

X 可以是 "dgCMatrix" 类的稀疏矩阵列表，在这种情况下需要反向索引，将存储的矩阵行映射到完整矩阵中的行(这与 k 相反，它将完整矩阵行映射到存储唯一的矩阵行)。 r 与 k 具有相同的维度，而 off 是一个列表，其元素数量与 k 的列数量相同。 r 和 off 作为属性提供给 X 。为了简单起见，让r和off表示彼此对应的单个列和元素：然后r[off[j]:(off[j+1]-1)]包含与存储矩阵的行j相对应的完整矩阵的行。反向索引对于稀疏矩阵的高效计算至关重要。请参阅示例代码，了解如何从前向索引矩阵 k 高效创建它们。

例子

  library(mgcv);library(Matrix)
  ## simulate some data creating a marginal matrix sequence...
  set.seed(0);n <- 4000
  dat <- gamSim(1,n=n,dist="normal",scale=2)
  dat$x4 <- runif(n)
  dat$y <- dat$y + 3*exp(dat$x4*15-5)/(1+exp(dat$x4*15-5))
  dat$fac <- factor(sample(1:20,n,replace=TRUE))
  G <- gam(y ~ te(x0,x2,k=5,bs="bs",m=1)+s(x1)+s(x4)+s(x3,fac,bs="fs"),
           fit=FALSE,data=dat,discrete=TRUE)
  p <- ncol(G$X)
  ## create a sparse version...
  Xs <- list(); r <- G$kd*0; off <- list()
  for (i in 1:length(G$Xd)) Xs[[i]] <- as(G$Xd[[i]],"dgCMatrix")
  for (j in 1:nrow(G$ks)) { ## create the reverse indices...
    nr <- nrow(Xs[[j]]) ## make sure we always tab to final stored row 
    for (i in G$ks[j,1]:(G$ks[j,2]-1)) {
      r[,i] <- (1:length(G$kd[,i]))[order(G$kd[,i])]
      off[[i]] <- cumsum(c(1,tabulate(G$kd[,i],nbins=nr)))-1
    }
  }
  attr(Xs,"off") <- off;attr(Xs,"r") <- r 

  par(mfrow=c(2,3))

  beta <- runif(p)
  Xb0 <- Xbd(G$Xd,beta,G$kd,G$ks,G$ts,G$dt,G$v,G$qc)
  Xb1 <- Xbd(Xs,beta,G$kd,G$ks,G$ts,G$dt,G$v,G$qc)
  range(Xb0-Xb1);plot(Xb0,Xb1,pch=".")

  bb <- cbind(beta,beta+runif(p)*.3)
  Xb0 <- Xbd(G$Xd,bb,G$kd,G$ks,G$ts,G$dt,G$v,G$qc)
  Xb1 <- Xbd(Xs,bb,G$kd,G$ks,G$ts,G$dt,G$v,G$qc)
  range(Xb0-Xb1);plot(Xb0,Xb1,pch=".")
  
  w <- runif(n)
  XWy0 <- XWyd(G$Xd,w,y=dat$y,G$kd,G$ks,G$ts,G$dt,G$v,G$qc) 
  XWy1 <- XWyd(Xs,w,y=dat$y,G$kd,G$ks,G$ts,G$dt,G$v,G$qc)
  range(XWy1-XWy0);plot(XWy1,XWy0,pch=".")

  yy <- cbind(dat$y,dat$y+runif(n)-.5)
  XWy0 <- XWyd(G$Xd,w,y=yy,G$kd,G$ks,G$ts,G$dt,G$v,G$qc) 
  XWy1 <- XWyd(Xs,w,y=yy,G$kd,G$ks,G$ts,G$dt,G$v,G$qc)
  range(XWy1-XWy0);plot(XWy1,XWy0,pch=".")

  A <- XWXd(G$Xd,w,G$kd,G$ks,G$ts,G$dt,G$v,G$qc)
  B <- XWXd(Xs,w,G$kd,G$ks,G$ts,G$dt,G$v,G$qc)
  range(A-B);plot(A,B,pch=".")

  V <- crossprod(matrix(runif(p*p),p,p))
  ii <- c(20:30,100:200)
  jj <- c(50:90,150:160)
  V[ii,jj] <- 0;V[jj,ii] <- 0
  d1 <- diagXVXd(G$Xd,V,G$kd,G$ks,G$ts,G$dt,G$v,G$qc)
  Vs <- as(V,"dgCMatrix")
  d2 <- diagXVXd(Xs,Vs,G$kd,G$ks,G$ts,G$dt,G$v,G$qc)
  range(d1-d2);plot(d1,d2,pch=".")

作者

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

参考

Wood, S.N., Li, Z., Shaddick, G. & Augustin N.H. (2017) Generalized additive models for gigadata: modelling the UK black smoke network daily data. Journal of the American Statistical Association. 112(519):1199-1210 doi:10.1080/01621459.2016.1195744

Li, Z & S.N. Wood (2019) Faster model matrix crossproducts for large generalized linear models with discretized covariates. Statistics and Computing. doi:10.1007/s11222-019-09864-2