Обратная функция сглаживания R Pracma Whittaker? - PullRequest
Новая
       12

Обратная функция сглаживания R Pracma Whittaker?

0 голосов
/ 12 сентября 2018

Я использую функцию Уиттекера Pracma здесь 

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?

1 Ответ

0 голосов
/ 19 сентября 2018

Похоже, вы были абсолютно правы, и есть обратное сглаживание Уиттекера (кроме, может быть, в конечных точках). Действительно, довольно просто отменить процедуру сглаживания, если предположить, что параметры 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)
}

Давайте применим это к примеру на странице справки, начиная с синусоидальной функции. 

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

Теперь мы пытаемся вернуть t1 из t3. 

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

Я нашел это довольно удивительным, и я не совсем уверен, что это работает так хорошо для всех угловых случаев и конечных точек, например, когда вы берете за t1 точную кривую синуса.

...