wtf0: Monotone Degree-Zero Weighted Trend Filter

In barbehenna/coleReg: COmpound Loss Emporium with Regression

Monotone Degree-Zero Weighted Trend Filter

Description

Estimate the mean of a homoscedastic sequence of independent Gaussian observations under squared error loss. The optimal separable estimator is estimated by minimizing a biased estimate of the risk under a monotone non-decreasing constraint. Under the monotonicity constraint, the biased risk estimate looks like a degree-zero weighted trend filter. The pooled adjacent violators algorithm is used to fit the estimator. The default value of 'bw' is the optimal asymptotic rate (up to a constant).

Usage

wtf0(x, s, bw = s * length(x)^(-1/6), w = "approx", knots = NULL, ...)

Arguments

x	Gaussian sequence
s	standard deviation
bw	scalar bandwidth for Gaussian kernel density estimate
w	sequence of weights (see details for specification)
knots	locations of knots
...	additional parameters passed to stats::density()

Details

When 'w' is "approx" the FFT transform is used to quickly generate weights along 'gr' in O(N log N). When 'w' is "exact" a kernel density estimate is used to generate weights along 'gr' in O(N^2). Other weights can be specified by passing a length 'length(x)-1' numeric vector.

If 'knots = NULL' the knots will be placed at the mid-point of each consecutive value of 'x': 'sort(x)'. This is an approximation, but one that works well in general and saves a lot of compute time. The correct knots are placed at the minimum of the density between each consecutive observation. Users can specify their own grid of knots by providing a length 'length(x)-1' numeric vector for 'knots'.

Value

theta_hat	estimated values of means of Gaussian sequence
x	original Gaussian sequence
s	known standard deviation
bw	bandwidth used for weights or NULL if user specified weights
w	sequence of weights used
knots	location of knots
risk.est	estimated risk of the fitted estimator
...	additional parameters passed to stats::density()

Examples

# basic usage
theta = rnorm(250)
x = theta + rnorm(250)
res = wtf0(x, s = 1)
mean((theta - res$theta_hat)^2)

# fit visualization
plot(x, x, main = "WTF0") # observed
points(theta ~ x, col = "green") # true
points(theta_hat ~ x, as.data.frame(res[1:2]), col = "red") # estimated

# alternate usage
## w exact vs approx (very close)
all.equal(
  wtf0(x, s = 1, w = "exact"),
  wtf0(x, s = 1, w = "approx")
)

## Left knots
res2 = wtf0(x, s = 1, bw = 0.2, knots = sort(x)[-250]+min(diff(sort(x)))/100)
mean((theta - res2$theta_hat)^2)

## empirical weights
res3 = wtf0(x, s = 1, w = rep(1/250, 249))
mean((theta - res3$theta_hat)^2)

## extra parameters (better approximation and different bandwidth)
## bw can only be a string when `w = "approx"` (see stats::density() for more
## details).
res4 = wtf0(x, s = 1, w = "approx", n = 1048, bw = "SJ")
mean((theta - res4$theta_hat)^2)

barbehenna/coleReg documentation built on May 8, 2022, 12:05 a.m.