机器算法验证 - 从 GAMLSS 中的系数和结创建样条 - 吾爱随笔录

从 GAMLSS 中的系数和结创建样条

机器算法验证样条无游戏

2022-03-29 21:42:55

在R包GAMLSS中，可以对随机变量进行建模 $Y$ 作为某个预测变量的平滑非参数函数 $x$ .

这种函数的一个选项是使用的惩罚样条y~pb(x)。这会输出一个系数和节点列表，结合一组基本样条曲线，得到一个平滑函数 $x$ .

给定系数和节点，如何重新创建样条函数？（最好不必编写我自己的 b 样条生成函数）。

例如：

library(gamlss)
set.seed(9876)
xstart <- 0
xend <- 100
datan <- 20000

seq <- seq(xstart, xend)
mean <- sapply(seq, function(x){0.5+0.2*sin(x/10)})
xs <- ceiling(runif(datan, xstart, xend))
ys <- sapply(xs, function(x){rnorm(1, mean = mean[x], sd = 0.1)})

m1 <- gamlss(ys~pb(xs))
plot(xs, ys)
lines(seq, mean, col="red")
lines(xs[order(xs)], fitted(m1)[order(xs)], col="green")

intercept <- m1$mu.coefficients[1] # 0.5495853
weight <- m1$mu.coefficients[2] # -0.0002851018
coefficients <- c(m1$mu.coefSmo[[1]]$coef) # c(-0.170704842, -0.066451626,  0.026591530,  0.119289203,  0.159657021,  0.149185418,  0.086505094,  0.003904402, -0.100156999, -0.188811997, -0.238717366, -0.237884900, -0.197802945, -0.090559794,  0.012576273,  0.101003289,  0.169210741,  0.181836117,  0.143883546,  0.061281663, -0.038608572, -0.136215586, -0.232871483)
knots <- m1$mu.coefSmo[[1]]$knots # c(-5.039,  0.010,  5.059, 10.108, 15.157, 20.206, 25.255, 30.304, 35.353, 40.402, 45.451, 50.500, 55.549, 60.598, 65.647, 70.696, 75.745, 80.794, 85.843, 90.892, 95.941 100.990 106.039)

我怎样才能获得绿色函数，只知道截距、权重、系数和结？我目前使用fitted(m1). 然而，这只是一个列表 $y$ 最初输入的列表的值 $x$ 值，它不是一个函数 $y$ 对于任何新的 $x$ .

2个回答

pb() 函数符合 Eilers 和 Marx (1996) 所述的 P 样条曲线：等间距节点上的 B 样条曲线和有限差分惩罚。在同一篇论文中，有一些代码块显示了如何拟合模型（如果我没记错的话，在附录中）。我不知道fitted.gamlss 方法的细节，但下面的代码应该会有所帮助（看图中的蓝线）。

为了计算 B 样条基，我使用 splines 包中的函数 splineDesign（如果我记得清楚的话，我在上面提到的参考中使用了相同的函数）。

要得到 $\hat{y}$ 对于一个新的价值 $x$ ，您可以只计算 splineDesign 函数的相应值并使用之前估计的系数（参见代码的最后一行和绿点）

# B-splines stuffs - observed xs
ndx   = 20
deg   = 3
xr    = max(xs)
xl    = min(xs)
xmax  = xr + 0.01 * (xr - xl)
xmin  = xl - 0.01 * (xr - xl)
dt    = (xmax - xmin) / ndx
knots = seq(xmin - deg * dt, xmax + deg * dt, by = dt)
B     = splineDesign(knots = knots, x = xs, ord = deg + 1, derivs = 0,outer.ok = TRUE)

# Compute smooth
ff    = intercept + weight * xs + B %*% coefficients 

# New x-value
xn    = 35
Bn    =  splineDesign(knots = knots, x = xn, ord = deg + 1, derivs = 0,outer.ok = TRUE)
fn    = intercept + weight * xn + Bn %*% coefficients 

# Plot Results
lines(xs[order(xs)], ff[order(xs)], col = 'blue', lty = 2, lwd = 2)
points(xn, fn, col = 'green', pch = 16, cex = 1.5)

这是一个生成新的函数调用 $y$ 任何值 $x$ .

predict(m1, newdata = data.frame(xs = new_x), data = data.frame(xs=xs))

请注意，我xs来自xs您的示例中的。通常，您需要修改它以predict使用原始名称向函数提供原始数据。

其它你可能感兴趣的问题

上一篇如何统计比较机器学习“回归”模型？下一篇L1 正则化强制凸函数的稀疏性