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