R smooth.construct.so.smooth.spec 皂膜平滑劑

說明

為皂膜平滑器設置基函數和擺動懲罰(Wood、Bravington 和 Hedley，2008)。皂膜平滑器基於將二維平滑構造為連接平滑變化的封閉邊界的肥皂膜的想法。除非進行非常嚴重的平滑處理，否則膠片會對數據產生扭曲。平滑設計為不平滑跨越邊界要素(例如半島)。

so版本設置完全流暢。 sf 版本僅設置邊界插值肥皂膜，而 sw 版本設置肥皂膜的擺動組件(邊界上為零)。後兩者對於用肥皂膜形成張量積非常有用，並且可以與 gamm 和 gamm4 一起使用。要使用它們來簡單地設置基礎，然後通過包裝器 smooth.construct2 或 smoothCon 進行調用。

用法

## S3 method for class 'so.smooth.spec'
smooth.construct(object,data,knots)
## S3 method for class 'sf.smooth.spec'
smooth.construct(object,data,knots)
## S3 method for class 'sw.smooth.spec'
smooth.construct(object,data,knots)

參數

細節

對於皂膜平滑，*必須*提供以下內容：

• k 每個邊界環平滑的基本尺寸。

• xt$bnd 平滑的邊界規範。

• 打結內部結的位置以使其光滑。

當在 GAM 中使用時，k 和 xt 通過 s 提供，而 knots 在 gam 的 knots 參數中提供。

xt 列表的 bnd 元素是列表(或數據幀)的列表，指定定義邊界的循環。每個邊界循環列表必須包含 2 列，給出定義邊界循環的點的坐標(當通過線段順序連接時)。循環不應相交(未檢查)。如果一個點位於奇數個邊界環的內部，則該點被視為位於感興趣區域中。每個邊界循環列表還可以包含列f，給出循環上的已知邊界條件。

xt 的 bndSpec 元素(如果非 NULL)應包含

• bs 要使用的循環平滑基礎的類型： "cc" 和 "cp" 之一。如果不是"cc"，則使用循環p-spline，並且必須提供參數m。

• knot.space 設置為 "even"，以在 "cc" 基礎上獲得均勻的結間距。

• m 1 或 2 元素數組，指定 "cp" 基礎和懲罰的順序。

目前，代碼不會處理超過一層的循環嵌套，也不會處理沒有外部封閉循環的單獨循環：如果存在已知的邊界條件(可識別性約束變得很尷尬)。

請注意，函數 locator 提供了一種以圖形方式定義邊界的簡單方法，在生成感興趣的域圖(右鍵單擊停止)後，使用類似 bnd <-as.data.frame(locator(type="l")) 的方法。如果真實邊界非常複雜，最好使用包圍真實邊界的更簡單的平滑邊界，它代表您不想平滑的主要邊界特征，但不遵循每個微小的細節。

模型建立和預測涉及評估被定義為偏微分方程解的基函數。使用稀疏矩陣方法在網格上對偏微分方程進行數值求解，並使用雙線性插值來獲取平滑域內任何位置的值。 PDE 解網格的維度可以通過作為 gam 公式中 s 的參數 xt 提供的列表的元素 nmax(默認 200)進行控製：它給出了最長使用的單元格數量網格側。

一點理論：肥皂膜光滑f(x,y)定義為

條件是 f=s 位於邊界曲線上，其中 s 是平滑函數(通常是循環懲罰回歸樣條)。函數 g 被定義為

其中g=0位於邊界曲線上，g(x_k,y_k)=c_k位於曲麵的‘knots’處； c_k 是模型係數。

在最簡單的情況下，f 的係數(邊界係數加上c_k)的估計是通過最小化

其中J_s通常是邊界平滑上的一些三次樣條類型擺動懲罰，J_f是(f_xx+f_yy)^2在邊界內部的積分。兩種懲罰都可以在模型係數中表示為二次形式。 \lambda 是平滑參數，可以通過 GCV、REML、AIC 等進行選擇。z 表示 f 的噪聲觀測值。

值

包含 object 的所有元素的列表加上

此外，還有通常由 smooth.construct 方法添加的所有元素。

警告

皂膜平滑器非常專業，並且比大多數平滑器需要更多設置(例如，您必須提供邊界和內部結，加上邊界平滑器基本尺寸)。值得一看參考。

例子

require(mgcv)

##########################
## simple test function...
##########################

fsb <- list(fs.boundary())
nmax <- 100
## create some internal knots...
knots <- data.frame(v=rep(seq(-.5,3,by=.5),4),
                    w=rep(c(-.6,-.3,.3,.6),rep(8,4)))
## Simulate some fitting data, inside boundary...
set.seed(0)
n<-600
v <- runif(n)*5-1;w<-runif(n)*2-1
y <- fs.test(v,w,b=1)
names(fsb[[1]]) <- c("v","w")
ind <- inSide(fsb,x=v,y=w) ## remove outsiders
y <- y + rnorm(n)*.3 ## add noise
y <- y[ind];v <- v[ind]; w <- w[ind] 
n <- length(y)

par(mfrow=c(3,2))
## plot boundary with knot and data locations
plot(fsb[[1]]$v,fsb[[1]]$w,type="l");points(knots,pch=20,col=2)
points(v,w,pch=".");

## Now fit the soap film smoother. 'k' is dimension of boundary smooth.
## boundary supplied in 'xt', and knots in 'knots'...
 
nmax <- 100 ## reduced from default for speed.
b <- gam(y~s(v,w,k=30,bs="so",xt=list(bnd=fsb,nmax=nmax)),knots=knots)

plot(b) ## default plot
plot(b,scheme=1)
plot(b,scheme=2)
plot(b,scheme=3)

vis.gam(b,plot.type="contour")

################################
# Fit same model in two parts...
################################

par(mfrow=c(2,2))
vis.gam(b,plot.type="contour")

b1 <- gam(y~s(v,w,k=30,bs="sf",xt=list(bnd=fsb,nmax=nmax))+
            s(v,w,k=30,bs="sw",xt=list(bnd=fsb,nmax=nmax)) ,knots=knots)
vis.gam(b,plot.type="contour")
plot(b1)
 
##################################################
## Now an example with known boundary condition...
##################################################

## Evaluate known boundary condition at boundary nodes...
fsb[[1]]$f <- fs.test(fsb[[1]]$v,fsb[[1]]$w,b=1,exclude=FALSE)

## Now fit the smooth...
bk <- gam(y~s(v,w,bs="so",xt=list(bnd=fsb,nmax=nmax)),knots=knots)
plot(bk) ## default plot

##########################################
## tensor product example...
##########################################

set.seed(9)
n <- 10000
v <- runif(n)*5-1;w<-runif(n)*2-1
t <- runif(n)
y <- fs.test(v,w,b=1)
y <- y + 4.2
y <- y^(.5+t)
fsb <- list(fs.boundary())
names(fsb[[1]]) <- c("v","w")
ind <- inSide(fsb,x=v,y=w) ## remove outsiders
y <- y[ind];v <- v[ind]; w <- w[ind]; t <- t[ind] 
n <- length(y)
y <- y + rnorm(n)*.05 ## add noise
knots <- data.frame(v=rep(seq(-.5,3,by=.5),4),
                    w=rep(c(-.6,-.3,.3,.6),rep(8,4)))

## notice NULL element in 'xt' list - to indicate no xt object for "cr" basis...
bk <- gam(y~ te(v,w,t,bs=c("sf","cr"),k=c(25,4),d=c(2,1),
          xt=list(list(bnd=fsb,nmax=nmax),NULL))+
          te(v,w,t,bs=c("sw","cr"),k=c(25,4),d=c(2,1),
	  xt=list(list(bnd=fsb,nmax=nmax),NULL)),knots=knots)

par(mfrow=c(3,2))
m<-100;n<-50 
xm <- seq(-1,3.5,length=m);yn<-seq(-1,1,length=n)
xx <- rep(xm,n);yy<-rep(yn,rep(m,n))
tru <- matrix(fs.test(xx,yy),m,n)+4.2 ## truth

image(xm,yn,tru^.5,col=heat.colors(100),xlab="v",ylab="w",
      main="truth")
lines(fsb[[1]]$v,fsb[[1]]$w,lwd=3)
contour(xm,yn,tru^.5,add=TRUE)

vis.gam(bk,view=c("v","w"),cond=list(t=0),plot.type="contour")

image(xm,yn,tru,col=heat.colors(100),xlab="v",ylab="w",
      main="truth")
lines(fsb[[1]]$v,fsb[[1]]$w,lwd=3)
contour(xm,yn,tru,add=TRUE)

vis.gam(bk,view=c("v","w"),cond=list(t=.5),plot.type="contour")

image(xm,yn,tru^1.5,col=heat.colors(100),xlab="v",ylab="w",
      main="truth")
lines(fsb[[1]]$v,fsb[[1]]$w,lwd=3)
contour(xm,yn,tru^1.5,add=TRUE)

vis.gam(bk,view=c("v","w"),cond=list(t=1),plot.type="contour")


#############################
# nested boundary example...
#############################

bnd <- list(list(x=0,y=0),list(x=0,y=0))
seq(0,2*pi,length=100) -> theta
bnd[[1]]$x <- sin(theta);bnd[[1]]$y <- cos(theta)
bnd[[2]]$x <- .3 + .3*sin(theta);
bnd[[2]]$y <- .3 + .3*cos(theta)
plot(bnd[[1]]$x,bnd[[1]]$y,type="l")
lines(bnd[[2]]$x,bnd[[2]]$y)

## setup knots
k <- 8
xm <- seq(-1,1,length=k);ym <- seq(-1,1,length=k)
x=rep(xm,k);y=rep(ym,rep(k,k))
ind <- inSide(bnd,x,y)
knots <- data.frame(x=x[ind],y=y[ind])
points(knots$x,knots$y)

## a test function

f1 <- function(x,y) {
  exp(-(x-.3)^2-(y-.3)^2)
}

## plot the test function within the domain 
par(mfrow=c(2,3))
m<-100;n<-100 
xm <- seq(-1,1,length=m);yn<-seq(-1,1,length=n)
x <- rep(xm,n);y<-rep(yn,rep(m,n))
ff <- f1(x,y)
ind <- inSide(bnd,x,y)
ff[!ind] <- NA
image(xm,yn,matrix(ff,m,n),xlab="x",ylab="y")
contour(xm,yn,matrix(ff,m,n),add=TRUE)
lines(bnd[[1]]$x,bnd[[1]]$y,lwd=2);lines(bnd[[2]]$x,bnd[[2]]$y,lwd=2)

## Simulate data by noisy sampling from test function...

set.seed(1)
x <- runif(300)*2-1;y <- runif(300)*2-1
ind <- inSide(bnd,x,y)
x <- x[ind];y <- y[ind]
n <- length(x)
z <- f1(x,y) + rnorm(n)*.1

## Fit a soap film smooth to the noisy data
nmax <- 60
b <- gam(z~s(x,y,k=c(30,15),bs="so",xt=list(bnd=bnd,nmax=nmax)),
         knots=knots,method="REML")
plot(b) ## default plot
vis.gam(b,plot.type="contour") ## prettier version

## trying out separated fits....
ba <- gam(z~s(x,y,k=c(30,15),bs="sf",xt=list(bnd=bnd,nmax=nmax))+
          s(x,y,k=c(30,15),bs="sw",xt=list(bnd=bnd,nmax=nmax)),
	  knots=knots,method="REML")
plot(ba)
vis.gam(ba,plot.type="contour")

作者

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

參考

Wood, S.N., M.V. Bravington and S.L. Hedley (2008) "Soap film smoothing", J.R.Statist.Soc.B 70(5), 931-955. doi:10.1111/j.1467-9868.2008.00665.x

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

也可以看看

Predict.matrix.soap.film