设置

让描述一个连续的预测变量（例如年龄）。让是一个随机变量（例如高度），它是的某个函数。$x$$Y$$x$

数据由个点组成，每个点都是和的组合（例如）。$n$$x$$y$$data_i = (age_i, height_i)$

的分布可能是非正态的，具有一定程度的偏斜。不仅的平均值，而且它的方差和偏度也可能是的函数（例如，成年人平均比儿童高，但身高也更加多样化，并且异常高的异常值也更多）。$Y$$Y$$x$

使用 GAMLSS 的相关方法

我知道，如果假设一个分布，则 GAMLSS可用于描述分布的参数，其中每个参数都建模为的函数。这些函数可以作为多项式或样条给出，可能使用一些惩罚来平滑。$Y$$x$

$$Y \sim \mathcal{D}(\mu, \sigma)\\ g_1(\mu) = s_1(x)\\ g_2(\sigma) = s_2(x)$$

我的问题

然而，这些代表所选分布的参数的函数。是否也可以在不指定特定分布的情况下获得数据的矩函数（例如均值、方差、偏度）？

$$mean(Y) = s_3(x)\\ var(Y) = s_4(x)\\ skew(Y) = s_5(x)$$

我猜一个移动平均线，同样，一个移动方差和移动偏度，将给出一个类似的函数。但我想利用GAMLSS的优化（例如 ML、GAIC 或 GCV），包括对过度拟合的一些惩罚。如果这存在。如果在不首先指定分布的情况下甚至有意义。$x$

例子

作为一个最小的工作示例，我们将生成我们知道时刻的数据。

demo.BSplines首先，我们从库中改编gamlss.demos以创建一个生成随机样条的函数。

library(gamlss)
library(gamlss.demo)
print(demo.BSplines)

tpower <- function(x, t, p) (x - t)^p * (x > t)

bbase <- function(x, xl=min(x), xr=max(x), nseg=10, deg=3) {
  dx <- (xr - xl)/nseg
  knots <- seq(xl - deg * dx, xr + deg * dx, by = dx)
  P <- outer(x, knots, tpower, deg)
  n <- dim(P)[2]
  D <- diff(diag(n), diff = deg + 1)/(gamma(deg + 1) * dx^deg)
  B <- (-1)^(deg + 1) * P %*% t(D)
  return(B)
}

bs.random <- function(nseg=5, bdeg=3, xlim=100) {
  x <- seq(0, xlim)
  B <- bbase(x, nseg = nseg, deg = bdeg)
  a <- runif(ncol(B))
  z <- B %*% a
  return(z)
}

让我们生成两个 B 样条。

set.seed(9876)
nseg <- 5
bdeg <- 3
xlim <- 100
datan <- 20000

mu <- bs.random(nseg=nseg, bdeg=bdeg, xlim=xlim)
sigma <- bs.random(nseg=nseg, bdeg=bdeg, xlim=xlim)
plot(NULL, xlim=c(0,100), ylim=c(0,1), xlab="x", ylab="y", main="Random B-splines")
lines(mu, col="blue")
lines(sigma, col="pink")

这些用作LogNormal分布的输入参数。我们现在按照这种偏斜分布对各种$y$$x$

xs <- ceiling(runif(datan, 0, xlim))
ys <- sapply(xs, function(x){rlnorm(1, meanlog = mu[x], sdlog = sigma[x])})

由于我们知道作为两个参数函数的LogNormal分布的均值、方差和偏度的表达式，我们可以直接确定这些。

seq <- seq(0, xlim)
mean <- sapply(seq, function(x){exp(mu[x]+sigma[x]^2/2)})
variance <- sapply(seq, function(x){(exp(sigma[x]^2)-1)*exp(2*mu[x]+sigma[x]^2)})
skewness <- sapply(seq, function(x){(exp(sigma[x]^2)+2)*sqrt(exp(sigma[x]^2)-1)})

plot(xs, ys, ylim = c(0, 4), xlab="x", ylab="y", main="Random Data with Moments")
lines(seq, mean, col="red")
lines(seq, variance, col="orange")
lines(seq, skewness, col="green")

服从LogNormal分布的情况下，找到一个检索这三个矩的过程。$Y$