R Pracma whittaker 平滑的反函数？

Question

我用的是Pracma的惠特克函数here

whittaker <- function(y, lambda = 1600, d = 2){
    #   Smoothing with a finite difference penalty
    #   y:      signal to be smoothed
    #   lambda: smoothing parameter (rough 50..1e4 smooth)
    #   d:      order of differences in penalty (generally 2)

    m <- length(y)
    E <- eye(m)
    D <- diff(E, lag = 1, differences = d)
    B <- E + (lambda * t(D) %*% D)
    z <- solve(B, y)

    return(z)
}

为此我需要找到惠特克平滑的反函数，是否存在任何逆惠特克平滑算法？即使是近似值也可能有用。

初步尝试

y = B*z            //solve(B,z)
z = B^{-1} y       // *B
y = B*z

所以我必须找出 B 即

E <- eye(length(y))
D <- diff(E, lag = 1 , differences 2)
B <- E + (lambda * t(D) %*% D)

所以

y<- B * z

据我所知，平滑是无损的（没有信息在转换中丢失）所以逆应该包含与初始数据相同的信息。

如前所述，创建惠特克平滑的反函数似乎并非不可能，但我希望 R 中已经存在这样的函数。

R中是否存在惠特克平滑的反函数？

Answer 1

看来您是完全正确的，并且与 Whittaker 平滑相反（可能在端点处除外）。假设参数 lambda 和 d 相同，恢复平滑过程确实很容易——而且你走在正确的轨道上。

library(pracma)
inv_whittaker <- function(z, lambda = 1600, d = 2) {
    m = length(z); E = eye(m)
    D = diff(E, lag = 1, differences = d)
    B = E + (lambda * t(D) %*% D)
    y = B %*% z
    return(y)
}

让我们将其应用于帮助页面上的示例，从 'jitted' 正弦函数开始。

xx = linspace(0, 10*pi, 1000)
t1 = sin(xx) + rnorm(1000)/10
t3 = whittaker(t1, lambda = 1600)

现在我们尝试从 t3 取回 t1。

t2 = inv_whittaker(t3, lambda = 1600)
max(abs(t1-t2))
## [1] 3.527845e-12

我发现这非常令人惊讶，我不确定它是否适用于所有极端情况和端点，例如当您采用 t1 精确的正弦曲线时。

R Pracma whittaker 平滑的反函数？

Inverse function of R Pracma whittaker smoothing?

math

r

smoothing