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