机器算法验证 - 核密度估计和边界偏差 - 吾爱随笔录

核密度估计和边界偏差

机器算法验证内核平滑密度估计偏差校正

2022-03-02 19:17:17

使用哪种核密度估计器来避免边界偏差？

考虑估计密度的任务 $f_0(x)$ 有界支持，并且当接近边界时，概率质量不会减少或变为零。为了简化问题，假设密度的界限是已知的。

为了集中思想，以均匀分布为例：

给定样本量 $N$ 独立同住抽签 $\mathcal U(0,1)$ 可以考虑应用核密度估计器

\hat{f} (y) = \frac{1}{n s} \sum_{i} K (\frac{x_{i} - y}{s})

$\hat f(y) = \frac{1}{ns}\sum_i K\left( \frac{x_i-y}{s} \right)$

具有正常内核和一些平滑参数 $s$ . 为了说明边界偏差，请考虑（在软件 R: A Language and Environment for Statistical Computing 中实现）：

N <- 10000
x <- runif(N)
s <- .045

M <- 100
y <- seq(0,1,length.out=M)
out <- rep(0,M)
for (i in 1:M)
    {
        weights <- dnorm((x-y[i])/s)
        out[i] <- mean(weights)/s
    }
plot(y,out,type="l",ylim=c(0,1.5))

生成以下图

显然，该方法在捕捉密度函数的真实值时存在问题 $f_0(x)$ 在 $x$ 靠近边界。

对数样条方法效果更好，但肯定不是没有一些边界偏差

library(logspline)
set.seed(1)
N <- 10000
x <- runif(N)
m <- logspline(x,lbound=0,ubound=1,knots=seq(0,1,length.out=21))
plot(m)

2个回答

如果您知道边界，那么 Silverman 伟大的小书（统计和数据分析的密度估计）中提到的一种方法是“反射技术”。一个简单地反映有关边界（或边界）的数据。（这就是@NickCox 在他的评论中提到的。）

# Generate numbers from a uniform distribution
  set.seed(12345)
  N <- 10000
  x <- runif(N)

# Reflect the data at the two boundaries
  xReflected <- c(-x, x, 2-x)

# Construct density estimate
  d <- density(xReflected, from=0, to=1)
  plot(d$x, 3*d$y, ylab="Probability density", xlab="x", ylim=c(0,1.1), las=1)

density请注意，在这种情况下，我们最终得到的数据点数量是 3 倍，因此我们需要将函数得出的密度乘以 3 。

下面是 100 个模拟的动画显示（如上），但具有真实密度和两个估计密度（一个来自原始数据，一个来自反射数据）。density仅使用原始数据时，边界附近存在偏差是非常清楚的。

我不知道它是否有趣（鉴于最初的问题和它已经收到的答案）但是，我想建议一种替代方法。将来它也可能对某人有用（至少我希望如此）:-)。

如果您担心密度平滑方法的边界效应，我建议使用 P 样条（参见 Eilers 和马克思，1991 - 作者在第 8 条中专门讨论了密度平滑中的边界偏差）。引用 Eilers 和马克思的话，

P样条密度平滑器不受边界效应的影响，例如内核平滑器。

通常，P 样条结合了 B 样条和有限差分惩罚。密度平滑问题是 GLM 的一个特例。所以我们只需要相应地参数化我们的平滑问题。

为了回答最初的问题，我将考虑以直方图方式分组的数据。我会用 $y_{i}$ 落在 bin/bar 中的观测值的计数（但推理也可以适应密度情况） $u_{i}$ . 为了平滑这些数据，我将使用以下成分：

更平滑：惠特克更平滑（P 样条的特殊情况，基是单位矩阵）
一阶差分惩罚
IWLS 算法最大化我的惩罚可能性（参考中的 eq 36）
$L = \sum_{i} y_{i} \log μ_{i} - \sum_{i} μ_{i} - λ \sum_{i} (Δ^{(1)} η_{i})^{2}$ $L = \sum_{i} y_{i} \log \mu_{i} - \sum_{i} \mu_{i} - \lambda \sum_{i} (\Delta^{(1)} \eta_{i})^{2}$ 和 $\mu_{i} = \exp(\eta_{i})$ .

结果由下面的代码生成，固定值为 $\lambda$ （我留下了一些评论，希望它更容易阅读）。正如您将在结果中注意到的那样， $\lambda$ 参数调节最终估计的平滑度。对于一个非常高 $\lambda$ 我们得到了一条相当平坦的线。

library(colorout)

# Simulate data
set.seed(1)
N = 10000
x = runif(N)

# Construct histograms
his = hist(x, breaks = 50, plot = F)
X = his$counts
u = his$mids

# Prepare basis (I-mat) and penalty (1st difference)
B = diag(length(X))
D1 = diff(B, diff = 1)
lambda = 1e6 # fixed but can be selected (e.g. AIC)
P = lambda * t(D1) %*% D1

# Smooth
tol = 1e-8
eta = log(X + 1)
for (it in 1:20) 
{
    mu = exp(eta)
    z = X - mu + mu * eta
    a = solve(t(B) %*% (c(mu) * B) + P, t(B) %*% z)
    etnew = B %*% a
    de = max(abs(etnew - eta))
    cat('Crit', it, de, '\n')
    if(de < tol) break
    eta = etnew
}

# Plot
plot(u, exp(eta), ylim = c(0, max(X)), type = 'l', col = 2)
lines(u, X, type = 'h')

最后，我希望我的建议足够清楚，并（至少部分地）回答原始问题。

其它你可能感兴趣的问题

上一篇为什么是变分推理而不是最大似然？下一篇时间序列的简单异常值检测